Skip to content

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

Pattern matching

parent 81c431ad
Hide whitespace changes
Inline Side-by-side
Showing with 2629 additions and 2166 deletions
src/main/java/com/elphel/imagej/common/DoubleFHT.java
View file @ 02845a72
This source diff could not be displayed because it is too large. You can view the blob instead.
src/main/java/com/elphel/imagej/orthomosaic/OrthoMap.java
View file @ 02845a72
......@@ -2823,6 +2823,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
final int width,
final int psize, // power of 2, such as 64
final double [] pattern, // [psize*psize]
final boolean convolve, // convolve, not correlate
final double phaseCoeff,
final int debugLevel) {
final int height = data.length/width;
......@@ -2849,20 +2850,22 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
final double [] filter = null; // probably not needed
if ((width== psize) && (height == psize)) {
testPhaseCorr(
data, // double [] data,
pattern, // double [] pattern,
window, // double [] wnd,
return testPhaseCorr(
data, // double [] data,
pattern, // double [] pattern,
convolve, //final boolean convolve, // convolve, not correlate
window, // double [] wnd,
phaseCoeff); // double phaseCoeff)
return data;
}
final int tilesX = (int) Math.ceil(width/(psize/2)) + 1;
final int tilesY = (int) Math.ceil(height/(psize/2)) + 1;
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final int dbg_tileX = -20;
final int dbg_tileY = -1;
// final int dbg_tileX = 20;
// final int dbg_tileY = 1;
final int dbg_x = -2668;
final int dbg_y = 256;
for (int passY = 0; passY < 2; passY++) { // to isolate threads results - no overlap
final int tileY0=passY;
final int ntilesY=(tilesY + 1 - passY)/ 2;
......@@ -2880,15 +2883,13 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
for (int nTile = ai.getAndIncrement(); nTile < ntiles; nTile = ai.getAndIncrement()){
int tileX = tileX0 - 1 + 2 * (nTile % ntilesX); // nTile % ntilesX;
int tileY = tileY0 - 1 + 2 * (nTile / ntilesX); // nTile / ntilesX;
boolean dbg_tile = (Math.abs(tileX - dbg_tileX) < 1) && (Math.abs(tileY - dbg_tileY) < 1);
// int px0 = (psize/2) *(tileX0 - 1 + 2 * tileX); // absolute in the original/result image
// int py0 = (psize/2) *(tileY0 - 1 + 2 * tileY);
int px0 = (psize/2) * tileX; // absolute in the original/result image
int py0 = (psize/2) * tileY;
int x0 = Math.max(0, -px0);
int y0 = Math.max(0, -py0);
int x1 = Math.min(psize, width- px0);
int y1 = Math.min(psize, height-py0);
boolean dbg_tile = (Math.abs((px0 + psize/2) - dbg_x) < psize/4) && (Math.abs((py0 + psize/2) - dbg_y) < psize/4);
if (dbg_tile) {
System.out.println("correlateWithPattern(): tileX="+tileX+", tileY="+tileY);
System.out.println("correlateWithPattern(): px0="+px0+", py0="+py0);
......@@ -2929,14 +2930,22 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
"input_corr_tx"+tileX+"_ty"+tileY,
rslt_titles);
}
// now cross-correlate with fat zero dtile and
// important to assign: dtile is modified, but not to the result!
dtile = doubleFHT.phaseCorrelatePattern ( // new
dtile, // double [] first,
patternFD, // double [] secondFD,
phaseCoeff, // double phaseCoeff,
filter, // double [] filter, // high/low pass filtering
null); // double [] first_save ) //null-OK
if (convolve) { // convolve, not correlate
dtile = doubleFHT.convolvePattern (
dtile, // double [] first,
patternFD, // double [] secondFD,
filter, // double [] filter, // high/low pass filtering
null); // double [] first_save ) //null-OK
} else {
// now cross-correlate with fat zero dtile and
// important to assign: dtile is modified, but not to the result!
dtile = doubleFHT.phaseCorrelatePattern ( // new
dtile, // double [] first,
patternFD, // double [] secondFD,
phaseCoeff, // double phaseCoeff,
filter, // double [] filter, // high/low pass filtering
null); // double [] first_save ) //null-OK
}
if (dbg_tile) {
String [] rslt_titles= {"corr","pattern"};
ShowDoubleFloatArrays.showArrays(
......@@ -2965,12 +2974,13 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
return dout;
}
public static void testPhaseCorr(
public static double [] testPhaseCorr(
double [] data_in,
double [] pattern_in,
boolean convolve, // convolve, not correlate
double [] wnd_in,
double phaseCoeff) {
phaseCoeff = 0.5;
// phaseCoeff = 0.5;
double [] data = data_in.clone();
double [] pattern = pattern_in.clone();
double [] wnd = wnd_in.clone();
......@@ -2987,7 +2997,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
size,
size,
true,
"input_corr",
convolve?"input_conv":"input_corr",
rslt_titles);
}
DoubleFHT doubleFHT = new DoubleFHT();
......@@ -2998,20 +3008,35 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
if (use_new) {
final double [] patternFD = pattern.clone();
doubleFHT.transformPattern(patternFD);
corr_out = doubleFHT.phaseCorrelatePattern ( // new
if (convolve) {
corr_out = doubleFHT.convolvePattern (
data, // double [] first,
patternFD, // double [] secondFD,
null, // double [] filter, // high/low pass filtering
null); // double [] first_save ) //null-OK
} else {
corr_out = doubleFHT.phaseCorrelatePattern ( // new
data, // double [] first,
patternFD, // double [] secondFD,
phaseCoeff, // double phaseCoeff,
null, // double [] filter, // high/low pass filtering
null); // double [] first_save ) //null-OK
}
} else {
corr_out = doubleFHT. phaseCorrelate ( // new
if (convolve) {
corr_out = doubleFHT. convolve ( // new
data, // double [] first,
pattern, // double [] second,
null);// double [] filter, // high/low pass filtering
} else {
corr_out = doubleFHT. phaseCorrelate ( // new
data, // double [] first,
pattern, // double [] second,
phaseCoeff, // double phaseCoeff,
null, // double [] filter, // high/low pass filtering
null, // double [] first_save,
null); // double [] second_save ) //null-OK
null); // double [] second_save ) //null-OK
}
}
{
......@@ -3021,10 +3046,11 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
size,
size,
true,
"output_corr"+phaseCoeff,
convolve?"output_conv":("output_corr"+phaseCoeff),
rslt_titles);
}
System.out.println("testPhaseCorr() done");
return corr_out;
}
public static void fillNaNs(
......@@ -3115,12 +3141,19 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
int corr_size = 256;
double phaseCoeff = 0.5;
int debugLevel = 1;
double min_corr = 0.0025; // real max >0.005 with scale -500000;
double adv_radius = 60.0; // pixels
int corr_radius = 30; // may use same as adv_radius; but not more
GenericJTabbedDialog gd = new GenericJTabbedDialog("Correlate image with pattern",1200,500);
gd.addStringField ("Pattern directory", pattern_dir, 180, "Absolute path including trailing \"/\".");
gd.addStringField ("Pattern filename", pattern_file,180, "Pattern filename.");
gd.addNumericField("Zoom-out factor", zoomout, 0,4,"x", "Reduce pattern resolution to match image.");
gd.addNumericField("Phase correlation coefficient", phaseCoeff, 3,7,"","1.0 - pure phase correlation, 0.0 - regular correlation.");
gd.addNumericField("Minimal correlation", min_corr, 3,7,"","Minimal correlation value to keep.");
gd.addNumericField("Adversarial distance", adv_radius, 1,6,"pix", "Suppress weaker if they have closer strong ones.");
gd.addNumericField("Correlation peak max radius", corr_radius, 0,4,"pix", "Limit correlation peak radius. Should be <= adversarial");
gd.addNumericField("Correlation size", corr_size, 0,4,"pix", "Should be power of 2.");
gd.addNumericField("Debug level", debugLevel, 0,4,"", "Debug level.");
gd.showDialog();
if (gd.wasCanceled()) return;
......@@ -3128,6 +3161,10 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
pattern_file= gd.getNextString();
zoomout= (int) gd.getNextNumber();
phaseCoeff= gd.getNextNumber();
min_corr= gd.getNextNumber();
adv_radius= gd.getNextNumber();
corr_radius= (int) gd.getNextNumber();
corr_size= (int) gd.getNextNumber();
debugLevel= (int) gd.getNextNumber();
float [] fpixels = (float[]) imp_src.getProcessor().getPixels();
......@@ -3160,44 +3197,103 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
patterns[n][i] = fpixels_pattern[i];
}
}
/*
float [] fpixels_pattern = (float[]) imp_pattern.getProcessor().getPixels();
double [ ] pattern = new double[fpixels_pattern.length];
for (int i = 0; i < fpixels_pattern.length; i++) {
pattern[i] = fpixels_pattern[i];
}
*/
double [][] corrs_out = new double[patterns.length][];
double [][] corrs_out = new double[patterns.length][];
double [][] convolve_out = new double[patterns.length][];
double [][] lim_corr = new double[patterns.length][];
double [][] corr_patterns = new double[patterns.length][];
for (int n = 0; n < patterns.length; n++) {
double [] corr_pattern = patternZoomCropPad(
corr_patterns[n] = patternZoomCropPad(
patterns[n], // double [] pattern,
pattern_size, // int pattern_size,
corr_size, // int size,
zoomout, // int zoomout,
true); // out_normalize); // boolean normalize)
false); // true); // out_normalize); // boolean normalize)
if (kernel != null) {
corr_pattern = convolveWithKernel(
corr_pattern, // final double [] data,
kernel, // final double [] kernel,
corr_size); // final int width)
corr_patterns[n] = convolveWithKernel(
corr_patterns[n], // final double [] data,
kernel, // final double [] kernel,
corr_size); // final int width)
}
corrs_out[n]= correlateWithPattern(
data, // final double [] data,
width, // final int width,
corr_size, // final int psize, // power of 2, such as 64
corr_pattern, //final double [] pattern, // [psize*psize]
phaseCoeff, // final double phaseCoeff,
debugLevel); // final int debugLevel) {
data, // final double [] data,
width, // final int width,
corr_size, // final int psize, // power of 2, such as 64
corr_patterns[n], //final double [] pattern, // [psize*psize]
false, // final boolean convolve, // convolve, not correlate
phaseCoeff, // final double phaseCoeff,
debugLevel); // final int debugLevel) {
}
ShowDoubleFloatArrays.showArrays(
corrs_out,
width,
height,
true,
removeKnownExtension(imp_src.getTitle())+"-PATTERN_CORR-0");
//data
double [] bestcorr = new double [data.length];
int [] bestpatt = new int [bestcorr.length];
ArrayList<Point> object_list = combineDirCorrs (
corrs_out, // final double [][] corrs,
width, // final int width,
bestcorr, // final double [] bestcorr
bestpatt, // final int [] bestpatt_in,
min_corr, // final double min_corr,
adv_radius, // final double adv_radius,
corr_radius); // final int corr_radius)
System.out.println("testPatternCorrelate(): Found "+object_list.size()+" candidates");
ShowDoubleFloatArrays.showArrays(
bestcorr,
width,
height,
removeKnownExtension(imp_src.getTitle())+"-BEST_CORR");
double [][] filtered_corrs = splitPattenCorrelations(
patterns.length, // final int num_patterns,
bestcorr, // final double [] bestcorr,
bestpatt); // final int [] bestpatt)
/*
for (int n = 0; n < patterns.length; n++) {
lim_corr[n] = corrs_out[n].clone();
for (int i = 0; i < lim_corr.length; i++) {
if (lim_corr[n][i] < min_corr) {
lim_corr[n][i] = 0;
}
}
}
String [] rslt_titles =new String [patterns.length+1];
double [][] rslt = new double [patterns.length + 1][]; // ]{data, corr_out};
*/
for (int n = 0; n < patterns.length; n++) {
convolve_out[n]= correlateWithPattern(
filtered_corrs[n],// corrs_out[n], // final double [] data,
width, // final int width,
corr_size, // final int psize, // power of 2, such as 64
corr_patterns[n], // final double [] pattern, // [psize*psize]
true, // final boolean convolve, // convolve, not correlate
phaseCoeff, // final double phaseCoeff,
debugLevel); // final int debugLevel) {
}
double [] convolved_merged = mergePattenConvolutions(
convolve_out); // final double [][] convs)
String [] rslt_titles =new String [patterns.length+2];
double [][] rslt = new double [patterns.length + 2][]; // ]{data, corr_out};
double [][] conv_rslt = new double [patterns.length + 2][]; // ]{data, corr_out};
rslt[0] = data;
conv_rslt[0] = data;
rslt_titles[0] = "original";
rslt[1] = bestcorr;
conv_rslt[1] = convolved_merged;
rslt_titles[1] = "combined";
for (int i = 0; i < patterns.length; i++) {
rslt[i+1] = corrs_out[i];
rslt_titles[i+1] = "corr_"+pattern_labels[i];
rslt[i+2] = filtered_corrs[i];
conv_rslt[i+2] = convolve_out[i];
rslt_titles[i+2] = "corr_"+pattern_labels[i];
}
ShowDoubleFloatArrays.showArrays(
rslt,
......@@ -3206,14 +3302,226 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
true,
removeKnownExtension(imp_src.getTitle())+"-PATTERN_CORR",
rslt_titles);
ShowDoubleFloatArrays.showArrays(
conv_rslt,
width,
height,
true,
removeKnownExtension(imp_src.getTitle())+"-PATTERN_CORR_CONVOLVE_min_corr"+min_corr,
rslt_titles);
System.out.println("testPatternCorrelate(): correlation done");
}
/**
* Combine results of the phase correlations with different patters
* (now with semi-obscured half-circles), using the best (highest value match),
* filter by correlation value and suppress weeaker maximums by stronger neighbors
* @param corrs array of 2D correlation results [num_pattern][pixel]
* @param width width of the correlation array
* @param bestcorr_in null or double array [corrs[0].length] will contain
* best correlation value from multiple patterns of corrs[][].
* @param bestbatt_in null or integer array [corrs[0].length] will contain
* number of the best pattern
* @param min_corr - minimal correlation value to use (lower - replace with 0)
* @param adv_radius - suppress weaker local maximums near the stronger one
* @param corr_radius - maximal correlation peak radius (zero out outside).
* Should be <= adv_radius to avoid thread races
* @return ArrayList of integer pairs - pixel index (x+y*width) and the pattern index
*/
public static ArrayList<Point> combineDirCorrs (
final double [][] corrs,
final int width,
final double [] bestcorr_in,
final int [] bestpatt_in,
final double min_corr,
final double adv_radius,
final int corr_radius) {
// final int num_patt = corrs.length;
final int height = corrs[0].length/width;
final double [] bestcorr = (bestcorr_in != null) ? bestcorr_in: (new double [corrs[0].length]);
final int [] bestpatt = (bestpatt_in != null) ? bestpatt_in: (new int [corrs[0].length]);
// before adversarial filtering
final double [] bestcorr0 = new double [corrs[0].length];
final int [] bestpatt0 = new int [corrs[0].length];
Arrays.fill(bestpatt, -1);
Arrays.fill(bestpatt0, -1); // maybe not needed
// final boolean [] local_max = new boolean [bestcorr.length];
final Thread[] threads = ImageDtt.newThreadArray();
final AtomicInteger ai = new AtomicInteger(0);
final TileNeibs tn = new TileNeibs(width, height);
// Combine patterns correlations
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
for (int ipix = ai.getAndIncrement(); ipix < bestcorr0.length; ipix = ai.getAndIncrement()) {
for (int n = 0; n < corrs.length; n++) {
if (corrs[n][ipix] > bestcorr0[ipix]) {
bestcorr0[ipix] = corrs[n][ipix];
bestpatt0[ipix] = n;
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
// mark all local max
ai.set(0);
final AtomicInteger ati = new AtomicInteger(0);
ArrayList<ArrayList<Point>> lists = new ArrayList<ArrayList<Point>>();
for (int i = 0; i < threads.length; i++) {
lists.add(new ArrayList<Point>());
}
final int iradius = (int) Math.ceil(adv_radius);
final double adv_radius2 = adv_radius * adv_radius;
final double corr_radius2 = corr_radius*corr_radius;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
public void run() {
int thread_num = ati.getAndIncrement();
final TileNeibs tn = new TileNeibs(width, height);
int corr_diameter = 2*corr_radius+1;
TileNeibs tnl = new TileNeibs(corr_diameter,corr_diameter);
int local_center = (2 * corr_radius + 2) * corr_radius;
boolean [] mask = new boolean [corr_diameter*corr_diameter];
for (int ipix = ai.getAndIncrement(); ipix < bestcorr0.length; ipix = ai.getAndIncrement()) if (bestcorr0[ipix] >= min_corr) {
if (ipix == -792148) {
System.out.println("combineDirCorrs(): ipix="+ipix);
}
is_max: {
for (int dir = 0; dir > TileNeibs.DIRS; dir++) {
int ipix1 = tn.getNeibIndex(ipix, dir);
if ((ipix1 >= 0) && (bestcorr0[ipix1] > bestcorr0[ipix])){
break is_max;
}
}
// compare with others in adversarial radius, if equal - use higher index ipix
for (int dy = -iradius; dy <= iradius; dy++) {
for (int dx = -iradius; dx <= iradius; dx++) {
double r2 = dx*dx + dy*dy;
if ((r2 < adv_radius2) && (r2 > 2.1)) {
int ipix1 = tn.getNeibIndex(ipix, dx, dy);
if (ipix1 >= 0) {
if (bestcorr0[ipix1] > bestcorr0[ipix]) {
break is_max;
} else if ((bestcorr0[ipix1] == bestcorr0[ipix]) && (ipix1 < ipix)) {
break is_max;
}
}
}
}
}
int patt = bestpatt0[ipix];
lists.get(thread_num).add(new Point(ipix,patt));
// Below should be no conflicts as no intersect
Arrays.fill(mask, false);
for (int dy = -corr_radius; dy <= corr_radius; dy++) {
for (int dx = -corr_radius; dx <= corr_radius; dx++) {
double r2 = dx*dx + dy*dy;
if (r2 < corr_radius2) {
int ipix1 = tn.getNeibIndex(ipix, dx, dy);
// TODO: mark mask with continuous pixels
if (ipix1 >= 0) {
int li = tnl.getNeibIndex(local_center, dx, dy); // should be always
if (li < 0) {
System.out.println("combineDirCorrs(): BUG1 li<0");
} else {
mask[li] = corrs[patt][ipix1] >= min_corr; // only consider the same pattern
}
}
}
}
}
int [] clusters = tnl.enumerateClusters(
mask, // boolean [] tiles,
null, // int [] num_clusters,
false); // boolean ordered)
int center_cluster = clusters[local_center];
// Copy continuous pixels from the selected pattern
for (int dy = -corr_radius; dy <= corr_radius; dy++) {
for (int dx = -corr_radius; dx <= corr_radius; dx++) {
int li = tnl.getNeibIndex(local_center, dx, dy); // should be always
if (li < 0) {
System.out.println("combineDirCorrs(): BUG2 li<0");
} else {
if (clusters[li] == center_cluster) {
int ipix1 = tn.getNeibIndex(ipix, dx, dy);
if (ipix1 >= 0) {
bestcorr[ipix1] = corrs[patt][ipix1];
bestpatt[ipix1] = patt;
}
}
}
}
}
}
}
}
};
}
ImageDtt.startAndJoin(threads);
// combine all lists
ArrayList<Point> list = new ArrayList<Point>();
for (int i = 0; i < lists.size(); i++) {
list.addAll(lists.get(i));
}
return list;
}
/**
* Split correlation filter results into per-pattern ones before convolution.
* Temporary solution as convolution is needed only for a few found maximums
* @param num_patterns number of patters (should match bestpatt maximal value +1)
* @param bestcorr best correlation values
* @param bestpatt best pattern index matching bestcorr
* @return [num_patterns][] arrays of per-pattern correlation results to be convolved
* with individual patterns
*/
public static double [][] splitPattenCorrelations(
final int num_patterns,
final double [] bestcorr,
final int [] bestpatt){
final double [][] corrs = new double [num_patterns][bestcorr.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 < bestcorr.length; ipix = ai.getAndIncrement()) if (bestcorr[ipix] > 0){
corrs[bestpatt[ipix]][ipix] = bestcorr[ipix];
}
}
};
}
ImageDtt.startAndJoin(threads);
return corrs;
}
public static double [] mergePattenConvolutions(
final double [][] convs){
final double [] combo_conv = new double [convs[0].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 < combo_conv.length; ipix = ai.getAndIncrement()){
for (int n = 0; n < convs.length; n++) {
combo_conv[ipix] += convs[n][ipix];
}
}
}
};
}
ImageDtt.startAndJoin(threads);
return combo_conv;
}
public static void testPatternGenerate() {
int half_size = 100;
boolean half_pix = false; // center between pixels
double radius = 32; // 32;
double radius = 30; // 32; // 32;
double edge= 15; // 4;
double radius_in= 15; // if 0 - skip
double edge_in = 8;
......@@ -3318,7 +3626,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
String [] pattern_titles = new String [patterns.length];
pattern_titles[0] = "full";
for (int n = 1; n < pattern_titles.length; n++) {
pattern_titles[n] = String.format("%f5.1 deg", 360.0*(n-1)/halves_number);
pattern_titles[n] = String.format("%5.1f deg", 360.0*(n-1)/halves_number);
}
// public static ImagePlus makeArrays(double[][] pixels, int width, int height, String title, String [] titles) {
......@@ -3382,7 +3690,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
}
if (normalize) {
// double scale = 1.0/sum_pix;
double scale = -10000.0/sum_pix; //TODO: calculate better
double scale = -500000.0/sum_pix; //TODO: calculate better
for (int i = 0; i < pout.length; i++) {
pout[i] *= scale;
}
......
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