Skip to content

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

Added temporal unsharp mask

parent e3bbc197
Hide whitespace changes
Inline Side-by-side
Showing with 278 additions and 11 deletions
src/main/java/com/elphel/imagej/cuas/CuasMotion.java
View file @ 147f6635
......@@ -274,6 +274,8 @@ public class CuasMotion {
int corr_offset = clt_parameters.imp.cuas_corr_offset;
int precorr_ra = clt_parameters.imp.cuas_precorr_ra;
int corr_ra_step = clt_parameters.imp.cuas_corr_step;
int temporal_um = clt_parameters.imp.cuas_temporal_um;
double tum_threshold = clt_parameters.imp.cuas_tum_threshold;// 0.003; // clt_parameters.imp.rln_sngl_rstr; // FIXME: ADD
double fat_zero = clt_parameters.imp.cuas_fat_zero;
double cent_radius = clt_parameters.imp.cuas_cent_radius;
......@@ -2290,6 +2292,159 @@ public class CuasMotion {
return out_pix;
}
/**
* Version of runningAverage) that tolerates NaN value. It modifies fpixels[][] replacing first/last scene values if they were NaNs
* @param fpixels
* @param ra_length
* @param width
* @return
*/
public static float [][] runningAverageNaN(
final float [][] fpixels,
final int ra_length,
final int width){
final int num_scenes = fpixels.length;
final int num_pixels = fpixels[0].length;
final int height = num_pixels/width;
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final float [][] out_pix = new float [num_scenes][num_pixels];
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int nSeq = ai.getAndIncrement(); nSeq < num_scenes; nSeq = ai.getAndIncrement()) {
Arrays.fill(out_pix[nSeq],Float.NaN);
}
}
};
}
ImageDtt.startAndJoin(threads);
ai.set(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
float [] line_ra = new float[width];
int [] line_ra_length = new int [width];
for (int nLine = ai.getAndIncrement(); nLine < height; nLine = ai.getAndIncrement()) {
Arrays.fill(line_ra, 0.0f);
Arrays.fill(line_ra_length, 0);
int pix0 = nLine* width;
// fix first and last pixels if they ar nulls;
for (int x = 0; x < width; x++) {
if (Float.isNaN(fpixels[0][pix0+x])) {
int i;
for (i = 1; Float.isNaN(fpixels[i][pix0+x]) && (i < num_scenes); i++);
fpixels[0][pix0+x] = fpixels[i][pix0+x];
}
if (Float.isNaN(fpixels[num_scenes-1][pix0+x])) {
int i;
for (i = num_scenes-1; Float.isNaN(fpixels[i][pix0+x]) && (i >=0); i--);
fpixels[num_scenes-1][pix0+x] = fpixels[i][pix0+x];
}
}
int head = 0;
int tail = 0;
for (; (head < num_scenes) || (tail < num_scenes);) {
if (head < num_scenes) {
float [] fpixels_head = fpixels[head];
int pix = pix0;
for (int x = 0; x < width; x++) {
float d = fpixels_head[pix++];
if (!Float.isNaN(d)) {
line_ra[x] += d;
line_ra_length[x] ++;
}
}
head++;
}
if ((tail < (head - ra_length)) || (tail >= (num_scenes - ra_length))) {
float [] fpixels_tail = fpixels[tail];
int pix = pix0;
for (int x = 0; x < width; x++) {
float d = fpixels_tail[pix++];
if (!Float.isNaN(d)) {
line_ra[x] -= d;
line_ra_length[x]--;
}
}
tail++;
}
int avg_len = head - tail;
int avg_pos = (head + tail) / 2;
if ((avg_pos < num_scenes) && (((avg_len + ra_length) % 2 == 0 ) || (head == 1))) { // update output array
int pix = pix0;
for (int x = 0; x < width; x++) {
out_pix[avg_pos][pix++] = line_ra[x]/line_ra_length[x]; // avg_len;// should never be 0 as the first and the last are not NaN;
}
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return out_pix;
}
public static float [][] removeOutliers(
final float [][] fpixels,
final float [][] fpixels_avg,
double threshold){
final float fthreshold = (float) threshold;
final int num_scenes = fpixels.length;
final int num_pixels = fpixels[0].length;
float [][] fpixels_filtered = new float [num_scenes][];
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 nSeq = ai.getAndIncrement(); nSeq < num_scenes; nSeq = ai.getAndIncrement()) {
float [] fpix = fpixels[nSeq];
float [] apix = fpixels_avg[nSeq];
fpixels_filtered[nSeq]= fpix.clone();
for (int pix = 0; pix < num_pixels; pix++) {
if (Math.abs(fpix[pix] - apix[pix]) > fthreshold) {
fpixels_filtered[nSeq][pix] = Float.NaN;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
ai.set(0);
return fpixels;
}
public static float [][] subtract(
float [][] fpixels,
float [][] fpixels_avg){
final int num_scenes = fpixels.length;
final int num_pixels = fpixels[0].length;
float [][] fpixels_diff = new float [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() {
for (int nSeq = ai.getAndIncrement(); nSeq < num_scenes; nSeq = ai.getAndIncrement()) {
float [] fpix = fpixels[nSeq];
float [] apix = fpixels_avg[nSeq];
float [] dpix = fpixels_diff[nSeq];
for (int pix = 0; pix < num_pixels; pix++) {
dpix[pix] = fpix[pix] - apix[pix];
}
}
}
};
}
ImageDtt.startAndJoin(threads);
ai.set(0);
return fpixels_diff;
}
......@@ -2794,8 +2949,6 @@ public class CuasMotion {
int corr_offset = clt_parameters.imp.cuas_corr_offset;
int corr_pairs = clt_parameters.imp.cuas_corr_pairs;
boolean half_step = clt_parameters.imp.cuas_half_step; // true;
int precorr_ra = clt_parameters.imp.cuas_precorr_ra;
int corr_ra_step = clt_parameters.imp.cuas_corr_step;
int seq_length = corr_offset + corr_pairs;
int corr_inc = half_step ? (corr_offset/2) : corr_offset;
......@@ -2838,7 +2991,7 @@ public class CuasMotion {
CLTParameters clt_parameters,
final boolean batch_mode,
QuadCLT parentCLT, //
final float [][] fpixels,
float [][] fpixels,
final CuasMotion cuasMotion,
String [] scene_titles, // recreate slice_titles from scene titles?
String [] slice_titles,
......@@ -2846,6 +2999,8 @@ public class CuasMotion {
int corr_offset = clt_parameters.imp.cuas_corr_offset;
int corr_pairs = clt_parameters.imp.cuas_corr_pairs;
int precorr_ra = clt_parameters.imp.cuas_precorr_ra;
int temporal_um = clt_parameters.imp.cuas_temporal_um;
double tum_threshold = clt_parameters.imp.cuas_tum_threshold;// 0.003; // clt_parameters.imp.rln_sngl_rstr; // FIXME: ADD
int corr_ra_step = clt_parameters.imp.cuas_corr_step;
double fat_zero = clt_parameters.imp.cuas_fat_zero;
double cent_radius = clt_parameters.imp.cuas_cent_radius;
......@@ -2887,7 +3042,7 @@ public class CuasMotion {
double lma_maxr = clt_parameters.imp.cuas_lma_maxr; // = 5.0; // Maximal radius (>3.8)
double lma_mink = clt_parameters.imp.cuas_lma_mink; // = 0.0; // Minimal K (overshoot) <0.007
double lma_maxk = clt_parameters.imp.cuas_lma_maxk; // = 5.0; // Minimal K (overshoot) > 3.8
boolean remove_isolated= clt_parameters.imp.cuas_isolated;
// boolean remove_isolated= clt_parameters.imp.cuas_isolated;
double input_range = clt_parameters.imp.cuas_input_range; // 5;
int iter_show = clt_parameters.imp.cuas_iter_show; //1; // Maximal enhancement iteration to show intermediate result (0 - none)
boolean corr2d_save_show = clt_parameters.imp.cuas_2d_save_show; //true;
......@@ -2910,7 +3065,96 @@ public class CuasMotion {
final int frame0 = start_frame + seq_length/2;
final int half_accum_range = corr_pairs/2;
int [] remain = new int [num_corr_samples];
boolean debug_tum= false;
String model_prefix = parentCLT.getImageName()+getParametersSuffix(clt_parameters,null);
if (temporal_um > 0) {
if (intermed_high && debug_tum) {
ImagePlus imp_src = ShowDoubleFloatArrays.makeArrays(
fpixels, // float[][] pixels,
cuasMotion.gpu_max_width, // int width,
cuasMotion.gpu_max_height, // int height,
model_prefix+"-SOURCE", //String title,
scene_titles); //String [] titles)
imp_src.getProcessor().setMinAndMax(-input_range/2, input_range/2);
if (!batch_mode) {
imp_src.show();
}
parentCLT.saveImagePlusInModelDirectory(imp_src); // ImagePlus imp)
}
float [][] fpixels_ra = runningAverage(
fpixels, // final float [][] fpixels,
temporal_um, // final int ra_length,
cuasMotion.gpu_max_width); // final int width);
if (intermed_high && debug_tum) {
ImagePlus imp_src_ra = ShowDoubleFloatArrays.makeArrays(
fpixels_ra, // float[][] pixels,
cuasMotion.gpu_max_width, // int width,
cuasMotion.gpu_max_height, // int height,
model_prefix+"-SOURCE-RA"+temporal_um, //String title,
scene_titles); //String [] titles)
imp_src_ra.getProcessor().setMinAndMax(-input_range/2, input_range/2);
if (!batch_mode) {
imp_src_ra.show();
}
parentCLT.saveImagePlusInModelDirectory(imp_src_ra); // ImagePlus imp)
}
if (tum_threshold > 0) {
float [][] fpixels_nan = removeOutliers(
fpixels, // float [][] fpixels,
fpixels_ra, // float [][] fpixels_avg,
tum_threshold); // float threshold)
if (intermed_high && debug_tum) {
ImagePlus imp_outliers = ShowDoubleFloatArrays.makeArrays(
fpixels_nan, // float[][] pixels,
cuasMotion.gpu_max_width, // int width,
cuasMotion.gpu_max_height, // int height,
model_prefix+"-SOURCE-OUTLIERS", //String title,
scene_titles); //String [] titles)
imp_outliers.getProcessor().setMinAndMax(-input_range/2, input_range/2);
if (!batch_mode) {
imp_outliers.show();
}
parentCLT.saveImagePlusInModelDirectory(imp_outliers); // ImagePlus imp)
}
fpixels_ra = runningAverageNaN(
fpixels_nan, // final float [][] fpixels,
temporal_um, // final int ra_length,
cuasMotion.gpu_max_width); // final int width);
if (intermed_high && debug_tum) {
ImagePlus imp_src_ra = ShowDoubleFloatArrays.makeArrays(
fpixels_ra, // float[][] pixels,
cuasMotion.gpu_max_width, // int width,
cuasMotion.gpu_max_height, // int height,
model_prefix+"-SOURCE-RA-NOOUTLIERS"+temporal_um, //String title,
scene_titles); //String [] titles)
imp_src_ra.getProcessor().setMinAndMax(-input_range/2, input_range/2);
if (!batch_mode) {
imp_src_ra.show();
}
parentCLT.saveImagePlusInModelDirectory(imp_src_ra); // ImagePlus imp)
}
// modify original fpixels
fpixels = subtract(
fpixels, // float [][] fpixels,
fpixels_ra); // float [][] fpixels_avg)
if (intermed_high) {
ImagePlus imp_diff = ShowDoubleFloatArrays.makeArrays(
fpixels, // float[][] pixels,
cuasMotion.gpu_max_width, // int width,
cuasMotion.gpu_max_height, // int height,
model_prefix+"-SOURCE-DIFF", //String title,
scene_titles); //String [] titles)
imp_diff.getProcessor().setMinAndMax(-input_range/2, input_range/2);
if (!batch_mode) {
imp_diff.show();
}
parentCLT.saveImagePlusInModelDirectory(imp_diff); // ImagePlus imp)
}
}
}
float [][] fpixels_ra = fpixels;
if (precorr_ra > 1) {
fpixels_ra = runningAverage(
......@@ -3346,8 +3590,6 @@ public class CuasMotion {
final int corr_offset = clt_parameters.imp.cuas_corr_offset;
final int corr_pairs = clt_parameters.imp.cuas_corr_pairs;
final int precorr_ra = clt_parameters.imp.cuas_precorr_ra;
final int corr_ra_step = clt_parameters.imp.cuas_corr_step;
final boolean half_step = clt_parameters.imp.cuas_half_step; // true;
final boolean smooth = clt_parameters.imp.cuas_smooth; // true;
......@@ -3398,7 +3640,7 @@ public class CuasMotion {
boolean intermed_low = clt_parameters.imp.cuas_intermed_low; //true;
boolean intermed_high = clt_parameters.imp.cuas_intermed_high; //true;
boolean intermed_giga = clt_parameters.imp.cuas_intermed_giga; //false;
// boolean intermed_giga = clt_parameters.imp.cuas_intermed_giga; //false;
boolean save_mono = clt_parameters.imp.cuas_save_mono; //true;
boolean save_color = clt_parameters.imp.cuas_save_color; //true;
......@@ -3641,8 +3883,10 @@ public class CuasMotion {
double velocity_scale = 1.0/corr_offset;
double [][][] targets60hz = new double [fpixels.length][][];
float [][] background = fpixels;
float [][] background = fpixels;
String ra_bg_suffix=(ra_background? "-RABG":"");
if (ra_background) {
background = runningAverage(
......@@ -3773,15 +4017,21 @@ public class CuasMotion {
String prefix) {
int corr_offset = clt_parameters.imp.cuas_corr_offset;
int corr_pairs = clt_parameters.imp.cuas_corr_pairs;
int precorr_ra = clt_parameters.imp.cuas_precorr_ra;
int corr_ra_step = clt_parameters.imp.cuas_corr_step;
int temporal_um = clt_parameters.imp.cuas_temporal_um;
double tum_threshold = clt_parameters.imp.cuas_tum_threshold;// 0.003; // clt_parameters.imp.rln_sngl_rstr; // FIXME: ADD
int precorr_ra = clt_parameters.imp.cuas_precorr_ra;
int corr_ra_step = clt_parameters.imp.cuas_corr_step;
double rstr = clt_parameters.imp.cuas_rstr;// 0.003; // clt_parameters.imp.rln_sngl_rstr; // FIXME: ADD
double lma_rms = clt_parameters.imp.cuas_lma_rms; // = 1.5; // Maximal RMS, regardless of amplitude
double lma_arms = clt_parameters.imp.cuas_lma_arms; // = 0.06; // Maximal absolute RMS, sufficient for any amplitude
double lma_rrms = clt_parameters.imp.cuas_lma_rrms; // = 0.15; // Maximal relative to A rms. OK is when (RMS < cuas_lma_arms) || (RMS < cuas_lma_rrms * A)
double lma_mina = clt_parameters.imp.cuas_lma_mina; // = 1.0; // Minimal A (amplitude)
double fat_zero = clt_parameters.imp.cuas_fat_zero;
String ps = ((prefix != null)?prefix:"") + "-OFFS"+corr_offset+"-PAIRS"+corr_pairs+"-RSTR"+rstr+"-RMS"+lma_rms+"-SRMS"+lma_arms+"-RRMS"+lma_rrms+"-AMP"+lma_mina+"-FZ"+fat_zero+"-PCRA"+precorr_ra+"-PS"+corr_ra_step;
String ps = ((prefix != null)?prefix:"") + "-OFFS"+corr_offset+"-PAIRS"+corr_pairs+"-RSTR"+rstr+"-RMS"+lma_rms+
"-SRMS"+lma_arms+"-RRMS"+lma_rrms+"-AMP"+lma_mina+"-FZ"+fat_zero+"-PCRA"+precorr_ra+"-PS"+corr_ra_step+
"-TUM"+temporal_um+"-TT"+tum_threshold;
return ps;
}
......
src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java
View file @ 147f6635
......@@ -707,6 +707,9 @@ min_str_neib_fpn 0.35
public boolean cuas_smooth = true; // used cosine window when averaging correlations
public int cuas_corr_pairs = 50; // number of correlation pairs to accumulate
public int cuas_corr_offset = 20; // offset between motion detection pairs
public int cuas_temporal_um = 100; // temporal "unsharp mask" - subtract running
public double cuas_tum_threshold = 5.0; // if >0, remove outliers frim the running average and recalculate RA
public int cuas_precorr_ra = 10; // rolling average before correlation
public int cuas_corr_step = 5; // correlation step when using rolling average
......@@ -2210,6 +2213,10 @@ min_str_neib_fpn 0.35
"The number of correlation pairs to accumulate.");
gd.addNumericField("Pairs offset", this.cuas_corr_offset, 0,3,"scenes",
"Offset between the correlation pairs");
gd.addNumericField("Temporal unsharp mask length", this.cuas_temporal_um, 0,3,"scenes",
"Subtract running average this long.");
gd.addNumericField("Temporal UM threshold", this.cuas_tum_threshold, 5,8,"",
"Remove outlier pixels that differ from the running average by more than this (to later re-calculate running average excluding those pixels.");
gd.addNumericField("Pre-correlation running average", this.cuas_precorr_ra, 0,3,"scenes",
"Smoothing input data by running average before correlation for motion vectors calculations. Target extraction wil still use individual scenes.");
......@@ -3323,6 +3330,9 @@ min_str_neib_fpn 0.35
this.cuas_corr_pairs = (int) gd.getNextNumber();
this.cuas_corr_offset = (int) gd.getNextNumber();
this.cuas_temporal_um = (int) gd.getNextNumber();
this.cuas_tum_threshold = gd.getNextNumber();
this.cuas_precorr_ra = (int) gd.getNextNumber();
this.cuas_corr_step = (int) gd.getNextNumber();
......@@ -4285,6 +4295,8 @@ min_str_neib_fpn 0.35
properties.setProperty(prefix+"cuas_smooth", this.cuas_smooth+""); // boolean
properties.setProperty(prefix+"cuas_corr_pairs", this.cuas_corr_pairs+""); // int
properties.setProperty(prefix+"cuas_corr_offset", this.cuas_corr_offset+""); // int
properties.setProperty(prefix+"cuas_temporal_um", this.cuas_temporal_um+""); // int
properties.setProperty(prefix+"cuas_tum_threshold", this.cuas_tum_threshold+""); // double
properties.setProperty(prefix+"cuas_precorr_ra", this.cuas_precorr_ra+""); // int
properties.setProperty(prefix+"cuas_corr_step", this.cuas_corr_step+""); // int
......@@ -5215,6 +5227,8 @@ min_str_neib_fpn 0.35
if (properties.getProperty(prefix+"cuas_smooth")!=null) this.cuas_smooth=Boolean.parseBoolean(properties.getProperty(prefix+"cuas_smooth"));
if (properties.getProperty(prefix+"cuas_corr_pairs")!=null) this.cuas_corr_pairs=Integer.parseInt(properties.getProperty(prefix+"cuas_corr_pairs"));
if (properties.getProperty(prefix+"cuas_corr_offset")!=null) this.cuas_corr_offset=Integer.parseInt(properties.getProperty(prefix+"cuas_corr_offset"));
if (properties.getProperty(prefix+"cuas_temporal_um")!=null) this.cuas_temporal_um=Integer.parseInt(properties.getProperty(prefix+"cuas_temporal_um"));
if (properties.getProperty(prefix+"cuas_tum_threshold")!=null) this.cuas_tum_threshold=Double.parseDouble(properties.getProperty(prefix+"cuas_tum_threshold"));
if (properties.getProperty(prefix+"cuas_precorr_ra")!=null) this.cuas_precorr_ra=Integer.parseInt(properties.getProperty(prefix+"cuas_precorr_ra"));
if (properties.getProperty(prefix+"cuas_corr_step")!=null) this.cuas_corr_step=Integer.parseInt(properties.getProperty(prefix+"cuas_corr_step"));
......@@ -6148,6 +6162,9 @@ min_str_neib_fpn 0.35
imp.cuas_smooth = this.cuas_smooth;
imp.cuas_corr_pairs = this.cuas_corr_pairs;
imp.cuas_corr_offset = this.cuas_corr_offset;
imp.cuas_temporal_um = this.cuas_temporal_um;
imp.cuas_tum_threshold = this.cuas_tum_threshold;
imp.cuas_precorr_ra = this.cuas_precorr_ra;
imp.cuas_corr_step = this.cuas_corr_step;
imp.cuas_half_step = this.cuas_half_step;
......
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