Skip to content

I
imagej-elphel
Commit c9ed0a3d authored by Andrey Filippov's avatar Andrey Filippov
Browse files
Options

trying filtering by SIFT-like features for better averaging

parent 1517acca
Show whitespace changes
Inline Side-by-side
Showing with 1155 additions and 6 deletions
src/main/java/com/elphel/imagej/correction/Eyesis_Correction.java
View file @ c9ed0a3d
......@@ -111,6 +111,7 @@ import com.elphel.imagej.jp4.JP46_Reader_camera;
import com.elphel.imagej.lwir.LwirReader;
import com.elphel.imagej.orthomosaic.OrthoMap;
import com.elphel.imagej.orthomosaic.ComboMatch;
import com.elphel.imagej.orthomosaic.OrangeTest;
import com.elphel.imagej.readers.ChangeImageResolution;
import com.elphel.imagej.readers.DumpImageMetadata;
import com.elphel.imagej.readers.ElphelTiffReader;
......@@ -165,7 +166,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
// private static final long serialVersionUID = -1507307664341265263L;
private Panel panel1, panel2, panel3, panel4, panel5, panel5a, panel6, panel7, panelPostProcessing1,
panelPostProcessing2, panelPostProcessing3, panelDct1, panelClt1, panelClt2, panelClt3, panelClt4,
panelClt5, panelClt5aux, panelClt_GPU, panelLWIR, panelLWIR16, panelLWIRWorld;
panelClt5, panelClt5aux, panelClt_GPU, panelLWIR, panelLWIR16, panelLWIRWorld, panelOrange;
JP46_Reader_camera JP4_INSTANCE = null;
// private deBayerScissors debayer_instance;
......@@ -532,7 +533,7 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
plugInFrame.addKeyListener(IJ.getInstance());
// int menuRows=4 + (ADVANCED_MODE?4:0) + (MODE_3D?3:0) + (DCT_MODE?6:0) + (GPU_MODE?1:0) +(LWIR_MODE?2:0);
int menuRows = 4 + (ADVANCED_MODE ? 4 : 0) + (MODE_3D ? 3 : 0) + (DCT_MODE ? 7 : 0) + (GPU_MODE ? 1 : 0)
+ (LWIR_MODE ? 3 : 0);
+ (LWIR_MODE ? 4 : 0);
plugInFrame.setLayout(new GridLayout(menuRows, 1));
panel6 = new Panel();
panel6.setLayout(new GridLayout(1, 0, 5, 5));
......@@ -861,6 +862,11 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
addButton("Properties compare", panelLWIRWorld, color_process);
plugInFrame.add(panelLWIRWorld);
panelOrange = new Panel();
panelOrange.setLayout(new GridLayout(1, 0, 5, 5)); // rows, columns, vgap, hgap
addButton("Test Orange", panelOrange, color_process);
addButton("Process Merged", panelOrange, color_process);
plugInFrame.add(panelOrange);
}
plugInFrame.pack();
//"LWIR batch"
......@@ -5788,7 +5794,12 @@ public class Eyesis_Correction implements PlugIn, ActionListener {
OrthoMap.testPatternCorrelate(imp_sel);
} else if (label.equals("Properties compare")) {
QuadCLTCPU.compareProperties();
} else if (label.equals("Test Orange")) {
OrangeTest.testOrange();
} else if (label.equals("Process Merged")) {
OrangeTest.processMerged();
}
//
}
//
public boolean debugInitOneScene() {
......
src/main/java/com/elphel/imagej/orthomosaic/ComboMatch.java
View file @ c9ed0a3d
......@@ -2524,10 +2524,10 @@ adjusted affines[1] for a pair: 1694564291_293695/1694564778_589341
{ 0, 0, sxy, sy2, 0, sy},
{ 0, 0, sx, sy, 0, s0}};
Matrix A = new Matrix(a);
A.print(15,4);
// A.print(15,4);
double [][] b = {{sxu,syu,su,sxv,syv,sv}};
Matrix B = new Matrix(b).transpose();
B.print(15,4);
// B.print(15,4);
double [] v = A.solve(B).getRowPackedCopy();
double [][] ab = {
{v[0], v[1]}, // |a00, a01|,
......@@ -2535,6 +2535,7 @@ adjusted affines[1] for a pair: 1694564291_293695/1694564778_589341
{v[4], v[5]}}; // | b0, b1|,
return ab;
}
public static boolean testReadTiff() {
String image_paths="";
......
src/main/java/com/elphel/imagej/orthomosaic/OrangeTest.java 0 → 100644
View file @ c9ed0a3d
package com.elphel.imagej.orthomosaic;
import java.awt.Rectangle;
import java.util.Arrays;
import java.util.concurrent.atomic.AtomicInteger;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.TileProcessor;
import ij.ImagePlus;
import ij.ImageStack;
public class OrangeTest {
public static double [][][] matching_points_initial= {
{ // 1697877488_162849 - 117 sfm 37.1
{ 110.500, 452.500},
{ 287.500, 121.500},
{1294.833, 545.833},
{1307.500, 162.167}
},
{ // 1697877490_930437 - 42 sfm 6.4
{ 410.167, 979.833},
{ 585.500, 649.500},
{1579.167,1126.833},
{1604.750, 737.000}
},
{ // 1697877491_613999 - 76 sfm 12.2 0, 76
{ 423.500, 1260.250},
{ 598.833, 937.500},
{1587.000, 1404.000},
{1613.167, 1017.167}
},
{ // 1697877491_997460 - 85 sfm 14.2 master
{ 464.000, 1551.000},
{ 617.167, 1222.500},
{1630.167, 1624.500},
{1630.250, 1239.500}
},
{ // 1697877492_314232 - 43 sfm 5.1
{ 637.000, 1830.000},
{ 763.167, 1489.167},
{1807.167, 1807.167},
{1772.833, 1420.500}
},
{ // 1697877492_614332 - 38 sfm 4.2
{ 731.000, 1952.500},
{ 828.500, 1606.500},
{1893.000, 1829.500},
{1826.500, 1453.000}
}
};
public static double [][][] matching_points_first_corr= { // /media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/02-first_point_correction/
{ // 1697877488_162849 - 117 sfm 37.1
{ 105.500, 452.500},
{ 289.500, 121.500},
{1298.833, 550.833},
{1305.500, 164.167}
},
{ // 1697877490_930437 - 42 sfm 6.4
{ 410.167, 981.833},
{ 585.500, 651.500},
{1579.167,1128.833},
{1606.750, 741.000}
},
{ // 1697877491_613999 - 76 sfm 12.2 0, 76
{ 419.500, 1260.250},
{ 596.833, 937.500},
{1585.000, 1404.000},
{1613.167, 1019.167}
},
{ // 1697877491_997460 - 85 sfm 14.2 MASTER
{ 464.000, 1551.000},
{ 617.167, 1222.500},
{1630.167, 1624.500},
{1630.250, 1239.500}
},
{ // 1697877492_314232 - 43 sfm 5.1
{ 637.000, 1828.000},
{ 765.167, 1489.167},
{1805.167, 1805.167},
{1770.833, 1422.500}
},
{ // 1697877492_614332 - 38 sfm 4.2
{ 733.000, 1950.500},
{ 832.500, 1604.500},
{1891.000, 1829.500},
{1826.500, 1455.000}
}
};
public static double [][][] matching_points_second_corr= { // /media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/02-first_point_correction/
{ // 1697877488_162849 - 117 sfm 37.1
{ 106.500, 453.500}, //+1,+1
{ 290.500, 123.500}, //+1,+2
{1299.833, 551.833}, //+1,+1
{1305.500, 164.667} // 0,+.5
},
{ // 1697877490_930437 - 42 sfm 6.4
{ 411.167, 981.833}, //+1,0
{ 585.500, 652.500}, // 0,+1
{1579.167,1128.833}, // 0, 0
{1606.750, 741.500} // 0,+.5
},
{ // 1697877491_613999 - 76 sfm 12.2 0, 76
{ 420.500, 1260.250}, //+1, 0
{ 596.833, 937.500}, // 0, 0
{1586.000, 1404.000}, //+1, 0
{1613.167, 1019.667} // 0,+.5
},
{ // 1697877491_997460 - 85 sfm 14.2 MASTER
{ 464.000, 1551.000},
{ 617.167, 1222.500},
{1630.167, 1624.500},
{1630.250, 1239.500}
},
{ // 1697877492_314232 - 43 sfm 5.1
{ 637.000, 1828.000}, // 0, 0
{ 765.167, 1489.167}, // 0, 0
{1805.167, 1805.167}, //+1,+1
{1770.833, 1422.500} // 0,+1
},
{ // 1697877492_614332 - 38 sfm 4.2
{ 732.000, 1950.500}, //-1, 0
{ 831.500, 1604.500}, //-1, 0
{1892.000, 1830.500}, //+1,+1
{1826.000, 1456.000} //-.5,+1
}
};
public static String [] src_files = {
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877488_162849-TERRAIN-MERGED.tiff",
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877490_930437-TERRAIN-MERGED.tiff",
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877491_613999-TERRAIN-MERGED.tiff",
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877491_997460-TERRAIN-MERGED.tiff",
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877492_314232-TERRAIN-MERGED.tiff",
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877492_614332-TERRAIN-MERGED.tiff"};
public static double [] value_offsets = {33.0, 0.0, 0.0, 0.0, 0.0, 0.0};
public static String merged_file = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/04-third-corr/merged_seq_scale2.0-slices260_x92_y280_w320_h120_zeronan.tiff";
public static double points_scale = 4.0;
public static double [][][] matching_points= { // /media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/02-first_point_correction/
{ // 1697877488_162849 - 117 sfm 37.1
{ 105.500, 453.500}, //-1, 0
{ 292.500, 123.500}, //+2, 0
{1299.833, 552.333}, // 0, +.5
{1304.500, 164.667} //-1, 0
},
{ // 1697877490_930437 - 42 sfm 6.4
{ 411.167, 982.833}, // 0,+1
{ 585.500, 653.500}, // 0,+1
{1578.167,1128.833}, //-1, 0
{1606.750, 741.500} // 0, 0
},
{ // 1697877491_613999 - 76 sfm 12.2 0, 76
{ 421.000, 1260.250}, //+.5,0
{ 596.833, 937.500}, // 0, 0
{1586.000, 1404.000}, // 0, 0
{1613.167, 1020.167} // 0,+.5
},
{ // 1697877491_997460 - 85 sfm 14.2 MASTER
{ 464.000, 1551.000},
{ 617.167, 1222.500},
{1630.167, 1624.500},
{1630.250, 1239.500}
},
{ // 1697877492_314232 - 43 sfm 5.1
{ 637.000, 1828.000}, // 0, 0
{ 764.667, 1489.167}, //-.5,0
{1806.167, 1806.167}, //+1,+1
{1770.833, 1422.000} //0,-.5
},
{ // 1697877492_614332 - 38 sfm 4.2
{ 732.000, 1950.000}, // 0,-.5
{ 830.500, 1604.500}, //-1, 0
{1893.000, 1830.500}, //+1, 0
{1826.000, 1455.000} // 0,-1
}
};
public static int [][] full_segments = {{ 0,117},{0,42},{0,76},{ 0,85},{ 0,43},{ 0,38}};
public static int [][] match_segments = {{ 0,117},{0, 7},{0,76},{62,85},{24,43},{20,38}}; // maybe improve to cut symmetrically?
public static int [][] single_segments = {{60,61},{0,1},{0,1},{0,1},{0,1},{0,1}};
public static void processMerged() {
int num_vert = 1; //4;
int vgap = 0; // 10;
double diff_time_rt = 2.9; // 5.5;
double diff_time_rxy = 2.9; // 5.5;
double space_rt = 2.9; // 5.5;
double space_rxy = 2.9; // 5.5;
double frac_outliers = 0.1;
boolean repeat_border= true;
boolean show_orig_dt = false;
boolean show_derivs_dt = true;
ImagePlus imp_merged = new ImagePlus(merged_file);
ImageStack stack = imp_merged.getImageStack();
final int width = imp_merged.getWidth();
// final int merged_height = imp_merged.getHeight();
// String [] stack_labels = stack.getSliceLabels();
final double [][] dpixels = new double [stack.getSize()][];
for (int i = 0; i < dpixels.length; i++) {
float [] fpixels = (float[]) stack.getPixels(i+1);
dpixels[i] = new double[fpixels.length];
for (int j = 0; j < fpixels.length; j++) {
dpixels[i][j] = fpixels[j];
}
}
if (show_orig_dt) {
double [][] out_data_dt = timeToY(
dpixels, // final double [][] multi_image,
width, // final int width,
num_vert, // final int num_vert,
vgap); // final int vgap)
int out_height = out_data_dt[0].length/width;
String [] titles = new String[out_data_dt.length];
for (int i = 0; i < titles.length; i++) {
if (num_vert == 1) {
titles[i]=String.format("row %d",i);
} else {
titles[i]=String.format("rows: %d-%d",i*num_vert,(i+1)*num_vert-1);
}
}
ShowDoubleFloatArrays.showArrays(
out_data_dt,
width,
out_height,
true,
"time_section-"+out_data_dt.length,
titles);
}
double [][][] kernel_diff_time= kernelDiffTime(
diff_time_rt, //double rt,
diff_time_rxy); // double rxy)
showKernel(
kernel_diff_time, // double [][][] kernel,
"kernel_diff_time_rt"+diff_time_rt+"_rxy"+diff_time_rxy); // String title)
double [][] time_derivs = convolve3d(
dpixels, // final double [][] data_in, // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
width, // final int width,
kernel_diff_time, // final double [][][] kernel, // (odd)x(odd)x(odd}
repeat_border); // final boolean repeat_border)
double [][][][] space_kernels= new double[5][][][];
// double space_rt = 5.5;
// double space_rxy = 5.5;
for (int dir = 0; dir < space_kernels.length; dir++) {
space_kernels[dir]= kernelXY(
dir, // int dir, // 0: +y, 1: +x+y, 2:+x, 3 +x-t, 4: 0: +center will be voting 1 of 10 (each plus and minus), weighted by abs(kernelDiffTime)
space_rt, // double rt,
space_rxy); // double rxy)
showKernel(
space_kernels[dir], // double [][][] kernel,
"space_kernels-dir"+dir+"_rt"+space_rt+"_rxy"+space_rxy); // String title)
}
double [][][] space_conv = new double[5][][];
double [][][] space_conv_temp = new double[5][][];
for (int dir = 0; dir < space_kernels.length; dir++) {
space_conv[dir] = convolve3d(
dpixels, // final double [][] data_in, // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
width, // final int width,
space_kernels[dir], // final double [][][] kernel, // (odd)x(odd)x(odd}
repeat_border); // final boolean repeat_border)
}
double [] scales = normalizeSpaceKernels(
space_conv); // double [][][] conv_data) // same data convolved with different kernels, may have NaNs
double average_space = scales[space_kernels.length]; // average std of all space kernels
double [] time_averages = normalizeTimeDerivs( // return {std, average abs()}
time_derivs, // final double [][] conv_data,
frac_outliers); // final double frac_outliers)
double time_std = time_averages[0]; // [1] - average of the absolute values;
for (int dir = 0; dir < space_kernels.length; dir++) {
scaleConvolutions(
space_conv[dir], // final double [][] conv_data, // same data convolved with different kernels, may have NaNs
scales[dir]); // final double scale)
}
scaleConvolutions(
time_derivs, // final double [][] conv_data, // same data convolved with different kernels, may have NaNs
average_space/time_std); // final double scale)
int [][] best_space_dir = new int [space_conv[0].length][];
double [][] best_space_value = getBestSpace(
space_conv, // final double [][][] space_conv,
best_space_dir); // final int [][] best_space) // null or int [nimg][]
if (show_derivs_dt) {
for (int dir = 0; dir < space_kernels.length; dir++) {
space_conv_temp[dir] = timeToY(
space_conv[dir], // final double [][] multi_image,
width, // final int width,
num_vert, // final int num_vert,
vgap); // final int vgap)
}
double [][] time_derivs_temp = timeToY(
time_derivs, // final double [][] multi_image,
width, // final int width,
num_vert, // final int num_vert,
vgap); // final int vgap)
int out_height = time_derivs_temp[0].length/width;
String [] titles = new String[time_derivs_temp.length];
for (int i = 0; i < titles.length; i++) {
if (num_vert == 1) {
titles[i]=String.format("row %d",i);
} else {
titles[i]=String.format("rows: %d-%d",i*num_vert,(i+1)*num_vert-1);
}
}
ShowDoubleFloatArrays.showArrays(
time_derivs_temp,
width,
out_height,
true,
"time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
titles);
for (int dir = 0; dir < space_kernels.length; dir++) {
ShowDoubleFloatArrays.showArrays(
space_conv_temp[dir],
width,
out_height,
true,
"space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
titles);
}
double [][] best_space_value_temp = timeToY(
best_space_value, // final double [][] multi_image,
width, // final int width,
num_vert, // final int num_vert,
vgap); // final int vgap)
ShowDoubleFloatArrays.showArrays(
best_space_value_temp,
width,
out_height,
true,
"best_space_value-rt"+space_rt+"-rxy"+space_rxy,
titles);
double [][] best_space_dir_dbg = new double [best_space_dir.length][best_space_dir[0].length];
for (int nimg = 0; nimg < best_space_dir.length; nimg++) {
Arrays.fill (best_space_dir_dbg[nimg], Double.NaN);
for (int i = 0; i < best_space_dir[nimg].length; i++) if (best_space_dir[nimg][i] >= 0){
best_space_dir_dbg[nimg][i] = best_space_dir[nimg][i];
}
}
double [][] best_space_dir_temp = timeToY(
best_space_dir_dbg, // final double [][] multi_image,
width, // final int width,
num_vert, // final int num_vert,
vgap); // final int vgap)
ShowDoubleFloatArrays.showArrays(
best_space_dir_temp,
width,
out_height,
true,
"best_space_dir-rt"+space_rt+"-rxy"+space_rxy,
titles);
}
return;
}
public static double [][] getBestSpace(
final double [][][] space_conv,
final int [][] best_space){ // null or int [nimg][]
final int num_kernels = space_conv.length; // may be modified, if there will be extra layers
final int num_images = space_conv[0].length;
final int num_pixels = space_conv[0][0].length;
final double [][] best_val = new double [num_images][num_pixels];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
// final AtomicInteger ati = new AtomicInteger(0);
for (int nimg = 0; nimg < num_images; nimg++) {
final int fimg = nimg;
Arrays.fill(best_val[fimg], Double.NaN);
if (best_space != null) {
best_space[fimg] = new int[space_conv[0][0].length];
Arrays.fill(best_space[fimg], -1);
}
ai.set(0);
// ati.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
// int nthread = ati.getAndIncrement();
for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
int best_dir = -1;
int best_sgn = 0;
double bestv = 0;
for (int dir = 0; dir < num_kernels; dir++) {
double d = space_conv[dir][fimg][ipix];
if (!Double.isNaN(d)) {
for (int sgn = 0; sgn<2; sgn++) {
if (sgn > 0) {
d = -d;
}
if (d > bestv) {
best_dir = dir;
best_sgn = sgn;
bestv = d;
}
}
}
}
if (best_dir >= 0) {
best_val[fimg][ipix] = bestv;
if (best_space != null) {
if (best_dir < 4) {
best_space[fimg][ipix] = best_dir + best_sgn * 4;
} else {
best_space[fimg][ipix] = 10 + best_sgn;
}
//best_space[fimg][ipix] = best_dir + best_sgn * num_kernels;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
return best_val;
}
public static double [] normalizeSpaceKernels(
double [][][] conv_data) { // same data convolved with different kernels, may have NaNs
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final AtomicInteger ati = new AtomicInteger(0);
final double [][] s0t = new double [threads.length][conv_data.length];
final double [][] sdt = new double [threads.length][conv_data.length];
final double [][] sdat = new double [threads.length][conv_data.length]; // will compare l2 and abs
final double [][] sd2t = new double [threads.length][conv_data.length];
final double [] s0 = new double [conv_data.length];
final double [] sd = new double [conv_data.length];
final double [] sda = new double [conv_data.length]; // will compare l2 and abs
final double [] sd2 = new double [conv_data.length];
final double [] std = new double [conv_data.length];
final double [] coeffs_std = new double [conv_data.length+1];
final double [] coeffs_abs = new double [conv_data.length+1];
for (int nkern = 0; nkern < conv_data.length; nkern++) {
final int fnkern = nkern;
for (int nimg = 0; nimg < conv_data[fnkern].length; nimg++) {
final int fimg = nimg;
ai.set(0);
ati.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
int nthread = ati.getAndIncrement();
for (int ipix = ai.getAndIncrement(); ipix < conv_data[fnkern][fimg].length; ipix = ai.getAndIncrement()) {
double d = conv_data[fnkern][fimg][ipix];
if (!Double.isNaN(d)) {
s0t[nthread][fnkern] += 1.0;
sdt[nthread][fnkern] += d;
sd2t[nthread][fnkern] += d*d;
sdat[nthread][fnkern] += Math.abs(d);
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
for (int nthread = 0; nthread < s0t.length; nthread++) {
s0 [nkern] += s0t [nthread][nkern];
sd [nkern] += sdt [nthread][nkern];
sd2[nkern] += sd2t[nthread][nkern];
sda[nkern] += sdat[nthread][nkern];
}
}
double avg_std = 0.0;
double avg_abs = 0.0;
for (int nkern = 0; nkern < conv_data.length; nkern++) {
sd[nkern] /= s0[nkern];
sda[nkern] /= s0[nkern];
sd2[nkern] /= s0[nkern];
std[nkern] = Math.sqrt(sd2[nkern] - sd[nkern] * sd[nkern]);
avg_std+= std[nkern];
avg_abs+= sda[nkern];
}
avg_std/=conv_data.length;
avg_abs/=conv_data.length;
for (int nkern = 0; nkern < conv_data.length; nkern++) {
coeffs_std[nkern] = avg_std/ std[nkern];
coeffs_abs[nkern] = avg_abs/ sda[nkern];
}
coeffs_std[conv_data.length] = avg_std; // last element - average
coeffs_abs[conv_data.length] = avg_abs;
return coeffs_std;
}
public static double [] normalizeTimeDerivs( // return {std, average abs()}
final double [][] conv_data,
final double frac_outliers) {
final int num_bins = 10000;
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final AtomicInteger ati = new AtomicInteger(0);
final double [] max_abs = new double [threads.length];
for (int nimg = 0; nimg < conv_data.length; nimg++) {
final double [] data = conv_data[nimg];
ai.set(0);
ati.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
int nthread = ati.getAndIncrement();
for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) if (!Double.isNaN(data[ipix])) {
max_abs[nthread] = Math.max(max_abs[nthread],Math.abs(data[ipix]));
}
}
};
}
ImageDtt.startAndJoin(threads);
}
double amax = 0;
for (int nthread = 0; nthread < max_abs.length; nthread++) {
amax=Math.max(amax, max_abs[nthread]);
}
// final double fmax = amax;
// final double scale = amax/num_bins;
final double step = amax/num_bins;
final double [][] histograms = new double [threads.length][num_bins];
for (int nimg = 0; nimg < conv_data.length; nimg++) {
final double [] data = conv_data[nimg];
ai.set(0);
ati.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
int nthread = ati.getAndIncrement();
for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) if (!Double.isNaN(data[ipix])) {
double ad = Math.abs(data[ipix]);
int bin = (int) Math.floor(ad/step);
if ((bin >=0) && (bin < num_bins)) {
histograms[nthread][bin] += 1.0;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
final double [] histogram = new double [num_bins];
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int bin = ai.getAndIncrement(); bin < num_bins; bin = ai.getAndIncrement()) {
for (int i = 0; i < histograms.length; i++) {
histogram[bin]+=histograms[i][bin];
}
}
}
};
}
ImageDtt.startAndJoin(threads);
double sum = 0;
for (int i = 0; i < num_bins; i++) {
sum+=histogram[i];
}
double threshold = sum * frac_outliers;
double run_sum = 0;
double prev_sum = 0;
double alim = 0;
for (int bin= num_bins-1; bin >=0; bin--) {
prev_sum = run_sum;
run_sum += histogram[bin];
if (run_sum >=threshold) {
alim = amax *(bin + (run_sum - threshold)/(run_sum - prev_sum))/num_bins;
break;
}
}
final double flim = alim;
// alim - absolute limit to remove outliers
final double [] s0t = new double [threads.length];
final double [] sdt = new double [threads.length];
final double [] sdat = new double [threads.length]; // will compare l2 and abs
final double [] sd2t = new double [threads.length];
for (int nimg = 0; nimg < conv_data.length; nimg++) {
final double [] data = conv_data[nimg];
ai.set(0);
ati.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
int nthread = ati.getAndIncrement();
for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) {
double d = data[ipix];
if (!Double.isNaN(d) && (Math.abs(d) < flim)) {
s0t[nthread] += 1.0;
sdt[nthread] += d;
sd2t[nthread] += d*d;
sdat[nthread] += Math.abs(d);
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
double s0=0, sd=0,sd2=0,sda=0;
for (int nthread = 0; nthread < s0t.length; nthread++) {
s0 += s0t [nthread];
sd += sdt [nthread];
sd2 += sd2t[nthread];
sda += sdat[nthread];
}
sd /= s0;
sda /= s0;
sd2 /= s0;
double std = Math.sqrt(sd2 - sd * sd);
return new double [] {std, sda};
}
public static void scaleConvolutions(
final double [][] conv_data, // same data convolved with different kernels, may have NaNs
final double scale) {
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
for (int nimg = 0; nimg < conv_data.length; nimg++) {
final double [] data = conv_data[nimg];
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < data.length; ipix = ai.getAndIncrement()) {
data[ipix] *= scale;
}
}
};
}
ImageDtt.startAndJoin(threads);
}
}
public static double [][] timeToY(
final double [][] multi_image,
final int width,
final int num_vert,
final int vgap){
final int height = multi_image[0].length/width;
final int num_img_in = multi_image.length - 1; //last slice is average
final int num_slices = (height+num_vert-1)/ num_vert;
final int out_height = num_img_in*num_vert + vgap*(num_vert-1);
final double [][] out_data = new double [num_slices][width * out_height];
for (int i = 0; i < out_data.length; i++) {
Arrays.fill(out_data[i], Double.NaN);
}
for (int nimg = 0; nimg <num_img_in; nimg++) {
for (int nrow = 0; nrow < height; nrow++) {
int nout_slice = nrow/num_vert;
int nsub = nrow%num_vert;
int offs_out = ((num_img_in + vgap)* nsub + nimg)*width;
// copy row
System.arraycopy(
multi_image[nimg],
nrow*width,
out_data[nout_slice],
offs_out,
width);
}
}
return out_data;
}
public static void showKernel(
double [][][] kernel,
String title) {
final int rT = (kernel.length-1)/2;
if (kernel.length!= (2 * rT+1)) {
throw new IllegalArgumentException("Not all kernel dimensions are odd");
}
int height = kernel[0].length;
int width = kernel[0][0].length;
String [] titles = new String [kernel.length];
double [][] data = new double [kernel.length][width*height];
for (int t = 0; t < kernel.length; t++) {
titles[t] = "t "+(t-rT);
for (int row = 0; row < kernel[0].length; row++) {
System.arraycopy(
kernel[t][row],
0,
data[t],
width*row,
width);
}
}
ShowDoubleFloatArrays.showArrays(
data,
width,
height,
true,
title,
titles);
}
public static double [][] convolve3d(
final double [][] data_in, // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
final int width,
final double [][][] kernel, // (odd)x(odd)x(odd}
final boolean repeat_border){
final int height = data_in[0].length/width;
final int rT = (kernel.length-1)/2;
final int rY = (kernel[0].length-1)/2;
final int rX = (kernel[0][0].length-1)/2;
if ((kernel.length!= (2 * rT+1)) || (kernel[0].length!= (2 * rY+1)) || (kernel[0][0].length!= (2 * rX+1))) {
throw new IllegalArgumentException("Not all kernel dimensions are odd");
}
final double [][] data_out = new double [data_in.length][data_in[0].length];
for (int i = 0; i < data_out.length; i++) {
Arrays.fill(data_out[i], Double.NaN);
}
final int lenXY = data_in[0].length;
final int lenT = data_in.length;
final int total_pix = lenXY * lenT;
final int lenTm1 = lenT-1;
final int lenYm1 = height-1;
final int lenXm1 = width-1;
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < total_pix; ipix = ai.getAndIncrement()) {
int iT = ipix / lenXY;
int iXY = ipix % lenXY;
if (!Double.isNaN(data_in[iT][iXY])) {
int iY = iXY / width;
int iX = iXY % width;
int t0 = iT-rT, t1 = iT+rT, y0 = iY-rY, y1 = iY + rY, x0= iX-rX, x1=iX+rX;
int kTindx = 2*rT, kYindx = 2*rY, kXindx = 2*rX;
boolean limited = (t0 < 0) || (t1 >= lenT) || (y0 < 0) || (y1 >= height) || (x0 < 0) || (x1 >= width);
if (limited) {
if (!repeat_border) {
if (t0 < 0) {
kTindx += t0; // negative
t0 = 0;
}
if (t1 >= lenT) t1 = lenT-1;
if (y0 < 0) {
kYindx += y0; // negative
y0 = 0;
}
if (y1 >= height) y1 = height-1;
if (x0 < 0) {
kXindx += x0; // negative
x0 = 0;
}
if (x1 >= width) x1 = width-1;
}
int kTi = kTindx;
double s = 0;
for (int tt = t0; tt <= t1; tt++) {
int t = Math.min(lenTm1, Math.max(tt, 0)); // if repeating
int kYi = kYindx;
for (int yy = y0; yy <= y1; yy++ ) {
int y = Math.min(lenYm1, Math.max(yy, 0)); // if repeating
int kXi = kXindx;
for (int xx = x0; xx <= x1; xx++ ) {
int x = Math.min(lenXm1, Math.max(xx, 0)); // if repeating
int xy = x + width * y;
double d = data_in[t][xy];
if (!Double.isNaN(d)) {
s += d*kernel[kTi][kYi][kXi];
}
kXi--;
}
kYi--;
}
kTi--;
}
data_out[iT][iXY] = s;
} else { // same faster, w/o tests
int kTi = kTindx;
double s = 0;
for (int t = t0; t <= t1; t++) {
int kYi = kYindx;
for (int y = y0; y <= y1; y++ ) {
int kXi = kXindx;
for (int x = x0; x <= x1; x++ ) {
int xy = x + width * y;
double d = data_in[t][xy];
if (!Double.isNaN(d)) {
s += d*kernel[kTi][kYi][kXi];
}
kXi--;
}
kYi--;
}
kTi--;
}
data_out[iT][iXY] = s;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return data_out;
}
public static double [][][] kernelDiffTime(
double rt,
double rxy) {
int irt = (int) Math.floor(rt);
int irxy = (int) Math.floor(rxy);
int lt = 2 * irt + 1;
int lxy = 2 * irxy + 1;
double [][][] kernel = new double [lt][lxy][lxy];
for (int it = 0; it < lt; it++) {
double t = it - irt;
double wt = -t * Math.cos(0.5*Math.PI* t / rt);
for (int iy = 0; iy < lxy; iy++) {
double y = iy- irxy;
double wty = wt * Math.cos(0.5*Math.PI* y / rxy);
for (int ix = 0; ix < lxy; ix++) {
double x = ix- irxy;
kernel[it][iy][ix] = wty * Math.cos(0.5*Math.PI* x / rxy);
}
}
}
double [][][] norm_kernel = kernelNormalize2(kernel);
return norm_kernel;
}
public static double [][][] kernelXY(
int dir, // 0: +y, 1: +x+y, 2:+x, 3 +x-t, 4: 0: +center will be voting 1 of 10 (each plus and minus), weighted by abs(kernelDiffTime)
double rt,
double rxy) {
double rxy2 = rxy*rxy;
int irt = (int) Math.floor(rt);
int irxy = (int) Math.floor(rxy);
int lt = 2 * irt + 1;
int lxy = 2 * irxy + 1;
double [][][] kernel = new double [lt][lxy][lxy];
for (int it = 0; it < lt; it++) {
double sum_cos = 0;
double sum_patt = 0;
double [][] cos_patt = new double [lxy][lxy];
double t = it - irt;
double wt = Math.cos(0.5*Math.PI* t / rt);
for (int iy = 0; iy < lxy; iy++) {
double y = iy- irxy;
// double wty = wt * Math.cos(0.5*Math.PI* y / rxy);
for (int ix = 0; ix < lxy; ix++) {
double x = ix- irxy;
double r2=(x*x+y*y)/rxy2;
double w = 0;
if (r2 < 1) {
double r = Math.sqrt(r2);
w = wt * Math.cos(0.5*Math.PI * r);
cos_patt[iy][ix] = w;
sum_cos += w;
switch (dir) {
case 0: w *= y; break;
case 1: w *= y - x; break;
case 2: w *= -x; break;
case 3: w *= -y - x; break;
case 4: w *= Math.cos(Math.PI * r); break; // center
}
sum_patt+=w;
}
kernel[it][iy][ix] = w;
}
}
double cos_corr = sum_patt/sum_cos; // balance average=0
for (int iy = 0; iy < lxy; iy++) {
for (int ix = 0; ix < lxy; ix++) {
kernel[it][iy][ix] -= cos_corr*cos_patt[iy][ix];
}
}
}
return kernelNormalize2(kernel);
}
public static double [][][] kernelNormalize2(double [][][] kernel){
double [][][] kernel_out = new double [kernel.length][kernel[0].length][kernel[0][0].length];
double s = 0;
for (int t = 0; t < kernel.length; t++) {
for (int y = 0; y < kernel[t].length; y++) {
for (int x = 0; x < kernel[t][y].length; x++) {
double d = kernel[t][y][x];
s+=d * d;
}
}
}
if (s == 0 ) {
throw new IllegalArgumentException("Kernel of all zeros");
}
double k = 1.0/Math.sqrt(s);
for (int t = 0; t < kernel.length; t++) {
for (int y = 0; y < kernel[t].length; y++) {
for (int x = 0; x < kernel[t][y].length; x++) {
kernel_out[t][y][x] = kernel[t][y][x] * k;
}
}
}
return kernel_out;
}
public static double [] applyWeights(
final double [][] data_in,
final double [][] weights) {
final double [] data_out = new double [data_in[0].length];
final double [] sum_weights = new double [data_in[0].length];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
for (int nimg = 0; nimg < weights.length; nimg++) {
final int fimg = nimg;
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < data_in[fimg].length; ipix = ai.getAndIncrement()) {
double w = weights[fimg][ipix];
if (w > 0 ) {
double d = data_in[fimg][ipix];
if (!Double.isNaN(d)) {
sum_weights[ipix] += w;
data_out[ipix] += w * d;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < data_out.length; ipix = ai.getAndIncrement()) {
if (sum_weights[ipix] > 0) {
data_out[ipix] /= sum_weights[ipix];
} else {
data_out[ipix] = Double.NaN;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return data_out;
}
public static void testOrange() {
boolean zero_is_nan=true;
int master = 3;
final double [] voffsets = value_offsets;
int [][] segments = match_segments; // full_segments; // single_segments;
double oscale = 2.0;
// Rectangle woi = new Rectangle(0,0,640,512); // full window;
// Rectangle woi = new Rectangle(92,223,320,240); // approx. overlap
Rectangle woi = new Rectangle(92,280,320,120); // minimal overlap
final Rectangle owoi = new Rectangle (
(int) Math.round(oscale *woi.x ),
(int) Math.round(oscale *woi.y ),
(int) Math.round(oscale *woi.width),
(int) Math.round(oscale *woi.height));
int [][] seg_offs_len = new int [segments.length][2];
int offs = 0;
for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
seg_offs_len[nimg][0] = offs;
int len = segments[nimg][1] - segments[nimg][0]; // length;
seg_offs_len[nimg][1] = len;
offs += len;
}
final int indx_avg = offs;
final int olen = owoi.width * owoi.height;
final int owidth = owoi.width;
final double [][] out_data = new double [offs+1][owoi.width * owoi.height];
final int [] num_valid = new int [owoi.width * owoi.height];
String [] titles = new String[out_data.length];
titles [indx_avg] = "average";
for (int i = 0; i < offs; i++) {
Arrays.fill(out_data[i], Double.NaN);
}
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
/*
faded_textures[nslice][0] = TileProcessor.fillNaNs(
faded_textures[nslice][0], // final double [] data,
null, // final boolean [] prohibit,
width, // int width,
// CAREFUL ! Remaining NaN is grown by unsharp mask filter ************* !
100, // 2*width, // 16, // final int grow,
0.7, // double diagonal_weight, // relative to ortho
100, // int num_passes,
0.03, // final double max_rchange, // = 0.01 - does not need to be accurate
THREADS_MAX); // final int threadsMax) // maximal number of threads to launch
*/
for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
final double voffs = +voffsets[nimg];
final int foffs = seg_offs_len[nimg][0];
ImagePlus imp_src = new ImagePlus(src_files[nimg]);
ImageStack stack = imp_src.getImageStack();
final int src_width = imp_src.getWidth();
final int src_height = imp_src.getHeight();
String [] stack_labels = stack.getSliceLabels();
final double [][] dpixels = new double [seg_offs_len[nimg][1]][];
for (int i = 0; i < seg_offs_len[nimg][1]; i++) {
titles[seg_offs_len[nimg][0] + i] =nimg+":"+stack_labels[segments[nimg][0]+i];
float [] fpixels = (float[]) stack.getPixels(segments[nimg][0]+i+1);
dpixels[i] = new double [fpixels.length];
for (int j = 0; j < dpixels[i].length; j++) {
dpixels[i][j] = fpixels[j];
if (zero_is_nan && (dpixels[i][j] == 0.0)) {
dpixels[i][j] = Double.NaN;
}
}
// fpixels[i] = (float[]) stack.getPixels(segments[nimg][0]+i+1);
}
final double [][] points0 = new double[matching_points[master].length][];
final double [][] points1 = new double[points0.length][]; // matching_points[nimg];
for (int i = 0; i < points0.length;i++) {
points0[i] = new double [] {matching_points[master][i][0]/points_scale, matching_points[master][i][1]/points_scale};
points1[i] = new double [] {matching_points[nimg ][i][0]/points_scale, matching_points[nimg ][i][1]/points_scale};
}
//points_scale
double [][] ab1 = ComboMatch.getAffineMatrix(
points0, // double [][] x,
points1); // double [][] y){
double k = 1.0/oscale;
final double [][] a1 = {{k*ab1[0][0],k*ab1[0][1]},{k*ab1[1][0],k*ab1[1][1]}};
// final double [] b1 = {k*ab1[2][0] + woi.x, k*ab1[2][1]+ woi.y}; // from output coordinates to image2 coordinates
// final double [] b1 = {ab1[2][0] + woi.x, ab1[2][1]+ woi.y}; // from output coordinates to image2 coordinates
final double [] b1 = {
ab1[2][0] + a1[0][0]*owoi.x + a1[0][1]*owoi.y,
ab1[2][1] + a1[1][0]*owoi.x + a1[1][1]*owoi.y}; // from output coordinates to image2 coordinates
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < olen; ipix = ai.getAndIncrement()) {
int y = ipix / owidth; // output image coordinates
int x = ipix % owidth; // 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 ix0 = (int) Math.floor(fx);
int iy0 = (int) Math.floor(fy);
int ix1 = ix0 + 1;
int iy1 = iy0 + 1;
if ((ix1 >=0) && (iy1 >= 0) && (ix0 < src_width) && (iy0 < src_height)) {
double dx = fx - ix0;
double dy = fy - iy0;
if (ix0 < 0) ix0 = 0;
if (iy0 < 0) iy0 = 0;
if (ix1 >= src_width) ix1 = src_width - 1;
if (iy1 >= src_height) iy1 = src_height - 1;
int indx00 = ix0 + iy0 * src_width;
int indx01 = ix1 + iy0 * src_width;
int indx10 = ix0 + iy1 * src_width;
int indx11 = ix1 + iy1 * src_width;
for (int islice = 0; islice < dpixels.length; islice++) {
int oslice = islice + foffs;
double d00 = dpixels[islice][indx00];
double d01 = dpixels[islice][indx01];
double d10 = dpixels[islice][indx10];
double d11 = dpixels[islice][indx11];
double d = d00*(1-dx)*(1-dy)+d01*dx*(1-dy)+d10*(1-dx)*dy+d11*dx*dy + voffs;
out_data[oslice][ipix]= d;
if (!Double.isNaN(d)) {
out_data[indx_avg][ipix] += d;
num_valid[ipix]++;
}
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
for (int ipix=0; ipix < olen; ipix++) {
if (num_valid[ipix] == 0) {
out_data[indx_avg][ipix] = Double.NaN;
} else {
out_data[indx_avg][ipix]/=num_valid[ipix];
}
}
ShowDoubleFloatArrays.showArrays(
out_data,
owoi.width,
owoi.height,
true,
"merged_seq_scale"+oscale+"-slices"+indx_avg+"_x"+woi.x+"_y"+woi.y+"_w"+woi.width+"_h"+woi.height+(zero_is_nan?"_zeronan":""),
titles);
return;
}
}
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