Skip to content

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

circular deconvolution kernels

parent 86efce83
Show whitespace changes
Inline Side-by-side
Showing with 399 additions and 6 deletions
src/main/java/com/elphel/imagej/common/DoubleFHT.java
View file @ 2ebcdb9c
......@@ -25,6 +25,7 @@ package com.elphel.imagej.common;
**
*/
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import ij.IJ;
......@@ -292,6 +293,60 @@ public class DoubleFHT {
return result;
}
/**
* Invert kernel for deconvolutions. Kernel is assumed to fade near edges
* @param data source kernel, square, power of 2 sides. Data is destroyed
* @param fat_zero add to denominator to prevent noise amplification
* @param highPassSigma high-pass sigma to create frequency filter or 0
* @param lowPassSigma low-pass sigma to create frequency filter or 0
* @return inverted kernel
*/
public double [] invert (
double [] data,
double fat_zero,
double highPassSigma,
double lowPassSigma
) {
double [] filter= null;
if ((highPassSigma> 0.0) || (lowPassSigma> 0.0)) {
updateMaxN(data);
createFrequencyFilter(highPassSigma, lowPassSigma); // for repetitive calls will reuse mask
filter =this.freqMask;
}
return invert (data, fat_zero, filter);
}
/**
* Invert kernel for deconvolutions. Kernel is assumed to fade near edges
* @param data source kernel, square, power of 2 sides. Data is destroyed
* @param fat_zero add to denominator to prevent noise amplification
* @param filter low-pass/high-pass filter or null
* @return inverted kernel
*/
public double [] invert (
double [] data, // destroyed
double fat_zero,
double [] filter
) {
updateMaxN(data);
int size = (int) Math.sqrt(data.length);
double nominator_FD_value = size*size; // fill FD of the nominator. what is the best value?
double [] nominator_FD = new double[data.length];
Arrays.fill(nominator_FD, nominator_FD_value);
swapQuadrants(data);
if (!transform(data,false)) return null; // direct FHT
double [] inverted = divide(nominator_FD,data, fat_zero);
if (filter != null) {
multiplyByReal(inverted, filter);
}
transform(inverted, true); // inverse transform
swapQuadrants(inverted);
return inverted;
}
/* calculates correlation, destroys original arrays */
public double[] correlate(double[] first, double[] second, double highPassSigma, double lowPassSigma,
......@@ -2320,6 +2375,8 @@ public class DoubleFHT {
return result;
}
public double[] setReal(double[] amp) { // only first half used
double[] result = new double[maxN * maxN];
int rowMod, colMod;
......
src/main/java/com/elphel/imagej/orthomosaic/ComboMatch.java
View file @ 2ebcdb9c
......@@ -282,9 +282,20 @@ public class ComboMatch {
double reversal_rad = 0; // 26;// 22; // cut at first direction reversal after
double frac_outliers = 0.4;
int default_kernel = 1; // 0; // 1; // default kernel choice 25-> 50// 1 for simulated, 0 - for extracted
boolean only_correlate = true;
boolean only_correlate = true; // false; // true;
int sub_pattern_index = 0;
double inv_fat_zero = 300;
double inv_lpf_sigma = .05;
int inv_scale_radial = 8; // scale up radial resolution
double inv_lim_rad = 50.0; // outside all 0
double inv_trans_width = 20.0; // lim_rad-trans_width - start reducing
double inv_blur_center = 1.0; // cosine blur (in original pixels) in the center (<blur_radius) area
double inv_blur_radius =100.0; // constant blur radius
double inv_blur_rate = 0.2; // blur increase outside of blur_radius
// String pattern_dir= "/media/elphel/SSD3-4GB/lwir16-proc/ortho_videos/debug/mines/pattern_25m_zoom1/synthetic/";
String pattern_dir= "/media/elphel/NVME/lwir16-proc/ortho_videos/mines_extract/evening_50m_sept12/debug_cross/new_synth/";
String [] pattern_files={
......@@ -301,7 +312,17 @@ public class ComboMatch {
"patterns_r30.0_e5.0_ir12.0_ie5.0_is-1.0_or45.0_oe20.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r30.0_e5.0_ir12.0_ie5.0_is-1.0_or50.0_oe40.0_os-1.0_h8_w0.6_s-60.0_200x200.tif", // best
"patterns_r30.0_e5.0_ir12.0_ie5.0_is-1.0_or55.0_oe50.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r30.0_e5.0_ir12.0_ie5.0_is0.0_or50.0_oe40.0_os-1.0_h8_w0.6_s-60.0_200x200.tif"
"patterns_r30.0_e5.0_ir12.0_ie5.0_is0.0_or50.0_oe40.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"mine2_multi_zoom-1-INV_PATTERN_fo0.4_lr80.0_tw60.0_rr0.0_128x128.tif",
"mine2_multi_zoom-1-INV_RADIAL_PATTERN_fo0.4_lr80.0_tw60.0_rr0.0_128x128.tif",
"patterns_r17.0_e15.0_ir10.0_ie8.0_is0.0_or50.0_oe30.0_os-1.0_h8_w0.6_s-60.0_200x200.tif", // mine2
"patterns_r17.0_e22.0_ir10.0_ie8.0_is0.0_or50.0_oe30.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r20.0_e25.0_ir10.0_ie8.0_is0.0_or50.0_oe25.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r20.0_e25.0_ir10.0_ie8.0_is0.0_or55.0_oe25.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r20.0_e25.0_ir10.0_ie8.0_is0.0_or45.0_oe20.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r20.0_e25.0_ir10.0_ie8.0_is0.0_or40.0_oe20.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r20.0_e25.0_ir10.0_ie8.0_is0.0_or35.0_oe15.0_os-1.0_h8_w0.6_s-60.0_200x200.tif",
"patterns_r20.0_e20.0_ir10.0_ie8.0_is0.0_or35.0_oe15.0_os-1.0_h8_w0.6_s-60.0_200x200.tif"
};
......@@ -328,6 +349,17 @@ public class ComboMatch {
gdo.addCheckbox ("Only correlate", only_correlate, "do not generate circular kernel.");
gdo.addNumericField("subpattern index", sub_pattern_index,0,4,"", "0 - full pattern, 1-8 - half-pattern.");
gdo.addMessage("Radial kernel inversion");
gdo.addNumericField("Inversion fat zero", inv_fat_zero, 3,7,"", "Add to amplitude to reduce noise.");
gdo.addNumericField("LPF sigma", inv_lpf_sigma, 3,7,"", "Filter inversion result.");
gdo.addNumericField("Scale up radial resolution", inv_scale_radial,0,4,"", "Increase radial resolution for artifacts reduction.");
gdo.addNumericField("Inverted kernel max radius", inv_lim_rad, 3,7,"pix", "Crop result kernel to this radius.");
gdo.addNumericField("Fade out width", inv_trans_width, 3,7,"pix", "Fade-out width for radial crop.");
gdo.addNumericField("Center area blur", inv_blur_center, 3,7,"pix", "Blur while integrating radial profile (center high-res area).");
gdo.addNumericField("Start radial blur", inv_blur_radius, 3,7,"pix", "Center area radius, blur more outside of it.");
gdo.addNumericField("Radial blur increase rate", inv_blur_rate, 3,7,"pix/pix", "Rate of the radial blur increase outside of the central area "+
"- pixel of (biderectional) blur for each radial pixel outside of center area.");
//only_correlate
gdo.showDialog();
if (gdo.wasCanceled()) return false;
......@@ -349,6 +381,15 @@ public class ComboMatch {
only_correlate = gdo.getNextBoolean();
sub_pattern_index= (int) gdo.getNextNumber();
inv_fat_zero = gdo.getNextNumber();
inv_lpf_sigma = gdo.getNextNumber();
inv_scale_radial = (int) gdo.getNextNumber();
inv_lim_rad = gdo.getNextNumber();
inv_trans_width = gdo.getNextNumber();
inv_blur_center = gdo.getNextNumber();
inv_blur_radius = gdo.getNextNumber();
inv_blur_rate = gdo.getNextNumber();
System.out.println("Will extract objects here, image zoom_level="+zool_lev_objects);
for (ObjectLocation ol: object_list) {
System.out.println(ol.name+": "+ol.xy_meters[0]+"/"+ol.xy_meters[1]);
......@@ -459,6 +500,65 @@ public class ComboMatch {
for (int j = 0; j < output_patterns[0].length; j++) {
output_patterns[object_stack.length][j] /= object_stack.length;
}
double [][] inverted_patterns = new double [output_patterns.length][];
for (int i = 0; i < output_patterns.length; i++) {
//pattern_file
String dbg_prefix = OrthoMap.removeKnownExtension(pattern_file)+"-"+i;
inverted_patterns[i] = ObjectLocation.invertRadialKernel(
output_patterns[i], // double [] kernel,
inv_fat_zero, // double fat_zero,
inv_lpf_sigma, // double lpf_sigma,
inv_scale_radial, // int scale_radial, // scale up radial resolution
inv_lim_rad, // double lim_rad, // outside all 0
inv_trans_width, // double trans_width, // lim_rad-trans_width - start reducing
inv_blur_center, // double blur_center, // cosine blur (in original pixels) in the center (<blur_radius) area
inv_blur_radius, // double blur_radius, // constant blur radius
inv_blur_rate, // double blur_rate, // blur increase outside of blur_radius
centers, // // double [][] centers,
object_stack, // double [][] src_cuts,
dbg_prefix, // String dbg_prefix,
1); // int debugLevel)
}
ImagePlus imp_inv_patt = ShowDoubleFloatArrays.makeArrays(
inverted_patterns,
corr_size,
corr_size,
OrthoMap.removeKnownExtension(imp_sel.getTitle())+"-INV_RADIAL_PATTERN"+settings_str+"_"+corr_size+"x"+corr_size+".tiff",
output_patt_titles); // test_titles,
imp_inv_patt.show();
// select kernel and correlatde with every cut image?
// Test invert patterns
/*
boolean test_invert = debugLevel > -1000;
if (test_invert) {
double invert_fz = 0;
double [][] inverted_patterns = new double [output_patterns.length][];
for (int i = 0; i < output_patterns.length; i++) {
inverted_patterns[i] = output_patterns[i].clone();
double [] dfr = new double [corr_size*corr_size];
dfr[corr_size/2*(corr_size+1)] = corr_size*corr_size;
double [][] data_pair = new double[][] {inverted_patterns[i],dfr};
inverted_patterns[i]=OrthoMap.deconvolvePair(
data_pair,
invert_fz,
1);
}
ImagePlus imp_inv_patt = ShowDoubleFloatArrays.makeArrays(
inverted_patterns,
corr_size,
corr_size,
OrthoMap.removeKnownExtension(imp_sel.getTitle())+"-INV_PATTERN"+settings_str+"_"+corr_size+"x"+corr_size+".tiff",
output_patt_titles); // test_titles,
imp_inv_patt.show();
}
*/
ImagePlus imp_out_patt = ShowDoubleFloatArrays.makeArrays(
output_patterns,
corr_size,
......@@ -684,7 +784,7 @@ public class ComboMatch {
}
}
if (pattern_match) {
ImagePlus imp_pat_match = maps_collection.patterMatchDualWrap (
ImagePlus imp_pat_match = maps_collection.patternMatchDualWrap (
gpu_pair, // int [] indices, // null or which indices to use (normally just 2 for pairwise comparison)
affines, // double [][][] affines, // null or [indices.length][2][3]
warp); // FineXYCorr warp)
......
src/main/java/com/elphel/imagej/orthomosaic/ObjectLocation.java
View file @ 2ebcdb9c
......@@ -5,6 +5,7 @@ import java.util.Arrays;
import java.util.Collections;
import java.util.Comparator;
import com.elphel.imagej.common.DoubleFHT;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.tileprocessor.Correlation2d;
import com.elphel.imagej.tileprocessor.TileNeibs;
......@@ -197,6 +198,229 @@ public class ObjectLocation {
return xys;
}
/**
* Invert kernel in frequency domain, then condition it knowing it is radial. Conditioning
* involves converting to 1d radial (linear interpolation with optionally increased radial
* resolution and simultaneous cosine blurring. Blurring may increase for the outer areas
* of the result kernel
* @param kernel direct radial kernel from extracted objects, centered at size/2,size/2
* @param fat_zero fat zero for frequency-domain inversion
* @param lpf_sigma low-pass sigma after frequency-domain inversion
* @param scale_radial resolution
* @param lim_rad limit inverted kernel by radius
* @param trans_width transition width for radial limiting (cosine)
* @param blur_center spread data in +/- blur_radius (cosine) among radial bins
* @param blur_radius radius of the constant radial blurring, outside increase blur
* @param blur_rate rate of blur increase (pix/pix) in the peripheral areas
* @param debugLevel debug level
* @return inverted kernel, same dimension. Normalize to have convolution of derect and inverted be 1.0
* in the center?
*/
public static double [] invertRadialKernel(
double [] kernel,
double fat_zero,
double lpf_sigma,
int scale_radial, // scale up radial resolution
double lim_rad, // outside all 0
double trans_width, // lim_rad-trans_width - start reducing
double blur_center, // cosine blur (in original pixels) in the center (<blur_radius) area
double blur_radius, // constant blur radius
double blur_rate, // blur increase outside of blur_radius
// just for testing - source images and corresponding centers to convolve
double [][] centers,
double [][] src_cuts,
String dbg_prefix,
int debugLevel
) {
int size = (int) Math.sqrt(kernel.length);
double [] kernel0 = kernel.clone();
double [] inverted0 = (new DoubleFHT()).invert (
kernel0, // double [] data,
fat_zero, // double fat_zero,
0, // double highPassSigma,
lpf_sigma); // double lowPassSigma);
int irad = (int) Math.ceil(scale_radial*lim_rad) +2;
double [] sw = new double [irad], swd = new double[irad];
// double maxrad2 = irad*irad;
double xc = size/2, yc=size/2;
for (int y = 0; y < size; y++) {
double y2 = (y - yc)*(y-yc);
for (int x = 0; x < size; x++) {
double r = Math.sqrt(y2+ (x-xc)*(x-xc));
double blur = (r < blur_radius) ? blur_center : ((r - blur_center) * blur_rate);
double low_r = (r - blur) * scale_radial;
if (((int) low_r) < sw.length) {
double d = inverted0[y*size + x];
double high_r = (r + blur) * scale_radial;
double s = Math.PI/(high_r-low_r);
int ilr = Math.max(0, (int) Math.ceil(low_r));
int ihr = Math.min(irad-1, (int) Math.floor(high_r));
for (int i = ilr; i <= ihr; i++) {
double w = Math.sin(s * (i - low_r));
sw[i] += w;
swd[i] += w * d;
}
}
}
}
for (int i = 0; i < sw.length; i++) {
if (sw[i] > 0) {
swd[i] /= sw[i];
}
}
// limit by lim_rad, trans_width
if (trans_width > 0) {
double lim_rad_scaled = lim_rad * scale_radial;
double trans_width_scaled = trans_width * scale_radial;
for (int i = (int) (lim_rad_scaled-trans_width_scaled); i < sw.length; i++) {
if (i < lim_rad_scaled) {
swd[i] *= 0.5*(1.0 + Math.cos(Math.PI*(i - (lim_rad_scaled-trans_width_scaled))/trans_width_scaled));
} else {
swd[i] = 0;
}
}
}
double [] inverted = new double [size*size];
double lr2 = lim_rad * lim_rad;
// generate output pattern
int oc = size/2;
for (int y = 0; y < size; y++) {
double dy = y - oc;
double y2 = dy*dy;
if (y2 <= lr2) {
for (int x = 0; x < size; x++) {
double dx = x - oc;
double r2 = y2 + dx*dx;
if (r2 <= lr2) {
double sr = Math.sqrt(r2)* scale_radial;
int isr = (int) Math.floor(sr);
double w1 = sr-isr;
double w0 = 1.0 - w1;
inverted[x + size * y] = w0 * swd[isr] + w1 * swd[isr + 1];
}
}
}
}
double [] convolved = (new DoubleFHT()).convolve(
kernel.clone(), // double[] first,
inverted.clone()); // double[] second)
double center_data = convolved[size*(size+1)/2];
double [] normalized = inverted.clone();
double k = 1.0/center_data;
for (int i = 0; i < normalized.length; i++) {
normalized[i]*=k;
}
if (debugLevel > -4) {
System.out.println("invertRadialKernel(): direct*inverted = "+center_data);
if (debugLevel >0) {
String dbg_title = dbg_prefix + "_invertRK_fz"+fat_zero+
"_lpf"+lpf_sigma+"_sr"+scale_radial+"_lr"+lim_rad+
"_tw"+trans_width+"_bc"+blur_center+"_br"+blur_radius+
"_brt"+blur_rate;
String [] dbg_titles= {"direct","inverted","radial","convolved","normalized"};
double [][] dbg_img = {kernel,inverted0,inverted,convolved,normalized};
ImagePlus imp_dbg = ShowDoubleFloatArrays.makeArrays(
dbg_img,
size,
size,
dbg_title,
dbg_titles); // test_titles,
// PointRoi roi = new PointRoi();
// roi.setOptions("label");
imp_dbg.show();
}
}
/*
public double[] convolve(double[] first, double[] second) {
return convolve(first, second, null);
}
*/
// testing by convolving with the source images
if (debugLevel>0) {
double [] wnd1d = new double[size];
boolean sq= false;
for (int i= 0; i < size; i++) {
wnd1d[i] = Math.sin(Math.PI*i/size);
if (sq) {
wnd1d[i]*=wnd1d[i];
}
}
double [] wnd = new double[size*size];
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
wnd[i*size+j] = wnd1d[i]*wnd1d[j];
}
}
int extr_size = (int) Math.sqrt(src_cuts[0].length);
double [][] test_src = new double [src_cuts.length][size*size];
for (int i = 0; i < src_cuts.length; i++) {
double xcent = extr_size/2+centers[i][0];
double ycent = extr_size/2+centers[i][1];
int x0 = (int) Math.round(xcent - size/2);
int y0 = (int) Math.round(ycent - size/2);
for (int y = 0; y < size; y++) {
int ys = y0+y;
if ((ys >= 0) && (ys < extr_size)) {
for (int x = 0; x < size; x++) {
int xs = x0+x;
if ((xs >= 0) && (xs < extr_size)) {
int indx = y*size + x;
test_src[i][indx] = src_cuts[i][ys * extr_size + xs];
}
}
}
}
OrthoMap.removeDC(test_src[i]);
for (int j = 0; j < test_src[i].length; j++) {
test_src[i][j] *= wnd[j];
}
}
double [][] test_conv = new double [src_cuts.length][];
for (int i = 0; i < src_cuts.length; i++) {
test_conv[i] = (new DoubleFHT()).convolve(
test_src[i].clone(), // double[] first,
normalized.clone()); // double[] second)
}
String dbg_title = dbg_prefix + "_conv_src_+fz"+fat_zero+
"_lpf"+lpf_sigma+"_sr"+scale_radial+"_lr"+lim_rad+
"_tw"+trans_width+"_bc"+blur_center+"_br"+blur_radius+
"_brt"+blur_rate;
String [] dbg_titles= new String[src_cuts.length] ; // {"direct","inverted","radial","convolved","normalized"};
for (int i = 0; i < src_cuts.length; i++) {
dbg_titles[i] = "img-"+i;
}
// double [][] dbg_img = {kernel,inverted0,inverted,convolved,normalized};
ImagePlus imp_dbg = ShowDoubleFloatArrays.makeArrays(
test_conv,
size,
size,
dbg_title,
dbg_titles); // test_titles,
// PointRoi roi = new PointRoi();
// roi.setOptions("label");
imp_dbg.show();
ImagePlus imp_dbg_src = ShowDoubleFloatArrays.makeArrays(
test_src,
size,
size,
dbg_prefix+"-src_cut",
dbg_titles); // test_titles,
// PointRoi roi = new PointRoi();
// roi.setOptions("label");
imp_dbg_src.show();
}
return normalized;
}
// Needs improvement. Probably look at neighbors when deleting outliers. Multipass and pull to neighbors?
public static double [] getRadialPattern(
final double [] data,
double [] xy_offs,
......
src/main/java/com/elphel/imagej/orthomosaic/OrthoMap.java
View file @ 2ebcdb9c
......@@ -1820,6 +1820,8 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
}
return kernel;
}
public static double [] deconvolvePair(
double [][] data,
double fat_zero, // 1000
......@@ -1842,6 +1844,14 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
doubleFHT.swapQuadrants(data[1]);
if (!doubleFHT.transform(data[1],false)) return null; // direct FHT
if (!doubleFHT.transform(data[0],false)) return null; // direct FHT {
if (debugLevel > 0) {
ShowDoubleFloatArrays.showArrays( // here every element
data[1],
fft_size,
fft_size,
"data[1]");
}
if (debugLevel > 0) {
double [] amp = doubleFHT.calculateAmplitude(data[0]);
ShowDoubleFloatArrays.showArrays(
......@@ -1850,7 +1860,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
fft_size,
"amp");
}
double [] deconv = doubleFHT.divide(data[1],data[0], fft_size);
double [] deconv = doubleFHT.divide(data[1],data[0], fat_zero);
if (debugLevel > 0) {
ShowDoubleFloatArrays.showArrays(
deconv,
......@@ -1860,7 +1870,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
}
// if (zero_phase) {
double [] div_amp = doubleFHT.calculateAmplitudeNoSwap(deconv);
double [] div_phase = doubleFHT.calculatePhaseNoSwap(deconv);
double [] div_phase = DoubleFHT.calculatePhaseNoSwap(deconv);
double [][] amp_phase = {div_amp, div_phase};
ShowDoubleFloatArrays.showArrays(
amp_phase,
......@@ -1903,6 +1913,8 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
}
public static void removeDC(
final double [] data) {
final Thread[] threads = ImageDtt.newThreadArray(QuadCLT.THREADS_MAX);
......
src/main/java/com/elphel/imagej/orthomosaic/OrthoMapsCollection.java
View file @ 2ebcdb9c
......@@ -1650,7 +1650,7 @@ public class OrthoMapsCollection implements Serializable{
return vf_out;
}
public ImagePlus patterMatchDualWrap (
public ImagePlus patternMatchDualWrap (
int [] indices, // null or which indices to use (normally just 2 for pairwise comparison)
double [][][] affines, // null or [indices.length][2][3]
FineXYCorr warp) { // use for a single pair only
......
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