Commit 98e487e9 authored by Andrey Filippov's avatar Andrey Filippov

Started manual pair matching by affine transform

parent cfaee44d
...@@ -109,6 +109,7 @@ import com.elphel.imagej.ims.EventLogger; ...@@ -109,6 +109,7 @@ import com.elphel.imagej.ims.EventLogger;
import com.elphel.imagej.ims.Imx5; import com.elphel.imagej.ims.Imx5;
import com.elphel.imagej.jp4.JP46_Reader_camera; import com.elphel.imagej.jp4.JP46_Reader_camera;
import com.elphel.imagej.lwir.LwirReader; import com.elphel.imagej.lwir.LwirReader;
import com.elphel.imagej.orthomosaic.ComboMatch;
import com.elphel.imagej.readers.ChangeImageResolution; import com.elphel.imagej.readers.ChangeImageResolution;
import com.elphel.imagej.readers.DumpImageMetadata; import com.elphel.imagej.readers.DumpImageMetadata;
import com.elphel.imagej.readers.EyesisTiff; import com.elphel.imagej.readers.EyesisTiff;
...@@ -841,6 +842,8 @@ public class Eyesis_Correction implements PlugIn, ActionListener { ...@@ -841,6 +842,8 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
addButton("Build World", panelLWIRWorld, color_process); addButton("Build World", panelLWIRWorld, color_process);
addButton("Test IMX5", panelLWIRWorld, color_process); addButton("Test IMX5", panelLWIRWorld, color_process);
addButton("Show mice", panelLWIRWorld, color_process); addButton("Show mice", panelLWIRWorld, color_process);
addButton("Set pair", panelLWIRWorld, color_process);
addButton("Warp pair", panelLWIRWorld, color_process);
plugInFrame.add(panelLWIRWorld); plugInFrame.add(panelLWIRWorld);
} }
...@@ -5672,6 +5675,14 @@ public class Eyesis_Correction implements PlugIn, ActionListener { ...@@ -5672,6 +5675,14 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
ImagePlus imp_src = WindowManager.getCurrentImage(); ImagePlus imp_src = WindowManager.getCurrentImage();
ImagePlus imp_mice = subtractAverage(imp_src); ImagePlus imp_mice = subtractAverage(imp_src);
imp_mice.show(); imp_mice.show();
} else if (label.equals("Set pair")) {
DEBUG_LEVEL = MASTER_DEBUG_LEVEL;
EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
ComboMatch.openTestPair();
} else if (label.equals("Warp pair")) {
DEBUG_LEVEL = MASTER_DEBUG_LEVEL;
EYESIS_CORRECTIONS.setDebug(DEBUG_LEVEL);
ComboMatch.testPair();
} }
} }
......
package com.elphel.imagej.orthomosaic;
import java.util.Arrays;
import java.util.concurrent.atomic.AtomicInteger;
import com.elphel.imagej.common.GenericJTabbedDialog;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.tileprocessor.ImageDtt;
import Jama.Matrix;
import ij.ImagePlus;
import ij.gui.PointRoi;
import ij.process.FloatPolygon;
public class ComboMatch {
public static ImagePlus imp_src1 = null;
public static ImagePlus imp_src2 = null;
public static boolean openTestPair() {
String [] image_paths = {
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/linked/linked_compare/02/1697877410_420287-RECT-PIX0.01-FLAT_CLN-GEO0.tiff",
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/linked/linked_compare/02/1697877412_004148-RECT-PIX0.01-FLAT_CLN-GEO0.tiff"
};
double [][][] image_points = {
{ { 487.0, 637.0},
{2237.5, 932.5},
{2004.5, 1929.5},
{ 320.5, 1951.500}
},{
{ 629.0, 179.0},
{2376.5, 479.5},
{2139.833, 1474.5},
{ 461.5, 1493.5}}};
boolean set_coords = true; // false;
GenericJTabbedDialog gd = new GenericJTabbedDialog("Set image pair",1090,900);
gd.addStringField ("First image path", image_paths[0], 120, "First image full path");
gd.addStringField ("Second image path", image_paths[1], 120, "Second image full path");
gd.addCheckbox ("Set image markers to following coordinates", set_coords, "Delete markers and set them according to data.");
for (int n = 0; n < image_points.length; n++) {
for (int i = 0; i < image_points[n].length; i++) {
gd.addMessage("image["+n+"] point "+i+":");
gd.addNumericField("x =", image_points[n][i][0], 3,7,"pix",
"Marker "+i+ " for image "+n+" X-coordinate");
gd.addNumericField("y =", image_points[n][i][1], 3,7,"pix",
"Marker "+i+ " for image "+n+" Y-coordinate");
}
}
gd.showDialog();
if (gd.wasCanceled()) return false;
image_paths[0] = gd.getNextString();
image_paths[1] = gd.getNextString();
set_coords = gd.getNextBoolean();
for (int n = 0; n < image_points.length; n++) {
for (int i = 0; i < image_points[n].length; i++) {
image_points[n][i][0] = gd.getNextNumber();
image_points[n][i][1] = gd.getNextNumber();
}
}
imp_src1 = new ImagePlus(image_paths[0]);
imp_src2 = new ImagePlus(image_paths[1]);
if (set_coords) {
PointRoi [] rois = new PointRoi[image_points.length];
for (int n = 0; n < image_points.length; n++) {
rois[n] = new PointRoi();
rois[n].setOptions("label");
for (int i = 0; i < image_points[n].length; i++) {
rois[n].addPoint(image_points[n][i][0],image_points[n][i][1],i+1);
}
}
imp_src1.setRoi(rois[0]);
imp_src2.setRoi(rois[1]);
}
imp_src1.show();
imp_src2.show();
if (!(imp_src1.getRoi() instanceof PointRoi) || !(imp_src2.getRoi() instanceof PointRoi)) {
System.out.println("At least one iomage does not have PointRoi");
return false;
}
FloatPolygon [] fp = {((PointRoi) imp_src1.getRoi()).getContainedFloatPoints(),
((PointRoi) imp_src2.getRoi()).getContainedFloatPoints()};
int expected_points = image_points[0].length;
for (int n = 0; n < fp.length; n++) {
if (fp[n].xpoints.length != expected_points) {
System.out.println("openTestPair(): expecting "+expected_points+" points, image "+n+"has "+fp[n].xpoints.length);
return false;
}
}
double [][][] points_xy= new double [fp.length][fp[0].xpoints.length][2];
for (int n = 0; n < points_xy.length; n++) {
for (int i = 0; i < points_xy[n].length; i++) {
points_xy[n][i][0] = fp[n].xpoints[i];
points_xy[n][i][1] = fp[n].ypoints[i];
}
}
System.out.println(String.format("%4s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s",
"N","X0","Y0","X1","Y1", "length1", "length2", "ratio"));
for (int i = 0; i < points_xy[0].length; i++) {
int ip = (i - 1 + points_xy[0].length) % points_xy[0].length;
double [] ll = new double[2];
for (int n = 0; n < points_xy.length; n++) {
double dx = points_xy[n][i][0]-points_xy[n][ip][0];
double dy = points_xy[n][i][1]-points_xy[n][ip][1];
ll[n] = Math.sqrt(dx*dx+dy*dy);
}
System.out.println(String.format("%4d\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t%8.6f",
i, points_xy[0][i][0], points_xy[0][i][1],points_xy[1][i][0], points_xy[1][i][1],
ll[0], ll[1], ll[1]/ll[0]));
}
return true;
}
public static boolean testPair() {
if (!(imp_src1.getRoi() instanceof PointRoi) || !(imp_src2.getRoi() instanceof PointRoi)) {
System.out.println("At least one iomage does not have PointRoi");
return false;
}
FloatPolygon [] fp = {((PointRoi) imp_src1.getRoi()).getContainedFloatPoints(),
((PointRoi) imp_src2.getRoi()).getContainedFloatPoints()};
int expected_points = 4; // image_points[0].length;
for (int n = 0; n < fp.length; n++) {
if (fp[n].xpoints.length != expected_points) {
System.out.println("openTestPair(): expecting "+expected_points+" points, image "+n+"has "+fp[n].xpoints.length);
return false;
}
}
double [][][] points_xy= new double [fp.length][fp[0].xpoints.length][2];
for (int n = 0; n < points_xy.length; n++) {
for (int i = 0; i < points_xy[n].length; i++) {
points_xy[n][i][0] = fp[n].xpoints[i];
points_xy[n][i][1] = fp[n].ypoints[i];
}
}
System.out.println(String.format("%4s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s",
"N","X0","Y0","X1","Y1", "length1", "length2", "ratio"));
for (int i = 0; i < points_xy[0].length; i++) {
int ip = (i - 1 + points_xy[0].length) % points_xy[0].length;
double [] ll = new double[2];
for (int n = 0; n < points_xy.length; n++) {
double dx = points_xy[n][i][0]-points_xy[n][ip][0];
double dy = points_xy[n][i][1]-points_xy[n][ip][1];
ll[n] = Math.sqrt(dx*dx+dy*dy);
}
System.out.println(String.format("%4d\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t%8.6f",
i, points_xy[0][i][0], points_xy[0][i][1],points_xy[1][i][0], points_xy[1][i][1],
ll[0], ll[1], ll[1]/ll[0]));
}
double [][] ab = getAffineMatrix(
points_xy[1], // double [][] x,
points_xy[0]); // double [][] y){
double [][] a = {{ab[0][0],ab[0][1]},{ab[1][0],ab[1][1]}};
double [] b = {ab[2][0],ab[2][1]};
System.out.println("Matrix A:");
System.out.println(String.format("%10.5f\t%10.5f", a[0][0], a[0][1]));
System.out.println(String.format("%10.5f\t%10.5f", a[1][0], a[1][1]));
System.out.println("Matrix B:");
System.out.println(String.format("%10.5f", b[0]));
System.out.println(String.format("%10.5f", b[1]));
final int [] widths = {imp_src1.getWidth(),imp_src2.getWidth()};
final int [] heights = {imp_src1.getHeight(),imp_src2.getHeight()};
int [][][] corners = new int[2][4][2];
for (int n =0; n < corners.length; n++) {
corners[n][0][0] = 0;
corners[n][0][1] = 0;
corners[n][1][0] = widths[n]-1;
corners[n][1][1] = 0;
corners[n][2][0] = widths[n]-1;
corners[n][2][1] = heights[n]-1;
corners[n][3][0] = 0;
corners[n][3][1] = heights[n]-1;
}
double [][] fcorners = new double[corners[1].length][2];
int [][] xy_minmax = {{0, widths[0]-1},{0, heights[0]-1}};
int [][] xy_minmax1 = new int[2][2];
for (int i = 0; i < corners[1].length; i++) {
fcorners[i][0] = corners[1][i][0] * a[0][0]+ corners[1][i][1] * a[0][1] + b[0];
fcorners[i][1] = corners[1][i][0] * a[1][0]+ corners[1][i][1] * a[1][1] + b[1];
xy_minmax[0][0] = Math.min(xy_minmax[0][0], (int) Math.floor(fcorners[i][0]));
xy_minmax[0][1] = Math.max(xy_minmax[0][1], (int) Math.ceil(fcorners[i][0]));
xy_minmax[1][0] = Math.min(xy_minmax[1][0], (int) Math.floor(fcorners[i][1]));
xy_minmax[1][1] = Math.max(xy_minmax[1][1], (int) Math.ceil(fcorners[i][1]));
if (i == 0) {
xy_minmax1[0][0] = (int) Math.floor(fcorners[i][0]);
xy_minmax1[0][1] = (int) Math.ceil(fcorners[i][0]);
xy_minmax1[1][0] = (int) Math.floor(fcorners[i][1]);
xy_minmax1[1][1] = (int) Math.ceil(fcorners[i][1]);
} else {
xy_minmax1[0][0] = Math.min(xy_minmax1[0][0], (int) Math.floor(fcorners[i][0]));
xy_minmax1[0][1] = Math.max(xy_minmax1[0][1], (int) Math.ceil(fcorners[i][0]));
xy_minmax1[1][0] = Math.min(xy_minmax1[1][0], (int) Math.floor(fcorners[i][1]));
xy_minmax1[1][1] = Math.max(xy_minmax1[1][1], (int) Math.ceil(fcorners[i][1]));
}
}
final int x0 = Math.max(0, -xy_minmax[0][0]);
final int y0 = Math.max(0, -xy_minmax[1][0]);
final int width = x0 + xy_minmax[0][1]+1;
final int height = y0 + xy_minmax[1][1]+1;
double [][] opix = new double [2][width*height];
for (int n = 0; n < opix.length; n++) {
Arrays.fill(opix[n], Double.NaN);
}
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final int width0= widths[0];
final int npix0 = width0*heights[0];
final float [] pixels0 = (float []) imp_src1.getProcessor().getPixels();
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < npix0; ipix = ai.getAndIncrement()) {
int y = ipix / width0;
int x = ipix % width0;
int p = x0 + x + (y0 + y)* width;
opix[0][p] = pixels0[ipix];
}
}
};
}
ImageDtt.startAndJoin(threads);
final int x10 = x0+ xy_minmax1[0][0]; // output image coordiantes
final int y10 = y0+ xy_minmax1[1][0]; // output image coordiantes
final int width1 = x0 + xy_minmax1[0][1] + 1;
final int height1 = y0 + xy_minmax1[0][1] + 1;
final int npix1 = width1 * height1;
double [][] ab1 = getAffineMatrix(
points_xy[0], // double [][] x,
points_xy[1]); // double [][] y){
final double [][] a1 = {{ab1[0][0],ab1[0][1]},{ab1[1][0],ab1[1][1]}};
final double [] b1 = {ab1[2][0]-x0,ab1[2][1]-y0}; // from output coordinates to image2 coordinates
System.out.println("Matrix A1:");
System.out.println(String.format("%10.5f\t%10.5f", a1[0][0], a1[0][1]));
System.out.println(String.format("%10.5f\t%10.5f", a1[1][0], a1[1][1]));
System.out.println("Matrix B1:");
System.out.println(String.format("%10.5f", b1[0]));
System.out.println(String.format("%10.5f", b1[1]));
ai.set(0);
final int img2_width = widths[1];
final int img2_height = heights[1];
final float [] pixels1 = (float []) imp_src2.getProcessor().getPixels();
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < npix1; ipix = ai.getAndIncrement()) {
int y = ipix / width1 + x10; // output image coordinates
int x = ipix % width1 + y10; // output image coordinates
int p = x + y * width;
double fx = a1[0][0] * x + a1[0][1] * y + b1[0];
double fy = a1[1][0] * x + a1[1][1] * y + b1[1];
int ix = (int) Math.floor(fx);
int iy = (int) Math.floor(fy);
double dx = fx - ix;
double dy = fy - iy;
if ((ix >= 0) && (iy >= 0) && (ix < (img2_width -1)) && (iy < (img2_height -1))) {
int indx0 = ix + iy * img2_width;
double d00 = pixels1[indx0];
double d01 = pixels1[indx0+1];
double d10 = pixels1[indx0+img2_width];
double d11 = pixels1[indx0+img2_width+1];
opix[1][p] = d00*(1-dx)*(1-dy)+d01*dx*(1-dy)+d10*(1-dx)*dy+d11*dx*dy;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
ShowDoubleFloatArrays.showArrays(
opix,
width,
height,
true,
"overlap");
return true;
}
public static double [][] getAffineMatrix(
double [][] x,
double [][] y){
if (x.length != y.length) {
System.out.println("getAffineMatrix(): x.length != y.length ("+x.length+" != "+y.length+")");
return null;
}
/*
Ui = A00 * Xi + A01* Yi + B0
Vi = A10 * Xi + A11* Yi + B1
A00^2*Xi^2 + A01^2*Yi^2 +B0^2 +Ui^2 + 2*A00*Xi*A01*Yi + 2*A00*Xi*B0 - 2*A00*Xi*Ui + 2*A01*Yi*B0 -2*A01*Yi*Ui - 2*B0*Ui
d/2dA00 = A00*Xi^2 + A01*Xi*Yi + B0 *Xi - XiUi = 0
d/2dA01 = A01*Yi^2 + A00*Xi*Yi + B0 *Yi - YiUi = 0
d/2dB0 = B0 + A00*Xi + A01*Yi - Ui = 0
d/2dA10 = A10*Xi^2 + A11*Xi*Yi + B1 *Xi - XiVi = 0
d/2dA11 = A11*Yi^2 + A10*Xi*Yi + B1 *Yi - YiVi = 0
d/2dB1 = B1 + A10*Xi + A11*Yi - Vi = 0
SX2*A00 + SXY*A01 + 0*A10 + 0*A11 + SX*B0 + 0*B1 = SXU
SXY*A00 + SY2*A01 + 0*A10 + 0*A11 + SY*B0 + 0*B1 = SYU
SX *A00 + SY *A01 + 0*A10 + 0*A11 + 1*B0 + 0*B1 = SU
0 *A00 + 0 *A01 + SX2*A10 + SXY*A11 + 0*B0 + SX*B1 = SXV
0 *A00 + 0 *A01 + SXY*A10 + SY2*A11 + 0*B0 + SY*B1 = SYV
0 *A00 + 0 *A01 + SX *A10 + SY *A11 + 0*B0 + 1*B1 = SV
*/
double sx=0, sy=0, sx2=0, sy2=0, sxy=0,sxu=0,syu=0,sxv=0,syv=0,su=0,sv=0;
for (int n = 0; n < x.length; n++) {
sx += x[n][0];
sy += x[n][1];
sx2+= x[n][0]*x[n][0];
sy2+= x[n][1]*x[n][1];
sxy+= x[n][0]*x[n][1];
sxu+= x[n][0]*y[n][0];
syu+= x[n][1]*y[n][0];
su+= y[n][0];
sxv+= x[n][0]*y[n][1];
syv+= x[n][1]*y[n][1];
sv+= y[n][1];
}
double s0 = x.length;
double [][] a = {
{sx2, sxy, 0, 0, sx, 0},
{sxy, sy2, 0, 0, sy, 0},
{ sx, sy, 0, 0, s0, 0},
{ 0, 0, sx2, sxy, 0, sx},
{ 0, 0, sxy, sy2, 0, sy},
{ 0, 0, sx, sy, 0, s0}};
Matrix A = new Matrix(a);
A.print(15,4);
double [][] b = {{sxu,syu,su,sxv,syv,sv}};
Matrix B = new Matrix(b).transpose();
B.print(15,4);
double [] v = A.solve(B).getRowPackedCopy();
double [][] ab = {
{v[0], v[1]}, // |a00, a01|,
{v[2], v[3]}, // |a10, a11|,
{v[4], v[5]}}; // | b0, b1|,
return ab;
}
}
...@@ -11624,6 +11624,10 @@ public class ImageDttCPU { ...@@ -11624,6 +11624,10 @@ public class ImageDttCPU {
* From Stephan Preibisch's Multithreading.java class. See: * From Stephan Preibisch's Multithreading.java class. See:
* http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD * http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
*/ */
public static Thread[] newThreadArray() {
return newThreadArray(THREADS_MAX);
}
public static Thread[] newThreadArray(int maxCPUs) { // USED in lwir public static Thread[] newThreadArray(int maxCPUs) { // USED in lwir
int n_cpus = Runtime.getRuntime().availableProcessors(); int n_cpus = Runtime.getRuntime().availableProcessors();
if (n_cpus>maxCPUs)n_cpus=maxCPUs; if (n_cpus>maxCPUs)n_cpus=maxCPUs;
......
...@@ -2673,12 +2673,14 @@ public class TexturedModel { ...@@ -2673,12 +2673,14 @@ public class TexturedModel {
} }
} }
// ------- // -------
boolean mixed_flat = false;
if (gsmth_enable && (!has_sfm || gsmth_sfm_gnd)) { if (gsmth_enable && (!has_sfm || gsmth_sfm_gnd)) {
double [] smooth_ground = scenes[ref_index].getSmoothGround( double [] smooth_ground = scenes[ref_index].getSmoothGround(
clt_parameters, clt_parameters,
has_sfm, // boolean sfm_mode, has_sfm, // boolean sfm_mode,
debugLevel); debugLevel);
if (has_sfm && gsmth_sfm_deviate) { if (has_sfm && gsmth_sfm_deviate) {
mixed_flat = true;
// Change! // Change!
// Use FG for FG, and ground plane as BG? // Use FG for FG, and ground plane as BG?
TileNeibs tn = new TileNeibs(tilesX, smooth_ground.length/tilesX); TileNeibs tn = new TileNeibs(tilesX, smooth_ground.length/tilesX);
...@@ -3189,7 +3191,7 @@ public class TexturedModel { ...@@ -3189,7 +3191,7 @@ public class TexturedModel {
String suffix ="-RECT"; String suffix ="-RECT";
suffix +="-PIX"+pix_size * hdr_whs[2]; suffix +="-PIX"+pix_size * hdr_whs[2];
if (gsmth_enable) { if (gsmth_enable) {
suffix +="-FLAT"; // flattened ground suffix +=mixed_flat ? "-FLAT_MIX":"-FLAT_CLN"; // flattened ground - mixed or clean
} }
scenes[ref_index].saveDoubleArrayInModelDirectory( scenes[ref_index].saveDoubleArrayInModelDirectory(
suffix+"-FULL", // String suffix, suffix+"-FULL", // String suffix,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment