Skip to content

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

experimenting with different modes

parent c9ed0a3d
Hide whitespace changes
Inline Side-by-side
Showing with 1737 additions and 213 deletions
src/main/java/com/elphel/imagej/common/ShowDoubleFloatArrays.java
View file @ be7873b5
...@@ -91,6 +91,8 @@ import ij.process.ImageProcessor; ...@@ -91,6 +91,8 @@ import ij.process.ImageProcessor;
return; return;
} else showArrays(pixels, width, height, title); } else showArrays(pixels, width, height, title);
} }
public static void showArrays(double[][] pixels, int width, int height, boolean asStack, String title, String [] titles) { public static void showArrays(double[][] pixels, int width, int height, boolean asStack, String title, String [] titles) {
int i,j; int i,j;
if (pixels==null) { if (pixels==null) {
...@@ -118,6 +120,99 @@ import ij.process.ImageProcessor; ...@@ -118,6 +120,99 @@ import ij.process.ImageProcessor;
return; return;
} else showArrays(pixels, width, height, titles); } else showArrays(pixels, width, height, titles);
} }
public static ImagePlus showArraysHyperstack(double[][] pixels, int width, int height, int slices, String title, String [] titles, boolean show) {
int i,j;
if (pixels==null) {
System.out.println("showDoubleFloatArrays.showArrays(): - pixel array is null");
}
float [] fpixels;
ImageStack array_stack=new ImageStack(width,height);
boolean not_empty = false;
for (i=0;i<pixels.length;i++) if (pixels[i]!=null) {
not_empty = true;
fpixels=new float[pixels[i].length];
for (j=0;j<fpixels.length;j++) fpixels[j]=(float)pixels[i][j];
if (i < titles.length) {
array_stack.addSlice(titles[i], fpixels);
} else {
array_stack.addSlice("slice-"+i, fpixels);
}
}
ImagePlus imp_stack = null;
if (not_empty) {
imp_stack = IJ.createImage(
title, // String title,
"32-bit,grayscale-mode", // ,label", // String type,
width, // int width,
height, // int height,
1, // int channels,
slices, // int slices,
pixels.length / slices); // frames)
imp_stack.setStack(array_stack);
imp_stack.getProcessor().resetMinAndMax();
if (show) {
imp_stack.show();
}
}
return imp_stack;
}
public static ImagePlus showArraysHyperstack(
double[][][] pixels,
int width,
String title,
String [] titles, // all slices*frames titles or just slice titles or null
String [] frame_titles, // frame titles or null
boolean show) {
int num_frames = pixels.length;
int num_slices = pixels[0].length;
double [][] dpixels = new double [num_frames*num_slices][];
// System.out.println("pixels.length="+pixels.length+" pixels[0].length="+pixels[0].length);
for (int f = 0; f < num_frames; f++) {
// System.out.println("f="+f);
for (int s = 0; s < num_slices; s ++) {
int indx = s + f * num_slices;
// System.out.println("f="+f+" s="+s+" indx="+indx);
dpixels[indx] = pixels[f][s];
}
}
String [] combo_titles;
if ((titles != null) && (titles.length == dpixels.length)) {
combo_titles = titles;
} else {
combo_titles = new String[dpixels.length];
if (titles == null) {
titles = new String [pixels[0].length];
for (int i = 0; i < titles.length; i++) {
titles[i] = ""+i;
}
}
if (frame_titles == null) {
frame_titles = new String [pixels.length];
for (int i = 0; i < titles.length; i++) {
frame_titles[i] = ""+i;
}
}
for (int f = 0; f < frame_titles.length; f++) {
for (int s = 0; s < titles.length; s++) {
combo_titles[s + f * titles.length] = frame_titles[f]+":"+titles[s];
}
}
}
return showArraysHyperstack(
dpixels, // double[][] pixels,
width, // int width,
pixels[0][0].length/width, // int height,
pixels[0].length, // int slices,
title, // String title,
combo_titles, // String [] titles);
show); // boolean show
}
public static void showComplex(double[][][] cpixels, String title) { public static void showComplex(double[][][] cpixels, String title) {
int height = cpixels.length; int height = cpixels.length;
int width = cpixels[0].length; int width = cpixels[0].length;
......
src/main/java/com/elphel/imagej/orthomosaic/OrangeTest.java
View file @ be7873b5
...@@ -4,12 +4,16 @@ import java.awt.Rectangle; ...@@ -4,12 +4,16 @@ import java.awt.Rectangle;
import java.util.Arrays; import java.util.Arrays;
import java.util.concurrent.atomic.AtomicInteger; import java.util.concurrent.atomic.AtomicInteger;
import com.elphel.imagej.common.DoubleGaussianBlur;
import com.elphel.imagej.common.ShowDoubleFloatArrays; import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.tileprocessor.ImageDtt; import com.elphel.imagej.tileprocessor.ImageDtt;
import com.elphel.imagej.tileprocessor.TileNeibs;
import com.elphel.imagej.tileprocessor.TileProcessor; import com.elphel.imagej.tileprocessor.TileProcessor;
import ij.ImagePlus; import ij.ImagePlus;
import ij.ImageStack; import ij.ImageStack;
import ij.Prefs;
import ij.io.FileSaver;
public class OrangeTest { public class OrangeTest {
public static double [][][] matching_points_initial= { public static double [][][] matching_points_initial= {
...@@ -128,7 +132,7 @@ public class OrangeTest { ...@@ -128,7 +132,7 @@ public class OrangeTest {
} }
}; };
public static int NUM_SPACE_KERNELS = 9;
public static String [] src_files = { 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/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/1697877490_930437-TERRAIN-MERGED.tiff",
...@@ -137,7 +141,9 @@ public class OrangeTest { ...@@ -137,7 +141,9 @@ public class OrangeTest {
"/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_314232-TERRAIN-MERGED.tiff",
"/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/1697877492_614332-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 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 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 String merged_file = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/04-third-corr/merged_weighted/merged_seq_scale2.0-slices260_x92_y280_w320_h120_zeronan.tiff";
public static String convolutions_dir = "/media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/04-third-corr/convolutions/";
public static double points_scale = 4.0; 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/ public static double [][][] matching_points= { // /media/elphel/SSD3-4GB/lwir16-proc/berdich3/debug/agent_orange/terrain-merged/02-first_point_correction/
...@@ -182,156 +188,559 @@ public class OrangeTest { ...@@ -182,156 +188,559 @@ public class OrangeTest {
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 [][] 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 int [][] single_segments = {{60,61},{0,1},{0,1},{0,1},{0,1},{0,1}};
public static double [][] patterns_integrate0 = {
{1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // first scale
{1.1, 1.1, 1.1, 1.1, 0.0, 0.0, 0.0, 0.0, 0.0}, // second scale
{1.2, 1.2, 1.2, 1.2, 0.0, 0.0, 0.0, 0.0, 0.4} // third scale
};
public static double [][] patterns_integrate1 = {
{1.0, 1.0, 1.0, 1.0, 0.65, 0.65, 0.65, 0.65, 0.0}, // first scale
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // second scale
{1.3, 1.3, 1.3, 1.3, 0.85, 0.85, 0.85, 0.85, 0.0} // third scale
};
public static double [][] patterns_integrate = { // only sharp
{1.0, 1.0, 1.0, 1.0, 0.65, 0.65, 0.65, 0.65, 0.0}, // first scale
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // second scale
{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0} // third scale
};
public static void processMerged() { public static void processMerged() {
boolean save_convolutions = true;
boolean convolve_and_exit = true;
boolean use_saved_convolutions = true; // false; // true;
boolean repeat_border= true;
boolean show_orig_dt = false; // true; // false;
boolean show_derivs_dt = true;
boolean show_kernels = true; // false; // true;
boolean full_convolves_only = true;
boolean show_pre_best = true; // false;
boolean show_best = true;
boolean show_best_temp = true;
boolean show_result = true;
int num_vert = 1; //4; int num_vert = 1; //4;
int vgap = 0; // 10; int vgap = 0; // 10;
double diff_time_rt = 2.9; // 5.5; double [][] time_rt_scales = {{2.9,2.9},{2.9,4.9},{4.9,4.9}, {9.9,2.9}};
double diff_time_rxy = 2.9; // 5.5; double [][] space_rt_scales = {{2.9,2.9},{2.9,4.9},{4.9,4.9}};
double [] blur3d = {2.9, 4.9}; // blur result weights before application
int num_tscales = time_rt_scales.length;
int num_sscales = space_rt_scales.length;
double space_rt = 2.9; // 5.5;
double space_rxy = 2.9; // 5.5;
double frac_outliers = 0.1; double frac_outliers = 0.1;
String conv_dir = save_convolutions ? convolutions_dir : null;
boolean repeat_border= true; // Create kernels
boolean show_orig_dt = false; double [][][][] kernels_diff_time= new double [num_tscales][][][];
boolean show_derivs_dt = true; for (int nk = 0; nk < num_tscales; nk++) {
kernels_diff_time[nk] = kernelDiffTime(
time_rt_scales[nk][0], //double rt,
time_rt_scales[nk][1]); // double rxy)
if (show_kernels) {
showKernel(
kernels_diff_time[nk], // double [][][] kernel,
"kernel_diff_time_num"+ nk + "_rt"+time_rt_scales[nk][0]+"_rxy"+time_rt_scales[nk][1]); // String title)
}
}
System.out.println("Created kernels");
double [][][][][] space_kernels= new double[num_sscales][NUM_SPACE_KERNELS][][][];
for (int nk = 0; nk < num_sscales; nk++) {
for (int dir = 0; dir < space_kernels[nk].length; dir++) {
space_kernels[nk][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_scales[nk][0], // double rt,
space_rt_scales[nk][1]); // double rxy)
}
if (show_kernels) {
showKernels(
space_kernels[nk], // double [][][][] kernels,
"space_kernels_num"+nk+"_rt"+space_rt_scales[nk][0]+"_rxy"+space_rt_scales[nk][1]); // String title)
}
}
double [][][] blur_kernel = kernelBlur3d(
blur3d[0], // double rt,
blur3d[1]); // double rxy)
if (show_kernels) {
showKernels(
new double [][][][] {blur_kernel}, // double [][][][] kernels,
"blur_kernel_rt"+blur3d[0]+"_rxy"+blur3d[1]); // String title)
}
// Read source data
ImagePlus imp_merged = new ImagePlus(merged_file); ImagePlus imp_merged = new ImagePlus(merged_file);
ImageStack stack = imp_merged.getImageStack(); ImageStack stack = imp_merged.getImageStack();
final int width = imp_merged.getWidth(); final int width = imp_merged.getWidth();
// final int merged_height = imp_merged.getHeight(); final int height = imp_merged.getHeight();
// String [] stack_labels = stack.getSliceLabels(); final double [][] dpixels_all = new double [stack.getSize()][];
final double [][] dpixels = new double [stack.getSize()][]; String [] titles_spatial_all = new String[stack.getSize()];
for (int i = 0; i < dpixels.length; i++) { for (int i = 0; i < stack.getSize(); i++) {
titles_spatial_all[i] = stack.getSliceLabel(i+1);
}
String [] titles_spatial = new String [titles_spatial_all.length - 2];
System.arraycopy(titles_spatial_all, 0, titles_spatial, 0, titles_spatial.length);
for (int i = 0; i < dpixels_all.length; i++) {
float [] fpixels = (float[]) stack.getPixels(i+1); float [] fpixels = (float[]) stack.getPixels(i+1);
dpixels[i] = new double[fpixels.length]; dpixels_all[i] = new double[fpixels.length];
for (int j = 0; j < fpixels.length; j++) { for (int j = 0; j < fpixels.length; j++) {
dpixels[i][j] = fpixels[j]; dpixels_all[i][j] = fpixels[j];
}
}
String [] titles_temporal = new String[dpixels_all[0].length/width]; // height
for (int t = 0; t < titles_temporal.length; t++) {
if (num_vert == 1) {
titles_temporal[t]=String.format("row %d",t);
} else {
titles_temporal[t]=String.format("rows: %d-%d",t*num_vert,(t+1)*num_vert-1);
} }
} }
if (show_orig_dt) { String [] titles_pattern = new String [num_sscales * 2 * NUM_SPACE_KERNELS];
double [][] out_data_dt = timeToY( for (int nscale = 0; nscale < num_sscales; nscale++) {
dpixels, // final double [][] multi_image, for (int nsign = 0; nsign<2; nsign++) {
for (int npatt = 0; npatt <NUM_SPACE_KERNELS; npatt++) {
int indx = npatt + nsign * NUM_SPACE_KERNELS + nscale * (NUM_SPACE_KERNELS * 2);
titles_pattern[indx] = String.format("%d: scale=%d, sign=%d, patt=%d", indx, nscale, nsign, npatt);
}
}
}
if (show_orig_dt) {
showOrigDt (
dpixels_all, // final double [][] multi_image,
width, // final int width, width, // final int width,
num_vert, // final int num_vert, num_vert, // final int num_vert,
vgap); // final int vgap) vgap); // final int vgap)
int out_height = out_data_dt[0].length/width; }
String [] titles = new String[out_data_dt.length]; final int remove_last = 2;
for (int i = 0; i < titles.length; i++) { final double [][] dpixels = new double[dpixels_all.length - remove_last][];
if (num_vert == 1) { System.arraycopy(dpixels_all, 0, dpixels, 0, dpixels.length);
titles[i]=String.format("row %d",i); final double [] dpixels_avg = dpixels_all[dpixels_all.length - 2];
} else { final double [] dpixels_weight = dpixels_all[dpixels_all.length - 1];
titles[i]=String.format("rows: %d-%d",i*num_vert,(i+1)*num_vert-1); System.out.println("Read source data");
// Convolve multi-scale temporal derivative kernels with source data
double [][][] time_derivs = new double[num_tscales][][];
double [][][][] space_conv = new double[num_sscales][NUM_SPACE_KERNELS][][];
String title_space_conv = "space_conv";
for (int nk = 0; nk < num_sscales; nk++) {
title_space_conv += "_"+space_rt_scales[nk][0]+"-"+space_rt_scales[nk][1];
}
String title_time_derivs = "time_derivs";
int num_temp_scales = num_tscales;
num_temp_scales = 3; //////////////////// Temporary
for (int nk = 0; nk < num_temp_scales; nk++) {
title_time_derivs += "_"+time_rt_scales[nk][0]+"-"+time_rt_scales[nk][1];
}
if (use_saved_convolutions) {
String dir = conv_dir;
if (!dir.endsWith(Prefs.getFileSeparator())){
dir+=Prefs.getFileSeparator();
}
double [][][] t_d = restoreTemporalConvolutions(dir + title_time_derivs+".tiff");
double [][][][] s_c = restoreSpaceConvolutions (dir + title_space_conv+".tiff");
System.arraycopy(t_d, 0, time_derivs, 0, t_d.length);
System.arraycopy(s_c, 0, space_conv, 0, s_c.length);
} else {
for (int nk = 0; nk < num_tscales; nk++) {
time_derivs[nk] = 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,
kernels_diff_time[nk], // final double [][][] kernel, // (odd)x(odd)x(odd}
repeat_border, // final boolean repeat_border)
full_convolves_only); // final boolean full_convolves_only)
}
// Convolve multi-scale space kernels with source data
System.out.println("Convolved temporal derivatives, created time_derivs");
double [][] space_scales = new double [num_sscales][];
for (int nk = 0; nk < num_sscales; nk++) {
for (int dir = 0; dir < space_kernels[nk].length; dir++) {
space_conv[nk][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[nk][dir], // final double [][][] kernel, // (odd)x(odd)x(odd}
repeat_border, // final boolean repeat_border)
full_convolves_only); // final boolean full_convolves_only)
}
space_scales[nk] = normalizeSpaceKernels(
space_conv[nk]); // double [][][] conv_data) // same data convolved with different kernels, may have NaNs
System.out.println("Convoved with kernels scale "+nk);
}
double average_space = space_scales[0][NUM_SPACE_KERNELS];
System.out.println("Space kernels normalization:");
for (int nk = 0; nk < num_sscales; nk++) {
System.out.println("Scale number "+nk);
for (int dir = 0; dir < space_kernels[nk].length; dir++) {
System.out.println("dir number "+dir+", scale "+space_scales[nk][dir]);
}
System.out.println("average value "+space_scales[nk][space_scales[nk].length-1]);
}
System.out.println("average_space = "+average_space);
// normalize all kernel groups for different scales together (match to the first group
for (int nk = 0; nk < num_sscales; nk++) {
for (int i = 0; i < NUM_SPACE_KERNELS; i++) {
space_scales[nk][i] *= average_space/space_scales[nk][NUM_SPACE_KERNELS];
}
for (int dir = 0; dir < NUM_SPACE_KERNELS; dir++) {
scaleConvolutions(
space_conv[nk][dir], // final double [][] conv_data, // same data convolved with different kernels, may have NaNs
space_scales[nk][dir]); // final double scale)
} }
} }
// normalize time derivatives to match those of the spatial ones
double [][] time_averages = new double [num_tscales][];
double [] time_std = new double [num_tscales];
for (int nk = 0; nk < num_tscales; nk++) {
time_averages[nk] = normalizeTimeDerivs( // return {std, average abs()}
time_derivs[nk], // final double [][] conv_data,
frac_outliers); // final double frac_outliers)
time_std[nk] = time_averages[nk][0]; // [1] - average of the absolute values;
scaleConvolutions(
time_derivs[nk], // final double [][] conv_data, // same data convolved with different kernels, may have NaNs
average_space/time_std[nk]); // final double scale)
}
System.out.println("Temporal kernels normalization:");
for (int nk = 0; nk < num_tscales; nk++) {
System.out.println("Scale number "+nk+" std="+time_averages[nk][0]+" abs="+time_averages[nk][1]);
}
saveSpaceConvolutions(
space_conv, // double [][][][] space_convolutions,
width, // int width,
title_space_conv, // String title,
conv_dir); // String save_dir) // do not save if null, but show
saveTemporalConvolutions(
time_derivs, // double [][][] time_derivs,
width, // int width,
title_time_derivs, // String title,
conv_dir); // String save_dir) // do not save if null, but show
if (convolve_and_exit) {
return;
}
}
// calculate time derivative with the selected time derivative kernel for each spatial convolution
System.out.println("Calculating space convolution derivatives with respect to time");
final int time_deriv_patt = 3; // convolve with a single derivative
double [][][][] space_conv_dt = new double [space_conv.length][space_conv[0].length][][];
for (int nscale = 0; nscale < space_conv.length; nscale++) {
for (int npatt = 0; npatt < space_conv[0].length; npatt++) {
space_conv_dt[nscale][npatt] = convolve3d(
space_conv[nscale][npatt], // final double [][] data_in, // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
width, // final int width,
kernels_diff_time[time_deriv_patt], // final double [][][] kernel, // (odd)x(odd)x(odd}
repeat_border, // final boolean repeat_border)
full_convolves_only); // final boolean full_convolves_only)
}
}
System.out.println("Done calculating space convolution derivatives with respect to time");
{
String title_space_dt = title_space_conv+"_dt_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1];
saveSpaceConvolutions(
space_conv_dt, // double [][][][] space_convolutions,
width, // int width,
title_space_dt, // String title,
null); // conv_dir); // String save_dir) // do not save if null, but show
}
// combine kernel derivatives fotr future exclusion from integration
boolean use_max = false; // true;
double [][] combo_dt = combineSelectedConvolutions(
space_conv_dt, // final double [][][][] space_conv_dt,
patterns_integrate, // final double [][] pattern_selection)
use_max); // final boolean use_max)
double [][] combo_dt1 = combineSelectedConvolutions(
space_conv_dt, // final double [][][][] space_conv_dt,
patterns_integrate, // final double [][] pattern_selection)
!use_max); // final boolean use_max)
{
String title_combo_dt = title_space_conv+"_combo_dt_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1]+(use_max?"-max":"-avg");
ShowDoubleFloatArrays.showArrays( ShowDoubleFloatArrays.showArrays(
out_data_dt, combo_dt,
width, width,
out_height, height,
true,
title_combo_dt+".tiff",
titles_spatial);
String title_combo_dt1 = title_space_conv+"_combo_dt_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1]+((!use_max)?"-max":"-avg");
ShowDoubleFloatArrays.showArrays(
combo_dt1,
width,
height,
true,
title_combo_dt1+".tiff",
titles_spatial);
}
double [][] weights_pre = weightsMinMax(
combo_dt, // final double [][] weights_in, // same data convolved with different kernels, may have NaNs
20.0, // final double min,
40.0, // final double max,
true); // final boolean invert)
{
String title_combo_dt = title_space_conv+"_weights-pre_"+time_rt_scales[time_deriv_patt][0]+"-"+time_rt_scales[time_deriv_patt][1]+(use_max?"-max":"-avg");
ShowDoubleFloatArrays.showArrays(
weights_pre,
width,
height,
true,
title_combo_dt+".tiff",
titles_spatial);
}
double [][] integrated_pattern_strengths = integrateSpaceConv(
space_conv, // final double [][][][] space_conv,
weights_pre); // final double [][] weights) // may be null [num_images][num_pixels]
if (show_best && show_pre_best) {
ShowDoubleFloatArrays.showArrays(
integrated_pattern_strengths,
width,
height,
true, true,
"time_section-"+out_data_dt.length, "integrated_pattern_strengths.tiff",
titles); titles_pattern);
}
int [][] best_pattern_frame = getBestConv(
space_conv); // final double [][][][] space_conv)
int patt_mod = 18; // = 54;
if (show_best && show_pre_best) {
double [][] best_pattern_frame_dbg = new double [best_pattern_frame.length][best_pattern_frame[0].length];
for (int nimg = 0; nimg < best_pattern_frame_dbg.length; nimg++) {
Arrays.fill(best_pattern_frame_dbg[nimg],Double.NaN);
for (int ipix = 0; ipix < best_pattern_frame_dbg[nimg].length; ipix++) if (best_pattern_frame[nimg][ipix] >= 0){
best_pattern_frame_dbg[nimg][ipix] = best_pattern_frame[nimg][ipix] % patt_mod;
}
}
ShowDoubleFloatArrays.showArrays(
best_pattern_frame_dbg,
width,
height,
true,
"best_pattern_frame_mod"+patt_mod+".tiff",
titles_pattern);
} }
double [][][] kernel_diff_time= kernelDiffTime( int [] best_pattern = getBestIntegratedConv(
diff_time_rt, //double rt, integrated_pattern_strengths); // double [][] integ_conv){ // [kern][pix]
diff_time_rxy); // double rxy) if (show_best && show_pre_best) {
double [] best_pattern_dbg = new double [best_pattern.length];
Arrays.fill(best_pattern_dbg,Double.NaN);
for (int ipix = 0; ipix < best_pattern_dbg.length; ipix++) if (best_pattern[ipix] >= 0){
best_pattern_dbg[ipix] = best_pattern[ipix] % patt_mod;
}
ShowDoubleFloatArrays.showArrays(
new double[][] {best_pattern_dbg},
width,
height,
true,
"best_pattern_mod"+patt_mod+".tiff",
titles_pattern);
}
double [][] patt_strength_abs = getPatternStrengths(
space_conv, // final double [][][][] space_conv,
best_pattern, // final int [] best_patt_integ, // best pattern index for integrated each pixel
null); // final int [][] best_patt_frame){ // best pattern index for each frame, each pixel. If !=null calculate strength relative to the best pattern, null - absolute
double [][] patt_strength_rel = getPatternStrengths(
space_conv, // final double [][][][] space_conv,
best_pattern, // final int [] best_patt_integ, // best pattern index for integrated each pixel
best_pattern_frame); // final int [][] best_patt_frame){ // best pattern index for each frame, each pixel. If !=null calculate strength relative to the best pattern, null - absolute
double strength_rel_min = 0.5;
double strength_rel_max = 0.9;
double [][] weights_mm = weightsMinMax(
patt_strength_rel, //final double [][] weights_in, // same data convolved with different kernels, may have NaNs
strength_rel_min, // final double min,
strength_rel_max, // final double max)
false); // final boolean invert) { // should be > min
// combine weights_pre (removed fg/bg borders with tempooral filtering) and weights_mm (by same pattern)
double [][] weights_mul = weightsMultiply(
weights_pre, // final double [][] weights0,
weights_mm); // final double [][] weights1)
// blur weights 3d with kernel
double [][] weights_blurred = convolve3d(
weights_mul, // final double [][] data_in, // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
width, // final int width,
blur_kernel, // final double [][][] kernel, // (odd)x(odd)x(odd}
repeat_border, // final boolean repeat_border)
full_convolves_only); // final boolean full_convolves_only)
if (show_result) {
double [][][] weights_dbg = {weights_pre, weights_mm, weights_mul, weights_blurred,dpixels};
String [] titles_weights= {"weights_pre", "weights_mm", "weights_mul", "weights_blurred","dpixels"};
ShowDoubleFloatArrays.showArraysHyperstack(
weights_dbg, // double[][][] pixels,
width, // int width,
"weights_combining", // String title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
titles_spatial, // String [] titles, // all slices*frames titles or just slice titles or null
titles_weights, // String [] frame_titles, // frame titles or null
true); // boolean show)
}
showKernel( double [][] applied_weight = applyWeights(
kernel_diff_time, // double [][][] kernel, dpixels, // final double [][] data_in,
"kernel_diff_time_rt"+diff_time_rt+"_rxy"+diff_time_rxy); // String title) weights_blurred); // weights_mm); // final double [][] weights)
double [][] time_derivs = convolve3d( if (show_result) {
dpixels, // final double [][] data_in, // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN double [][] result_img = {dpixels_avg, applied_weight[0], applied_weight[1]};
width, // final int width, String [] titles_result = {"average", "weighted", "sum.weights"};
kernel_diff_time, // final double [][][] kernel, // (odd)x(odd)x(odd} ShowDoubleFloatArrays.showArrays(
repeat_border); // final boolean repeat_border) result_img,
double [][][][] space_kernels= new double[5][][][]; width,
// double space_rt = 5.5; height,
// double space_rxy = 5.5; true,
for (int dir = 0; dir < space_kernels.length; dir++) { "integrated_frames_min"+strength_rel_min+"_max"+strength_rel_max+".tiff",
space_kernels[dir]= kernelXY( titles_result);
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()} if (show_best) {
time_derivs, // final double [][] conv_data, ShowDoubleFloatArrays.showArrays(
frac_outliers); // final double frac_outliers) patt_strength_abs,
double time_std = time_averages[0]; // [1] - average of the absolute values; width,
height,
true,
"patt_strength_abs.tiff",
titles_spatial);
ShowDoubleFloatArrays.showArrays(
patt_strength_rel,
width,
height,
true,
"patt_strength_rel.tiff",
titles_spatial);
}
for (int dir = 0; dir < space_kernels.length; dir++) { if (show_best_temp) {
scaleConvolutions( double [][] patt_strength_abs_temp = timeToY(
space_conv[dir], // final double [][] conv_data, // same data convolved with different kernels, may have NaNs patt_strength_abs, // final double [][] multi_image,
scales[dir]); // final double scale) width, // final int width,
} num_vert, // final int num_vert,
scaleConvolutions( vgap); // final int vgap)
time_derivs, // final double [][] conv_data, // same data convolved with different kernels, may have NaNs
average_space/time_std); // final double scale) double [][] patt_strength_rel_temp = timeToY(
patt_strength_rel, // final double [][] multi_image,
width, // final int width,
num_vert, // final int num_vert,
vgap); // final int vgap)
ShowDoubleFloatArrays.showArrays(
patt_strength_abs_temp,
width,
height,
true,
"patt_strength_abs_temp.tiff",
titles_temporal);
ShowDoubleFloatArrays.showArrays(
patt_strength_rel_temp,
width,
height,
true,
"patt_strength_rel_temp.tiff",
titles_temporal);
return;
}
double [][][][] space_conv_temp = new double[num_sscales][NUM_SPACE_KERNELS][][];
double [][][] time_derivs_temp = new double[num_tscales][][];
int [][] best_space_dir = new int [space_conv[0].length][];
// Find the best space kernel for each time and XY
int [][] best_space_dir = new int [space_conv[0][0].length][];
double [][] best_space_value = getBestSpace( double [][] best_space_value = getBestSpace(
space_conv, // final double [][][] space_conv, space_conv, // final double [][][][] space_conv
best_space_dir); // final int [][] best_space) // null or int [nimg][] best_space_dir); // final int [][] best_space) // null or int [nimg][]
System.out.println("Generated best_space_value");
if (show_derivs_dt) { if (show_derivs_dt) {
for (int dir = 0; dir < space_kernels.length; dir++) { for (int nk = 0; nk < num_sscales; nk++) {
space_conv_temp[dir] = timeToY( for (int dir = 0; dir < space_conv[nk].length; dir++) {
space_conv[dir], // final double [][] multi_image, space_conv_temp[nk][dir] = timeToY(
space_conv[nk][dir], // final double [][] multi_image,
width, // final int width,
num_vert, // final int num_vert,
vgap); // final int vgap)
// System.out.println("space_conv_temp["+nk+"]["+dir+"].length="+space_conv_temp[nk][dir].length);
}
}
// double [][][] time_derivs_temp = new double[num_tscales][][];
for (int nk = 0; nk < num_tscales; nk++) {
time_derivs_temp[nk] = timeToY(
time_derivs[nk], // final double [][] multi_image,
width, // final int width, width, // final int width,
num_vert, // final int num_vert, num_vert, // final int num_vert,
vgap); // final int vgap) vgap); // final int vgap)
}
int out_height = time_derivs_temp[0][0].length/width;
int num_slices = time_derivs_temp[0].length;
String [] frame_titles = new String[num_tscales];
for (int nk = 0; nk < num_tscales; nk++) {
frame_titles[nk]=String.format("kern %d: ",nk);
}
} // show time representations of temporal derivatives with different kernels
ShowDoubleFloatArrays.showArraysHyperstack(
double [][] time_derivs_temp = timeToY( time_derivs_temp, // double[][][] pixels,
time_derivs, // final double [][] multi_image, width, // int width,
width, // final int width, title_time_derivs+"-temp",// String title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
num_vert, // final int num_vert, titles_temporal, // String [] titles, // all slices*frames titles or just slice titles or null
vgap); // final int vgap) frame_titles, // String [] frame_titles, // frame titles or null
int out_height = time_derivs_temp[0].length/width; true); // boolean show)
String [] titles = new String[time_derivs_temp.length]; System.out.println("Displayed time_derivs_temp");
for (int i = 0; i < titles.length; i++) {
if (num_vert == 1) { double [][][] space_conv_temp_combo = new double [space_conv_temp.length*space_conv_temp[0].length][][];
titles[i]=String.format("row %d",i); System.out.println("space_conv_temp_combo.length="+space_conv_temp_combo.length);
} else { String [] frame_titles_combo = new String[space_conv_temp_combo.length];
titles[i]=String.format("rows: %d-%d",i*num_vert,(i+1)*num_vert-1); for (int nk = 0; nk < space_conv_temp.length; nk++) {
for (int dir = 0; dir < space_conv_temp[0].length; dir++) {
int indx = dir + nk*space_conv_temp[0].length;
space_conv_temp_combo[indx] = space_conv_temp[nk][dir];
frame_titles_combo[indx] = String.format("scale%d-dir%d",nk,dir);
// System.out.println("nk="+nk+" dir="+dir+" indx="+indx+" frame_titles_combo="+frame_titles_combo[indx]);
// System.out.println("space_conv_temp_combo["+indx+"].length="+space_conv_temp_combo[indx].length);
} }
} }
ShowDoubleFloatArrays.showArrays( System.out.println("title_space_conv="+title_space_conv);
time_derivs_temp, // System.out.println("space_conv_temp_combo.length="+space_conv_temp_combo.length);
width, //space_conv-dir3-rt5.5-rxy5.5.tif
out_height, ShowDoubleFloatArrays.showArraysHyperstack(
true, space_conv_temp_combo, // double[][][] pixels,
"time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy, width, // int width,
titles); title_space_conv, // String title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
for (int dir = 0; dir < space_kernels.length; dir++) { titles_temporal, // String [] titles, // "space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
ShowDoubleFloatArrays.showArrays( frame_titles_combo, // String [] frame_titles) // frame titles or null,
space_conv_temp[dir], true); // boolean show)
width, System.out.println("Displayed space_conv_temp_combo");
out_height,
true,
"space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
titles);
}
double [][] best_space_value_temp = timeToY( double [][] best_space_value_temp = timeToY(
best_space_value, // final double [][] multi_image, best_space_value, // final double [][] multi_image,
...@@ -339,13 +748,18 @@ public class OrangeTest { ...@@ -339,13 +748,18 @@ public class OrangeTest {
num_vert, // final int num_vert, num_vert, // final int num_vert,
vgap); // final int vgap) vgap); // final int vgap)
String title_best_space_value = "best_space_value";
for (int nk = 0; nk < num_sscales; nk++) {
title_best_space_value += "_"+space_rt_scales[nk][0]+"-"+space_rt_scales[nk][1];
}
ShowDoubleFloatArrays.showArrays( ShowDoubleFloatArrays.showArrays(
best_space_value_temp, best_space_value_temp,
width, width,
out_height, out_height,
true, true,
"best_space_value-rt"+space_rt+"-rxy"+space_rxy, title_best_space_value,
titles); titles_temporal);
System.out.println("Displayed best_space_value_temp");
double [][] best_space_dir_dbg = new double [best_space_dir.length][best_space_dir[0].length]; 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++) { for (int nimg = 0; nimg < best_space_dir.length; nimg++) {
...@@ -361,80 +775,563 @@ public class OrangeTest { ...@@ -361,80 +775,563 @@ public class OrangeTest {
num_vert, // final int num_vert, num_vert, // final int num_vert,
vgap); // final int vgap) vgap); // final int vgap)
String title_bdest_space_dir = "best_space_dir";
for (int nk = 0; nk < num_sscales; nk++) {
title_bdest_space_dir += "_"+space_rt_scales[nk][0]+"-"+space_rt_scales[nk][1];
}
ShowDoubleFloatArrays.showArrays( ShowDoubleFloatArrays.showArrays(
best_space_dir_temp, best_space_dir_temp,
width, width,
out_height, out_height,
true, true,
"best_space_dir-rt"+space_rt+"-rxy"+space_rxy, title_bdest_space_dir,
titles); titles_temporal);
System.out.println("Displayed best_space_dir_temp");
} }
return; return;
} }
public static double [][] getBestSpace( public static ImagePlus saveSpaceConvolutions(
final double [][][] space_conv, double [][][][] space_convolutions,
final int [][] best_space){ // null or int [nimg][] int width,
final int num_kernels = space_conv.length; // may be modified, if there will be extra layers String title,
final int num_images = space_conv[0].length; String save_dir) { // do not save if null, but show
final int num_pixels = space_conv[0][0].length; int num_scales = space_convolutions.length; // number of scales for kernel sets
int num_kernels = space_convolutions[0].length; // number of kernels for each scale
int num_times = space_convolutions[0][0].length; // number of time slices
// int num_pix = space_convolutions[0][0][0].length; // number of pixels in each slice
int num_slices = num_scales * num_kernels * num_times;
// System.out.println("num_scales="+num_scales+" num_kernels="+num_kernels+" num_times="+num_times);
String [] titles = new String[num_slices];
double [][][] data = new double [num_scales * num_kernels][][];
int indx = 0;
for (int s = 0; s < num_scales; s++) {
for (int k = 0; k < num_kernels; k++) {
for (int t = 0; t < num_times; t++) {
titles [indx++] = String.format("%d:%d:%d", s,k,t);
// System.out.println("titles[" + (indx - 1) + "]=" + titles[indx - 1]);
}
data[s * num_kernels + k] = space_convolutions[s][k];
}
}
ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
data, // double[][][] pixels,
width, // int width,
title, // String title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
titles, // String [] titles, // "space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
null, // String [] frame_titles) // frame titles or null,
(save_dir==null)); // boolean show)
if (save_dir != null) {
if (!save_dir.endsWith(Prefs.getFileSeparator())) {
save_dir += Prefs.getFileSeparator();
}
String file_path = save_dir+title+".tiff";
FileSaver fs=new FileSaver(imp);
fs.saveAsTiff(file_path);
}
return imp;
}
public static double [][][][] restoreSpaceConvolutions(String path){
ImagePlus imp = new ImagePlus(path);
final ImageStack stack = imp.getImageStack();
final int num_slices = stack.getSize(); // .length;
final String [] labels = new String[num_slices];
System.arraycopy(
stack.getSliceLabels(),
0,
labels,
0,
num_slices);
int nscales = 0; // number of scales for kernel sets
int nkernels = 0; // number of kernels for each scale
int ntimes = 0; // number of time slices
final int num_pixels = ((float[]) stack.getPixels(1)).length;
for (int i = 0; i < labels.length; i++) {
String [] tokens = labels[i].split(":");
int s = Integer.parseInt(tokens[0]);
int k = Integer.parseInt(tokens[1]);
int t = Integer.parseInt(tokens[2]);
nscales = Math.max(nscales, s);
nkernels = Math.max(nkernels, k);
ntimes = Math.max(ntimes, t);
}
final int num_scales = nscales+1; // number of scales for kernel sets
final int num_kernels = nkernels+1; // number of kernels for each scale
final int num_times = ntimes+1; // number of time slices
final double [][][][] data = new double [num_scales][num_kernels][num_times][num_pixels];
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 slice = ai.getAndIncrement(); slice < num_slices; slice = ai.getAndIncrement()) {
int s= slice / (num_kernels * num_times);
int kt = slice % (num_kernels * num_times);
int k = kt / num_times;
int t = kt % num_times;
float [] pixels = (float []) stack.getPixels(slice+1);
for (int i = 0; i < num_pixels; i++) {
data[s][k][t][i]=pixels[i];
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return data;
}
public static ImagePlus saveTemporalConvolutions(
double [][][] time_derivs,
int width,
String title,
String save_dir) { // do not save if null, but show
int num_scales = time_derivs.length; // number of scales for kernel sets
int num_times = time_derivs[0].length; // number of time slices
// int num_pix = time_derivs[0][0].length; // number of pixels in each slice
int num_slices = num_scales * num_times;
String [] titles = new String[num_slices];
int indx = 0;
for (int s = 0; s < num_scales; s++) {
for (int t = 0; t < num_times; t++) {
titles [indx++] = String.format("%d:%d", s,t);
}
}
ImagePlus imp = ShowDoubleFloatArrays.showArraysHyperstack(
time_derivs, // double[][][] pixels,
width, // int width,
title, // String title, "time_derivs-rt"+diff_time_rt+"-rxy"+diff_time_rxy,
titles, // String [] titles, // "space_conv-dir"+dir+"-rt"+space_rt+"-rxy"+space_rxy,
null, // String [] frame_titles) // frame titles or null,
(save_dir==null)); // boolean show)
if (save_dir != null) {
if (!save_dir.endsWith(Prefs.getFileSeparator())) {
save_dir += Prefs.getFileSeparator();
}
String file_path = save_dir+title+".tiff";
FileSaver fs=new FileSaver(imp);
fs.saveAsTiff(file_path);
}
return imp;
}
public static double [][][] restoreTemporalConvolutions(String path){
ImagePlus imp = new ImagePlus(path);
final ImageStack stack = imp.getImageStack();
final int num_slices = stack.getSize(); // .length;
final String [] labels = new String[num_slices];
System.arraycopy(
stack.getSliceLabels(),
0,
labels,
0,
num_slices);
int nscales = 0; // number of scales for kernel sets
int ntimes = 0; // number of time slices
final int num_pixels = ((float[]) stack.getPixels(1)).length;
for (int i = 0; i < labels.length; i++) {
// System.out.println(i+":"+labels[i]);
String [] tokens = labels[i].split(":");
int s = Integer.parseInt(tokens[0]);
int t = Integer.parseInt(tokens[1]);
nscales = Math.max(nscales, s);
ntimes = Math.max(ntimes, t);
}
final int num_scales = nscales+1; // number of scales for kernel sets
final int num_times = ntimes+1; // number of time slices
final double [][][] data = new double [num_scales][num_times][num_pixels];
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 slice = ai.getAndIncrement(); slice < num_slices; slice = ai.getAndIncrement()) {
int s= slice / num_times;
int t = slice % num_times;
float [] pixels = (float []) stack.getPixels(slice+1);
for (int i = 0; i < num_pixels; i++) {
data[s][t][i]=pixels[i];
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return data;
}
public static void showOrigDt (
final double [][] multi_image,
final int width,
final int num_vert,
final int vgap){
double [][] out_data_dt = timeToY(
multi_image, // 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);
}
public static double [][] getBestSpace(
final double [][][][] space_conv,
final int [][] best_space){ // null or int [nimg][]
final int num_scales = space_conv.length; // may be modified, if there will be extra layers
final int num_kernels = space_conv[0].length; // may be modified, if there will be extra layers
final int num_images = space_conv[0][0].length;
final int num_pixels = space_conv[0][0][0].length;
final double [][] best_val = new double [num_images][num_pixels]; final double [][] best_val = new double [num_images][num_pixels];
final Thread[] threads = ImageDtt.newThreadArray(); final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0); final AtomicInteger ai = new AtomicInteger(0);
// final AtomicInteger ati = new AtomicInteger(0);
for (int nimg = 0; nimg < num_images; nimg++) { for (int nimg = 0; nimg < num_images; nimg++) {
final int fimg = nimg; final int fimg = nimg;
Arrays.fill(best_val[fimg], Double.NaN); Arrays.fill(best_val[fimg], Double.NaN);
if (best_space != null) { if (best_space != null) {
best_space[fimg] = new int[space_conv[0][0].length]; best_space[fimg] = new int[num_pixels];
Arrays.fill(best_space[fimg], -1); Arrays.fill(best_space[fimg], -1);
} }
ai.set(0); ai.set(0);
// ati.set(0); // ati.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
public void run() { public void run() {
// int nthread = ati.getAndIncrement();
for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) { for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
int best_scale = -1;
int best_dir = -1; int best_dir = -1;
int best_sgn = 0; int best_sgn = 0;
double bestv = 0; double bestv = 0;
for (int dir = 0; dir < num_kernels; dir++) { for (int nscale = 0; nscale < num_scales; nscale++) {
double d = space_conv[dir][fimg][ipix]; for (int dir = 0; dir < num_kernels; dir++) {
if (!Double.isNaN(d)) { double d = space_conv[nscale][dir][fimg][ipix];
for (int sgn = 0; sgn<2; sgn++) { if (!Double.isNaN(d)) {
if (sgn > 0) { for (int sgn = 0; sgn<2; sgn++) {
d = -d; if (sgn > 0) {
d = -d;
}
if (d > bestv) {
best_scale = nscale;
best_dir = dir;
best_sgn = sgn;
bestv = d;
}
}
}
}
}
if (best_dir >= 0) {
best_val[fimg][ipix] = bestv;
if (best_space != null) {
if (best_dir < 8) {
best_space[fimg][ipix] = best_dir + best_sgn * 8 + best_scale * 16 ;
} else {
best_space[fimg][ipix] = 16* num_scales + best_sgn + 2 * best_scale;
}
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
return best_val;
}
public static int [][] getBestConv(
final double [][][][] space_conv){
final int num_scales = space_conv.length; // may be modified, if there will be extra layers
final int num_kernels = space_conv[0].length; // may be modified, if there will be extra layers
final int num_images = space_conv[0][0].length;
final int num_pixels = space_conv[0][0][0].length;
final int num_signs = 2;
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final int [][] best_conv = new int [num_images][num_pixels];
for (int nimg = 0; nimg < num_images; nimg++) {
final int fimg = nimg;
Arrays.fill(best_conv[fimg], -1);
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
int best_scale = -1;
int best_dir = -1;
int best_sgn = 0;
double bestv = 0;
for (int nscale = 0; nscale < num_scales; nscale++) {
for (int dir = 0; dir < num_kernels; dir++) {
double d = space_conv[nscale][dir][fimg][ipix];
if (!Double.isNaN(d)) {
for (int sgn = 0; sgn<2; sgn++) {
if (sgn > 0) {
d = -d;
}
if (d > bestv) {
best_scale = nscale;
best_dir = dir;
best_sgn = sgn;
bestv = d;
}
} }
if (d > bestv) { }
best_dir = dir; }
best_sgn = sgn; }
bestv = d; if (best_dir >= 0) {
best_conv[fimg][ipix] = best_dir + num_kernels * best_sgn + best_scale * num_kernels * num_signs;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
return best_conv;
}
public static double [][] integrateSpaceConv(
final double [][][][] space_conv,
final double [][] weights) { // may be null [num_images][num_pixels]
final int num_scales = space_conv.length; // may be modified, if there will be extra layers
final int num_kernels = space_conv[0].length; // may be modified, if there will be extra layers
final int num_images = space_conv[0][0].length;
final int num_pixels = space_conv[0][0][0].length;
final int num_signs = 2;
final double [][] integrated_convs = new double [num_scales*num_signs*num_kernels][num_pixels];
final double [] sum_weights = new double [num_pixels];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
for (int nimg = 0; nimg < num_images; 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 < num_pixels; ipix = ai.getAndIncrement()) {
double w = (weights == null)? 1.0 : weights[fimg][ipix] ;
sum_weights[ipix] += w;
for (int nscale = 0; nscale < num_scales; nscale++) {
for (int dir = 0; dir < num_kernels; dir++) {
double d = space_conv[nscale][dir][fimg][ipix];
if (!Double.isNaN(d)) {
for (int sgn = 0; sgn<num_signs; sgn++) {
if (sgn > 0) {
d = -d;
}
if (d > 0) {
if (weights != null) {
d *= weights[fimg][ipix];
}
int indx = dir + num_kernels * sgn + nscale * num_kernels * num_signs;
integrated_convs[indx][ipix] += d*d;
}
}
}
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
// normalize
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()) {
if (sum_weights[ipix] > 0 ) {
for (int i = 0; i < integrated_convs.length; i++) {
if (integrated_convs[i][ipix] >0) {
integrated_convs[i][ipix] = Math.sqrt(integrated_convs[i][ipix])/ sum_weights[ipix];
}
}
} else {
for (int i = 0; i < integrated_convs.length; i++) {
integrated_convs[i][ipix] = Double.NaN;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return integrated_convs;
}
public static int [] getBestIntegratedConv(
double [][] integ_conv){ // [kern][pix]
final int num_patt = integ_conv.length;
final int num_pixels = integ_conv[0].length;
final int [] best_index = new int [num_pixels];
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 < num_pixels; ipix = ai.getAndIncrement()) {
int best_patt = -1;
double bestv = 0;
for (int ipatt = 0; ipatt < num_patt; ipatt++) {
double d = integ_conv[ipatt][ipix];
if (d > bestv) { // cares of NaN
bestv = d;
best_patt = ipatt;
}
}
if (best_patt >=0) {
best_index[ipix] = best_patt;
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return best_index;
}
public static double [][] combineSelectedConvolutions(
final double [][][][] space_conv_dt,
final double [][] pattern_selection,
final boolean use_max) {
final int num_scales = space_conv_dt.length; // may be modified, if there will be extra layers
final int num_patterns = space_conv_dt[0].length; // may be modified, if there will be extra layers
final int num_images = space_conv_dt[0][0].length;
final int num_pixels = space_conv_dt[0][0][0].length;
if ((pattern_selection.length != num_scales) || (pattern_selection[0].length != num_patterns)) {
throw new IllegalArgumentException("pattern_selection does not match space_conv_dt in dimensions");
}
final double [][] combo_dt = new double [num_images][num_pixels];
for (int nimg = 0; nimg < num_images;nimg++) {
Arrays.fill(combo_dt[nimg], Double.NaN);
}
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 < num_pixels; ipix = ai.getAndIncrement()) {
for (int nimg = 0; nimg< num_images; nimg++) {
int num_valid = 0;
double s2 = 0;
for (int nscale = 0; nscale < num_scales; nscale++) {
for (int npatt = 0; npatt < num_patterns; npatt++) if (pattern_selection[nscale][npatt] > 0){
double d = space_conv_dt[nscale][npatt][nimg][ipix];
if (!Double.isNaN(d)) {
d*= pattern_selection[nscale][npatt]; // scaling d before squaring
d*=d;
if (use_max) {
s2 = Math.max(s2,d);
num_valid = 1;
} else {
s2 += d;
num_valid++;
} }
} }
} }
} }
if (best_dir >= 0) { if (num_valid > 0) {
best_val[fimg][ipix] = bestv; combo_dt[nimg][ipix] = Math.sqrt(s2/num_valid);
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; };
}
ImageDtt.startAndJoin(threads);
return combo_dt;
}
public static double [][] getPatternStrengths(
final double [][][][] space_conv,
final int [] best_patt_integ, // best pattern index for integrated each pixel
final int [][] best_patt_frame){ // best pattern index for each frame, each pixel. If !=null calculate strength relative to the best pattern, null - absolute
final boolean relative = (best_patt_frame != null); // false - just strength, true - relative to the best pattern for each frame/pix
// final int num_scales = space_conv.length; // may be modified, if there will be extra layers
final int num_patterns = space_conv[0].length; // may be modified, if there will be extra layers
final int num_images = space_conv[0][0].length;
final int num_pixels = space_conv[0][0][0].length;
final int num_signs = 2;
final double [][] pattern_strengths = new double [num_images][num_pixels];
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 < num_pixels; ipix = ai.getAndIncrement()) {
int indx = best_patt_integ[ipix];
if (indx >= 0) {
int npatt = indx % num_patterns;
int r = indx/ num_patterns;
int sgn = r % num_signs;
int nscale = r / num_signs;
for (int nimg = 0; nimg < num_images; nimg++) {
double s = space_conv[nscale][npatt][nimg][ipix];
if (!Double.isNaN(s)) {
if (sgn > 0) {
s = -s;
}
if (s > 0) {
if (relative) {
int indx_best =best_patt_frame[nimg][ipix];
if (indx_best >=0) {
int npatt_best = indx_best % num_patterns;
int r_best = indx_best / num_patterns;
int sgn_best = r_best % num_signs;
int nscale_best = r_best / num_signs;
double s_best = space_conv[nscale_best][npatt_best][nimg][ipix];
if (sgn_best > 0) {
s_best = -s_best;
}
pattern_strengths[nimg][ipix] = s/s_best;
}
} else {
pattern_strengths[nimg][ipix] = s;
}
} }
//best_space[fimg][ipix] = best_dir + best_sgn * num_kernels;
} }
} }
} }
} }
}; }
} };
ImageDtt.startAndJoin(threads); }
} ImageDtt.startAndJoin(threads);
return best_val; return pattern_strengths;
} }
// TODO: Normalize strengths with power and upper limit or fraction
public static double [] normalizeSpaceKernels( public static double [] normalizeSpaceKernels(
double [][][] conv_data) { // same data convolved with different kernels, may have NaNs double [][][] conv_data) { // same data convolved with different kernels, may have NaNs
...@@ -654,6 +1551,72 @@ public class OrangeTest { ...@@ -654,6 +1551,72 @@ public class OrangeTest {
} }
} }
public static double [][] weightsMultiply(
final double [][] weights0,
final double [][] weights1) {
final int num_frames = weights0.length;
final int num_pixels = weights0[0].length;
final double [][] weights = new double[num_frames][num_pixels];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < num_pixels; ipix = ai.getAndIncrement()){ // if (w_frame[ipix] > min){
for (int nimg = 0; nimg < num_frames; nimg++) {
double w = weights0[nimg][ipix]*weights1[nimg][ipix];
if (w > 0) { // NaN-aware
weights[nimg][ipix] = w;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return weights;
}
/*
for (int nimg = 0; nimg < num_frames; nimg++) {
final double [] w_frame = weights_in[nimg];
final double [] w_out = weights[nimg];
*/
public static double [][] weightsMinMax(
final double [][] weights_in, // same data convolved with different kernels, may have NaNs
final double min,
final double max,
final boolean invert) { // should be > min
final double scale = 1.0/(max-min);
final int num_frames = weights_in.length;
final int num_pixels = weights_in[0].length;
final double [][] weights = new double[num_frames][num_pixels];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
for (int nimg = 0; nimg < num_frames; nimg++) {
final double [] w_frame = weights_in[nimg];
final double [] w_out = weights[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 < num_pixels; ipix = ai.getAndIncrement()) if (!Double.isNaN(w_frame[ipix])){ // if (w_frame[ipix] > min){
double w = (Math.min(Math.max(w_frame[ipix] , min), max) - min) *scale;
w_out[ipix] = invert ? (1.0-w) : w;
}
}
};
}
ImageDtt.startAndJoin(threads);
}
return weights;
}
...@@ -717,12 +1680,53 @@ public class OrangeTest { ...@@ -717,12 +1680,53 @@ public class OrangeTest {
title, title,
titles); titles);
} }
public static void showKernels(
double [][][][] kernels,
String title) {
final int num_kernels = kernels.length;
final int num_times = kernels[0].length;
final int height = kernels[0][0].length;
final int width = kernels[0][0][0].length;
final int num_slices = num_kernels*num_times;
final int rT = (num_times - 1)/2;
if (kernels[0].length!= (2 * rT+1)) {
throw new IllegalArgumentException("Not all kernel dimensions are odd");
}
String [] titles = new String [num_slices];
double [][] data = new double [num_slices][width*height];
for (int k = 0; k < num_kernels; k++) {
for (int t = 0; t < num_times; t++) {
int slice = t+ num_times * k;
titles[slice] = "k "+k+" t "+(t-rT);
for (int row = 0; row < height; row++) {
System.arraycopy(
kernels[k][t][row],
0,
data[slice],
width*row,
width);
}
}
}
ShowDoubleFloatArrays.showArraysHyperstack(
data, // double[][] pixels,
width, // int width,
height, // int height,
num_times, // int slices,
title, // String title,
titles, // String [] titles,
true); // boolean show)
}
public static double [][] convolve3d( 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 double [][] data_in, // fill NaN before running, will treat NaN as 0 on input, skip if center is NaN
final int width, final int width,
final double [][][] kernel, // (odd)x(odd)x(odd} final double [][][] kernel, // (odd)x(odd)x(odd}
final boolean repeat_border){ final boolean repeat_border,
final boolean full_convolves_only){
final int height = data_in[0].length/width; final int height = data_in[0].length/width;
final int rT = (kernel.length-1)/2; final int rT = (kernel.length-1)/2;
final int rY = (kernel[0].length-1)/2; final int rY = (kernel[0].length-1)/2;
...@@ -756,6 +1760,9 @@ public class OrangeTest { ...@@ -756,6 +1760,9 @@ public class OrangeTest {
boolean limited = (t0 < 0) || (t1 >= lenT) || (y0 < 0) || (y1 >= height) || (x0 < 0) || (x1 >= width); boolean limited = (t0 < 0) || (t1 >= lenT) || (y0 < 0) || (y1 >= height) || (x0 < 0) || (x1 >= width);
if (limited) { if (limited) {
if (!repeat_border) { if (!repeat_border) {
if (full_convolves_only) {
continue;
}
if (t0 < 0) { if (t0 < 0) {
kTindx += t0; // negative kTindx += t0; // negative
t0 = 0; t0 = 0;
...@@ -772,48 +1779,56 @@ public class OrangeTest { ...@@ -772,48 +1779,56 @@ public class OrangeTest {
} }
if (x1 >= width) x1 = width-1; if (x1 >= width) x1 = width-1;
} }
int kTi = kTindx; no_nan_limited: {
double s = 0; int kTi = kTindx;
for (int tt = t0; tt <= t1; tt++) { double s = 0;
int t = Math.min(lenTm1, Math.max(tt, 0)); // if repeating for (int tt = t0; tt <= t1; tt++) {
int kYi = kYindx; int t = Math.min(lenTm1, Math.max(tt, 0)); // if repeating
for (int yy = y0; yy <= y1; yy++ ) { int kYi = kYindx;
int y = Math.min(lenYm1, Math.max(yy, 0)); // if repeating for (int yy = y0; yy <= y1; yy++ ) {
int kXi = kXindx; int y = Math.min(lenYm1, Math.max(yy, 0)); // if repeating
for (int xx = x0; xx <= x1; xx++ ) { int kXi = kXindx;
int x = Math.min(lenXm1, Math.max(xx, 0)); // if repeating for (int xx = x0; xx <= x1; xx++ ) {
int xy = x + width * y; int x = Math.min(lenXm1, Math.max(xx, 0)); // if repeating
double d = data_in[t][xy]; int xy = x + width * y;
if (!Double.isNaN(d)) { double d = data_in[t][xy];
s += d*kernel[kTi][kYi][kXi]; if (!Double.isNaN(d)) {
} s += d*kernel[kTi][kYi][kXi];
kXi--; } else if (full_convolves_only) {
} break no_nan_limited;
kYi--; }
kXi--;
}
kYi--;
}
kTi--;
} }
kTi--; data_out[iT][iXY] = s;
} }
data_out[iT][iXY] = s;
} else { // same faster, w/o tests } else { // same faster, w/o tests
int kTi = kTindx; no_nan: {
double s = 0; int kTi = kTindx;
for (int t = t0; t <= t1; t++) { double s = 0;
int kYi = kYindx; for (int t = t0; t <= t1; t++) {
for (int y = y0; y <= y1; y++ ) { int kYi = kYindx;
int kXi = kXindx; for (int y = y0; y <= y1; y++ ) {
for (int x = x0; x <= x1; x++ ) { int kXi = kXindx;
int xy = x + width * y; for (int x = x0; x <= x1; x++ ) {
double d = data_in[t][xy]; int xy = x + width * y;
if (!Double.isNaN(d)) { double d = data_in[t][xy];
s += d*kernel[kTi][kYi][kXi]; if (!Double.isNaN(d)) {
} s += d*kernel[kTi][kYi][kXi];
kXi--; } else if (full_convolves_only) {
} break no_nan;
kYi--; }
kXi--;
}
kYi--;
}
kTi--;
} }
kTi--; data_out[iT][iXY] = s;
} }
data_out[iT][iXY] = s;
} }
} }
} }
...@@ -824,6 +1839,33 @@ public class OrangeTest { ...@@ -824,6 +1839,33 @@ public class OrangeTest {
return data_out; return data_out;
} }
public static double [][][] kernelBlur3d(
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 = 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 [][][] kernelDiffTime( public static double [][][] kernelDiffTime(
double rt, double rt,
double rxy) { double rxy) {
...@@ -849,8 +1891,19 @@ public class OrangeTest { ...@@ -849,8 +1891,19 @@ public class OrangeTest {
return norm_kernel; return norm_kernel;
} }
/**
* Adding more dirs:
* 0-3 - tilts (0:+y, 1:+x+y, 2:+x, 3:+x-y),
* 4-7 - lines (4:along y, 5: along x+y, 6:along x, 7: along x-y)
* 8 - center,
* @param dir
* @param rt
* @param rxy
* @return
*/
public static double [][][] kernelXY( 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)
int dir, // 0: +y, 1: x+y, 2:x, 3 +x-y, 4: 0: +center will be voting 1 of 10 (each plus and minus), weighted by abs(kernelDiffTime)
double rt, double rt,
double rxy) { double rxy) {
double rxy2 = rxy*rxy; double rxy2 = rxy*rxy;
...@@ -858,6 +1911,7 @@ public class OrangeTest { ...@@ -858,6 +1911,7 @@ public class OrangeTest {
int irxy = (int) Math.floor(rxy); int irxy = (int) Math.floor(rxy);
int lt = 2 * irt + 1; int lt = 2 * irt + 1;
int lxy = 2 * irxy + 1; int lxy = 2 * irxy + 1;
double sqrt05 = Math.sqrt(0.5);
double [][][] kernel = new double [lt][lxy][lxy]; double [][][] kernel = new double [lt][lxy][lxy];
for (int it = 0; it < lt; it++) { for (int it = 0; it < lt; it++) {
double sum_cos = 0; double sum_cos = 0;
...@@ -882,7 +1936,11 @@ public class OrangeTest { ...@@ -882,7 +1936,11 @@ public class OrangeTest {
case 1: w *= y - x; break; case 1: w *= y - x; break;
case 2: w *= -x; break; case 2: w *= -x; break;
case 3: w *= -y - x; break; case 3: w *= -y - x; break;
case 4: w *= Math.cos(Math.PI * r); break; // center case 4: w *= Math.cos(Math.PI * x / rxy); break;
case 5: w *= Math.cos(Math.PI * sqrt05 * (x - y) / rxy); break;
case 6: w *= Math.cos(Math.PI * y / rxy); break;
case 7: w *= Math.cos(Math.PI * sqrt05 * (x + y) / rxy); break;
case 8: w *= Math.cos(Math.PI * r); break; // center
} }
sum_patt+=w; sum_patt+=w;
} }
...@@ -929,11 +1987,18 @@ public class OrangeTest { ...@@ -929,11 +1987,18 @@ public class OrangeTest {
return kernel_out; return kernel_out;
} }
public static double [] applyWeights( public static double [][] applyWeights(
final double [][] data_in, final double [][] data_in,
final double [][] weights) { final double [][] weights) {
final double [] data_out = new double [data_in[0].length]; final int num_frames = data_in.length;
final double [] sum_weights = new double [data_in[0].length]; final int num_pixels = data_in[0].length;
if ((weights.length != num_frames) || (weights[0].length != num_pixels)) {
throw new IllegalArgumentException("Weights do not match data");
}
final double [] data_out = new double [num_pixels];
final double [] sum_weights = new double [num_pixels];
final Thread[] threads = ImageDtt.newThreadArray(); final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0); final AtomicInteger ai = new AtomicInteger(0);
for (int nimg = 0; nimg < weights.length; nimg++) { for (int nimg = 0; nimg < weights.length; nimg++) {
...@@ -973,13 +2038,15 @@ public class OrangeTest { ...@@ -973,13 +2038,15 @@ public class OrangeTest {
}; };
} }
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
return data_out; return new double [][] { data_out, sum_weights};
} }
////
public static void testOrange() { public static void testOrange0() {
// double shrink= 8.0; // multiply by oscale * 2
// double fade_sigma = 4.0;
boolean zero_is_nan=true; boolean zero_is_nan=true;
int master = 3; int master = 3;
final double [] voffsets = value_offsets; final double [] voffsets = value_offsets;
...@@ -1014,29 +2081,31 @@ public class OrangeTest { ...@@ -1014,29 +2081,31 @@ public class OrangeTest {
final Thread[] threads = ImageDtt.newThreadArray(); final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0); final AtomicInteger ai = new AtomicInteger(0);
/* // final AtomicInteger anans = new AtomicInteger(0);
faded_textures[nslice][0] = TileProcessor.fillNaNs( // final boolean [] nan_pixels = new boolean [olen];
faded_textures[nslice][0], // final double [] data, // final double [] weights = new double [olen];
null, // final boolean [] prohibit, // final double [] sum_weights = new double [olen];
width, // int width, // final TileNeibs tn = new TileNeibs (owoi.width, owoi.height);
// CAREFUL ! Remaining NaN is grown by unsharp mask filter ************* ! // final int grow_nans = (int) Math.round(shrink * oscale * 2);
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++) { for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
// Arrays.fill(nan_pixels, false);
final double voffs = +voffsets[nimg]; final double voffs = +voffsets[nimg];
final int foffs = seg_offs_len[nimg][0]; final int foffs = seg_offs_len[nimg][0];
ImagePlus imp_src = new ImagePlus(src_files[nimg]); ImagePlus imp_src = new ImagePlus(src_files[nimg]);
ImageStack stack = imp_src.getImageStack(); ImageStack stack = imp_src.getImageStack();
final int src_width = imp_src.getWidth(); final int src_width = imp_src.getWidth();
final int src_height = imp_src.getHeight(); final int src_height = imp_src.getHeight();
String [] stack_labels = stack.getSliceLabels(); final int num_slices = stack.getSize(); // .length;
final String [] stack_labels = new String[num_slices];
System.arraycopy(
stack.getSliceLabels(),
0,
stack_labels,
0,
num_slices);
final double [][] dpixels = new double [seg_offs_len[nimg][1]][]; final double [][] dpixels = new double [seg_offs_len[nimg][1]][];
for (int i = 0; i < seg_offs_len[nimg][1]; i++) { 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]; titles[seg_offs_len[nimg][0] + i] =nimg+":"+stack_labels[segments[nimg][0]+i];
...@@ -1048,7 +2117,6 @@ public class OrangeTest { ...@@ -1048,7 +2117,6 @@ public class OrangeTest {
dpixels[i][j] = Double.NaN; 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 [][] points0 = new double[matching_points[master].length][];
final double [][] points1 = new double[points0.length][]; // matching_points[nimg]; final double [][] points1 = new double[points0.length][]; // matching_points[nimg];
...@@ -1064,20 +2132,17 @@ public class OrangeTest { ...@@ -1064,20 +2132,17 @@ public class OrangeTest {
points1); // double [][] y){ points1); // double [][] y){
double k = 1.0/oscale; 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 [][] 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 = { final double [] b1 = {
ab1[2][0] + a1[0][0]*owoi.x + a1[0][1]*owoi.y, 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 ab1[2][1] + a1[1][0]*owoi.x + a1[1][1]*owoi.y}; // from output coordinates to image2 coordinates
ai.set(0); ai.set(0);
// anans.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) { for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() { threads[ithread] = new Thread() {
public void run() { public void run() {
for (int ipix = ai.getAndIncrement(); ipix < olen; ipix = ai.getAndIncrement()) { for (int ipix = ai.getAndIncrement(); ipix < olen; ipix = ai.getAndIncrement()) {
int y = ipix / owidth; // output image coordinates int y = ipix / owidth; // output image coordinates
int x = 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 fx = a1[0][0] * x + a1[0][1] * y + b1[0];
double fy = a1[1][0] * x + a1[1][1] * y + b1[1]; double fy = a1[1][0] * x + a1[1][1] * y + b1[1];
int ix0 = (int) Math.floor(fx); int ix0 = (int) Math.floor(fx);
...@@ -1106,6 +2171,9 @@ public class OrangeTest { ...@@ -1106,6 +2171,9 @@ public class OrangeTest {
if (!Double.isNaN(d)) { if (!Double.isNaN(d)) {
out_data[indx_avg][ipix] += d; out_data[indx_avg][ipix] += d;
num_valid[ipix]++; num_valid[ipix]++;
} else {
// nan_pixels[ipix] = true;
// anans.getAndIncrement();
} }
} }
} }
...@@ -1114,6 +2182,38 @@ public class OrangeTest { ...@@ -1114,6 +2182,38 @@ public class OrangeTest {
}; };
} }
ImageDtt.startAndJoin(threads); ImageDtt.startAndJoin(threads);
/*
if (anans.get() == 0) {
tn.growSelection(
grow_nans, // int grow, // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
nan_pixels, // final boolean [] tiles,
null); // final boolean [] prohibit)
for (int i = 0; i < nan_pixels.length; i++) {
weights[i] = nan_pixels[i]? 0.0 : 1.0;
}
(new DoubleGaussianBlur()).blurDouble(
weights, // double[] pixels,
tn.getSizeX(), // int width,
tn.getSizeY(), // int height,
fade_sigma, // double sigmaX,
fade_sigma, // double sigmaY,
0.01); // double accuracy);
} else {
Arrays.fill(weights, 1.0);
}
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()) {
}
}
};
}
ImageDtt.startAndJoin(threads);
*/
} }
for (int ipix=0; ipix < olen; ipix++) { for (int ipix=0; ipix < olen; ipix++) {
...@@ -1133,5 +2233,334 @@ public class OrangeTest { ...@@ -1133,5 +2233,334 @@ public class OrangeTest {
return; return;
} }
public static double [][] getFadingWeights(
final double [][] data, // input data from one file, with NaNs
final int width,
final int shrink, // 2 for 8 directions
final double fade_sigma, // gaussian blur
final boolean nan_only){
final int hshrink = shrink/2;
final int num_scenes = data.length;
final int num_pixels = data[0].length;
final int height = num_pixels/width;
final double [][] weights = new double [num_scenes][num_pixels];
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() {
boolean [] nan_pixels = new boolean [num_pixels];
final TileNeibs tn = new TileNeibs (width, height);
for (int nscene = ai.getAndIncrement(); nscene < num_scenes; nscene = ai.getAndIncrement()) {
int num_nans = 0;
Arrays.fill(nan_pixels, false);
Arrays.fill(weights[nscene], 1.0);
for (int ipix = 0; ipix < num_pixels; ipix++) if (Double.isNaN(data[nscene][ipix])){
nan_pixels[ipix] = true; // invert
num_nans++;
}
if (num_nans > 0) {
// grow here, before borders
tn.growSelection(
shrink, // int grow, // grow tile selection by 1 over non-background tiles 1: 4 directions, 2 - 8 directions, 3 - 8 by 1, 4 by 1 more
nan_pixels, // final boolean [] tiles,
null); // final boolean [] prohibit)
for (int i = 0; i < nan_pixels.length; i++) if (nan_pixels[i]){
weights[nscene][i] = 0.0;
}
}
if (!nan_only) {
for (int n = 0; n < hshrink; n++) {
int indx0 = n * width;
int indx1 = (height-n-1) * width;
int n1 = width-n-1;
for (int i = 0; i < width; i++) {
weights[nscene][indx0 + i] = 0;
weights[nscene][indx1 + i] = 0;
}
for (int i = hshrink; i < height-hshrink; i++) {
weights[nscene][i*width + n] = 0;
weights[nscene][i*width + n1] = 0;
}
}
}
(new DoubleGaussianBlur()).blurDouble(
weights[nscene], // double[] pixels,
tn.getSizeX(), // int width,
tn.getSizeY(), // int height,
fade_sigma, // double sigmaX,
fade_sigma, // double sigmaY,
0.01); // double accuracy);
}
}
};
}
ImageDtt.startAndJoin(threads);
return weights;
}
public static double [][] warpImageSequenceBilinear(
final double [][] data_in,
final double [][] weights, // may be null?
final int width,
final double [][] ab1, // {{a00, a01},{a10,a11},{b0,b1}}
final Rectangle owoi,
final double offs, // add to data out
final double oscale) {
final int num_scenes = data_in.length;
final int num_pixels = data_in[0].length;
final int height = num_pixels/width;
final int owidth = owoi.width;
final int oheight = owoi.height;
final int num_opixels = owidth*oheight;
final double [][] data_out = new double [num_scenes+2][num_opixels]; // {..., average, sum_weight}
final int indx_average = num_scenes;
final int indx_sumweights = num_scenes+1;
for (int nscene = 0; nscene <= indx_average; nscene++) {
Arrays.fill(data_out[nscene], Double.NaN); // [indx_sumweights] = 0
}
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
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 = {
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
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < num_opixels; ipix = ai.getAndIncrement()) {
int y = ipix / owidth; // output image coordinates
int x = ipix % owidth; // output image coordinates
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 < width) && (iy0 < height)) {
double dx = fx - ix0;
double dy = fy - iy0;
if (ix0 < 0) ix0 = 0;
if (iy0 < 0) iy0 = 0;
if (ix1 >= width) ix1 = width - 1;
if (iy1 >= height) iy1 = height - 1;
int indx00 = ix0 + iy0 * width;
int indx01 = ix1 + iy0 * width;
int indx10 = ix0 + iy1 * width;
int indx11 = ix1 + iy1 * width;
for (int nscene = 0; nscene < num_scenes; nscene++) {
double d00 = data_in[nscene][indx00];
double d01 = data_in[nscene][indx01];
double d10 = data_in[nscene][indx10];
double d11 = data_in[nscene][indx11];
double w00 = weights[nscene][indx00];
double w01 = weights[nscene][indx01];
double w10 = weights[nscene][indx10];
double w11 = weights[nscene][indx11];
double d = d00*(1-dx)*(1-dy)+d01*dx*(1-dy)+d10*(1-dx)*dy+d11*dx*dy + offs;
double w = w00*(1-dx)*(1-dy)+w01*dx*(1-dy)+w10*(1-dx)*dy+w11*dx*dy;
data_out[nscene][ipix]= d;
if (!Double.isNaN(d) && (w > 0)) {
double d_old = data_out[indx_average][ipix];
if (Double.isNaN(d_old)) {
data_out[indx_average][ipix] = d;
data_out[indx_sumweights][ipix] = w;
} else {
double w_old = data_out[indx_sumweights][ipix];
double w_new = w_old + w;
data_out[indx_average][ipix] = (d_old*w_old + d * w)/w_new;
data_out[indx_sumweights][ipix] = w_new;
}
}
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return data_out;
}
public static void combineAverages(
final double [] new_data,
final double [] new_weight,
final double [] avg_data,
final double [] sum_weight) {
final int num_pixels = new_data.length;
if ((new_weight.length != num_pixels) || (avg_data.length != num_pixels) || (sum_weight.length != num_pixels)) {
throw new IllegalArgumentException("Not all data arrays are the same length");
}
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 < num_pixels; ipix = ai.getAndIncrement()) {
double d = new_data[ipix];
double w = new_weight[ipix];
if (!Double.isNaN(d) && (w > 0)) {
double d_old = avg_data[ipix];
if (Double.isNaN(d_old)) {
avg_data[ipix] = d;
sum_weight[ipix] = w;
} else {
double w_old = sum_weight[ipix];
double w_new = w_old + w;
avg_data[ipix] = (d_old*w_old + d * w)/w_new;
sum_weight[ipix] = w_new;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
}
public static void testOrange() {
int shrink_src2= 16; // shrink source by 8 pixels
double fade_sigma = 4.0;
final boolean nan_only = false; // shrink from borders also
boolean show_debug = true;
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 num_scenes = offs;
final int olen = owoi.width * owoi.height;
// final int owidth = owoi.width;
final double [][] out_data = new double [num_scenes+2][]; // [owoi.width * owoi.height];
final int indx_average = num_scenes;
final int indx_sumweights = num_scenes+1;
out_data[indx_average] = new double [olen];
Arrays.fill(out_data[indx_average], Double.NaN); // [indx_sumweights] = 0
out_data[indx_sumweights] = new double [olen];
String [] titles = new String[out_data.length];
titles [indx_average] = "average";
titles [indx_sumweights] = "weight";
for (int nimg = 0; nimg < seg_offs_len.length; nimg++) {
ImagePlus imp_src = new ImagePlus(src_files[nimg]);
ImageStack stack = imp_src.getImageStack();
final int src_width = imp_src.getWidth();
final int num_slices = stack.getSize(); // .length;
final String [] stack_labels = new String[num_slices];
System.arraycopy(
stack.getSliceLabels(),
0,
stack_labels,
0,
num_slices);
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;
}
}
}
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};
}
double [][] ab1 = ComboMatch.getAffineMatrix(
points0, // double [][] x,
points1); // double [][] y){
final double [][] weights_in = getFadingWeights(
dpixels, // final double [][] data, // input data from one file, with NaNs
src_width, // final int width,
shrink_src2, // final int shrink, // 2 for 8 directions
fade_sigma, // final double fade_sigma, // gaussian blur
nan_only); // final boolean nan_only)
final double [][] data_warped = warpImageSequenceBilinear(
dpixels, // final double [][] data_in,
weights_in, // final double [][] weights, // may be null?
src_width, // final int width,
ab1, // final double [][] ab1, // {{a00, a01},{a10,a11},{b0,b1}}
owoi, // final Rectangle owoi,
voffsets[nimg], // final double offs, // add to data out
oscale); // final double oscale)
if (show_debug) {
String [] dbg_titles = new String[data_warped.length];
for (int i = 0; i < seg_offs_len[nimg][1]; i++) {
dbg_titles[ i] =i+":"+stack_labels[segments[nimg][0]+i];
}
dbg_titles[dbg_titles.length-2]="average";
dbg_titles[dbg_titles.length-1]="weight";
ShowDoubleFloatArrays.showArrays(
data_warped,
owoi.width,
owoi.height,
true,
"merged_seq_frag_"+imp_src.getTitle()+"-scale"+oscale+"-slices"+data_warped.length+"_x"+woi.x+"_y"+woi.y+"_w"+woi.width+"_h"+woi.height+(zero_is_nan?"_zeronan":""),
titles);
}
// copy individual images
System.arraycopy(
data_warped,
0,
out_data,
seg_offs_len[nimg][0], // scene offset
seg_offs_len[nimg][1]); // number of scenes in a file
// combine averages and weights
combineAverages(
data_warped[data_warped.length-2], // final double [] new_data,
data_warped[data_warped.length-1], // final double [] new_weight,
out_data[indx_average], // final double [] avg_data,
out_data[indx_sumweights]); // final double [] sum_weight)
}
ShowDoubleFloatArrays.showArrays(
out_data,
owoi.width,
owoi.height,
true,
"merged_seq_scale"+oscale+"-slices"+indx_average+"_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