Skip to content

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

Cleaned up pattern matching, implemented re-centering for correlations, 

working snapshot
parent c8f99e6d
Show whitespace changes
Inline Side-by-side
Showing with 1871 additions and 1419 deletions
src/main/java/com/elphel/imagej/common/DoubleFHT.java
View file @ 640b8df3
......@@ -293,6 +293,22 @@ public class DoubleFHT {
return result;
}
public double [] getFrequencyFilter (
double [] data, // just to find out size
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 filter;
}
/**
* Invert kernel for deconvolutions. Kernel is assumed to fade near edges
* @param data source kernel, square, power of 2 sides. Data is destroyed
......@@ -314,7 +330,6 @@ public class DoubleFHT {
filter =this.freqMask;
}
return invert (data, fat_zero, filter);
}
/**
......@@ -401,8 +416,6 @@ public class DoubleFHT {
}
public double[] phaseCorrelate(double[] first, double phaseCoeff, double high_pass, double low_pass, // high/low
// pass
// filtering
double[] fht_save) { // null-OK
updateMaxN(first);
double[] filter = null;
......@@ -503,7 +516,10 @@ public class DoubleFHT {
* @return phase correlation, first still contains original transformed!
*/
public double[] phaseCorrelatePattern( // new
double[] first, double[] secondFD, double phaseCoeff, double[] filter, // high/low pass filtering
double[] first,
double[] secondFD,
double phaseCoeff,
double[] filter, // high/low pass filtering
double[] first_save) { // null-OK
if (first.length != secondFD.length) {
IJ.showMessage("Error", "Correlation arrays should be the same size");
......@@ -537,7 +553,9 @@ public class DoubleFHT {
*/
public double[] convolvePattern( // new
double[] first, double[] secondFD, double[] filter, // high/low pass filtering
double[] first,
double[] secondFD,
double[] filter, // high/low pass filtering
double[] first_save) { // null-OK
if (first.length != secondFD.length) {
IJ.showMessage("Error", "Correlation arrays should be the same size");
......@@ -595,11 +613,16 @@ public class DoubleFHT {
return first;
}
public double[] convolve(double[] first, double[] second) {
public double[] convolve(
double[] first,
double[] second) {
return convolve(first, second, null);
}
public double[] convolve(double[] first, double[] second, double[] filter) { // high/low pass filtering
public double[] convolve(
double[] first,
double[] second,
double[] filter) { // high/low pass filtering
// System.out.println("correlate");
if (first.length != second.length) {
IJ.showMessage("Error", "Correlation arrays should be the same size");
......
src/main/java/com/elphel/imagej/orthomosaic/ComboMatch.java
View file @ 640b8df3
......@@ -544,13 +544,13 @@ public class ComboMatch {
}
imp_obj.setRoi(roi);
// calculate statistics for cor_ret[i]
double [][] stats = new double [object_stack.length][];
CorrelationPeakStats [] stats = new CorrelationPeakStats [object_stack.length];
double search_rad = 15.0;
double frac_max = 0.5;
double other_rad = 25;
for (int i = 0; i < object_stack.length; i++) {
double [] a_cent= {corr_size/2+centers[i][0],corr_size/2+centers[i][1]};
stats[i]=ObjectLocation.getMaxLocArea(
stats[i]=new CorrelationPeakStats(
corr_ret[i], // double [] data, // square data
a_cent, // centers[i], // double [] cent_xy, // if null, use center of the square
search_rad, // double radius, // search for maximum within this radius
......@@ -567,7 +567,8 @@ public class ComboMatch {
ObjectLocation ol = object_list.get(i);
System.out.println(String.format("%2d (%4d/%4d): %7.5f %7.3f %7.3f %7.3f %7.3f %7.3f %17s",
i, ol.getPixels()[0], ol.getPixels()[1],
stats[i][0],stats[i][1],stats[i][2],stats[i][3],stats[i][4]+centers[i][0],stats[i][5]+centers[i][1],
stats[i].best_d,stats[i].eff_rad,stats[i].elong,stats[i].dist,
stats[i].cent_offs[0]+centers[i][0],stats[i].cent_offs[1]+centers[i][1],
ol.getName()));
}
System.out.println();
......@@ -868,7 +869,8 @@ public class ComboMatch {
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]
null); // warp); // FineXYCorr warp)
null, // warp); // FineXYCorr warp)
null); // double [][] ground_planes) TODO: add calculation of the ground plane for single images
// imp_pat_match.show();
}
......@@ -895,11 +897,16 @@ public class ComboMatch {
boolean batch_mode = true; // false; // true;
boolean ignore_prev_rms = true;
Rectangle woi = new Rectangle(); // used to return actual woi from correlateOrthoPair()
double [][] ground_planes = null;
for (int zi = 0; zi < zooms.length; zi++) {
zoom_lev = zooms[zi];
if (zoom_lev >=1000) {
break;
}
boolean show_vf = render_match || pattern_match;
if (render_match || pattern_match) {
ground_planes = new double [gpu_pair.length][];
}
// will modify affines[1], later add jtj, weight, smth. else?
FineXYCorr warp = maps_collection.correlateOrthoPair(
clt_parameters, // CLTParameters clt_parameters,
......@@ -914,7 +921,8 @@ public class ComboMatch {
affines, // double [][][] affines, // on top of GPS offsets
woi, // Rectangle woi,
zoom_lev, // int zoom_lev,
show_vf, // boolean show_vf,
ground_planes, // double [][] ground_planes, // null or double[2] - will return ground planes
debugLevel); // final int debugLevel)
if (warp == null) {
System.out.println("Failed correlateOrthoPair()");
......@@ -951,7 +959,8 @@ public class ComboMatch {
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)
warp, // FineXYCorr warp)
ground_planes); // double [][] ground_planes); // null or
// imp_pat_match.show();
}
if (render_match) {
......
src/main/java/com/elphel/imagej/orthomosaic/CorrelationPeakStats.java 0 → 100644
View file @ 640b8df3
package com.elphel.imagej.orthomosaic;
import java.util.Arrays;
import com.elphel.imagej.tileprocessor.TileNeibs;
import Jama.EigenvalueDecomposition;
import Jama.Matrix;
public class CorrelationPeakStats {
public double best_d;
public double eff_rad;
public double elong;
public double dist;
public double [] cent_offs;
public double max_other;
public CorrelationPeakStats(
double [] data, // square data
double [] cent_xy, // if null, use center of the square
double radius, // search for maximum within this radius
double frac_max,
double other_radius,
int debugLevel) {
double [] stats = getStatsNaN(
data, // square data
cent_xy, // if null, use center of the square
radius, // search for maximum within this radius
frac_max,
other_radius,
debugLevel);
best_d = stats[0];
eff_rad = stats[1];
elong = stats[2];
dist = stats[3];
cent_offs = new double [] {stats[4], stats[5]};
max_other = stats[6];
}
int [] getIntOffset() {
return new int[] {(int) Math.round(cent_offs[0]), (int) Math.round(cent_offs[1])};
}
/**
* Search for maximum within specified radius from the center or specified offset from the center
* @param data square data
* @param cent_xy offset from the square center to search
* @param radius search for maximum within this radius
* @param frac_max fraction of maximum to measure area
* @param other_radius measure maximal disconnected pixel fraction of the maximum within this radius
* @return {maximum, area, dist, x, y}. No interpolation yet, each returned value is integer.
* Returns null if no local maximum within area. x,y are offsets from the (provided) center
*/
public static double [] getStatsNaN(
double [] data, // square data
double [] cent_xy, // if null, use center of the square
double radius, // search for maximum within this radius
double frac_max,
double other_radius,
int debugLevel) {
double [] rslt = getStats(
data, // square data
cent_xy, // if null, use center of the square
radius, // search for maximum within this radius
frac_max,
other_radius,
debugLevel);
if (rslt != null) {
return rslt;
} else {
return new double[] {Double.NaN,Double.NaN,Double.NaN,Double.NaN,Double.NaN,Double.NaN,Double.NaN};
}
}
public static double [] getStats(
double [] data, // square data
double [] cent_xy, // if null, use center of the square
double radius, // search for maximum within this radius
double frac_max,
double other_radius, // may be 0 if not needed
int debugLevel) {
boolean debug = debugLevel > 1;
int size = (int) Math.sqrt(data.length);
if (cent_xy == null) {
cent_xy = new double [] {size/2, size/2};
}
double maxr2 = radius * radius;
double best_d = 0;
int best_indx = -1;
int min_x = Math.max(1, (int) Math.floor(cent_xy[0]-radius));
int min_y = Math.max(1, (int) Math.floor(cent_xy[1]-radius));
int max_x = Math.min(size-2, (int) Math.ceil(cent_xy[0]+radius));
int max_y = Math.min(size-2, (int) Math.ceil(cent_xy[1]+radius));
for (int y = min_y; y <= max_y; y++) { // do not search on very edges
double dy = (y-cent_xy[1]);
double y2 = dy*dy;
if (y2 < maxr2) {
for (int x = min_x; x <= max_x; x++) {
double dx = x - cent_xy[0];
double r2 = y2 + dx*dx;
if (r2 < maxr2) {
int indx = y * size + x;
double d = data[indx];
if (d > best_d) {
best_indx = indx;
best_d = d;
}
}
}
}
}
if (best_indx < 0) {
return null;
}
// is it local max?
if ((data[best_indx - 1] > best_d) || (data[best_indx + 1] > best_d) ||
(data[best_indx - size] > best_d) || (data[best_indx + size] > best_d)) {
return null; // on the edge, not a local max
}
boolean [] above_thresh = new boolean [data.length];
double thresh = best_d * frac_max;
for (int i = 0; i < data.length; i++) {
above_thresh[i] = data[i] > thresh;
}
int [] clusters = (new TileNeibs(size,size)).enumerateClusters(
above_thresh, // boolean [] tiles,
null, // int [] num_clusters,
false); // boolean ordered)
int center_cluster = clusters[best_indx];
int best_x = best_indx % size;
int best_y = best_indx / size;
double xc = best_x - cent_xy[0];
double yc = best_y - cent_xy[1];
double s0=0, sx=0, sy=0, sx2 = 0, sy2=0, sxy = 0;
for (int i = 0; i < clusters.length; i++) {
if (clusters[i] == center_cluster) {
double y = i / size - (yc + cent_xy[1]); //(yc + cent_xy[1]) - absolute, from (0,0)
double x = i % size - (xc + cent_xy[0]);
double w = data[i]-thresh;
s0 += w;
sx += w * x;
sy += w * y;
sx2 += w * x * x;
sy2 += w * y * y;
sxy += w * x * y;
}
}
double cxx = sx2 - sx * sx / s0, cyy= sy2 - sy * sy / s0, cxy = sxy - sx * sy / s0;
/*
* sum(Mi*(Xi-avg(X))^2) = SX2 - SX^2/S0
* sum(Mi*(Yi-avg(Y))^2) = SY2 - SY^2/S0
* sum(Mi*(Xi-avg(X)*(Yi-avg(Y))) = SXY - SX*SY / S0
*/
Matrix covar = new Matrix(new double[][] {{cxx, cxy},{cxy,cyy}});
double [] cent_offs = {sx/s0 + xc,sy/s0 + yc};
EigenvalueDecomposition eig = covar.eig();
double [] eigval = {eig.getD().get(0, 0),eig.getD().get(1, 1)};
Arrays.sort(eigval); // ascending
double elong = Math.sqrt(eigval[1]/eigval[0]);
double eff_rad = Math.sqrt(Math.sqrt(eigval[1]*eigval[0]));
double dist = Math.sqrt(cent_offs[0]*cent_offs[0] + cent_offs[1]*cent_offs[1]);
double max_other = 0.0;
if (other_radius > 0) { // only calculate if asked for
double or2 = other_radius*other_radius;
min_x = Math.max(1, (int) Math.floor(best_x-other_radius));
min_y = Math.max(1, (int) Math.floor(best_y-other_radius));
max_x = Math.min(size-2, (int) Math.ceil(best_x+other_radius));
max_y = Math.min(size-2, (int) Math.ceil(best_y+other_radius));
for (int y = min_y; y <= max_y; y++) { // do not search on very edges
double dy = (y-best_y);
double y2 = dy*dy;
if (y2 < or2) {
for (int x = min_x; x <= max_x; x++) {
double dx = x - best_x;
double r2 = y2 + dx*dx;
if (r2 < or2) {
int indx = y * size + x;
if (clusters[indx] != center_cluster) { // only disconnected from the central area
double d = data[indx];
if (d > max_other) {
max_other = d;
}
}
}
}
}
}
max_other /= best_d;
}
if (debug) {
System.out.println("\ncenter offset ["+(sx/s0)+","+(sy/s0)+"] , from center: ["+cent_offs[0]+","+cent_offs[1]+"]");
System.out.println("Covariance matrix:");
covar.print(8, 6);
System.out.println("eig.getV()");
eig.getV().print(8, 6);
System.out.println("eig.getD()");
eig.getD().print(8, 6);
System.out.println(String.format("best_d=%7.5f, rad=%6.3f, elong=%6.3f, dist=%6.3f, dx=%6.3f, dy=%6.3f, mo=%5.3f",
best_d, eff_rad, elong, dist, cent_offs[0], cent_offs[1], max_other));
}
return new double [] {best_d, eff_rad, elong, dist, cent_offs[0], cent_offs[1], max_other} ;
}
}
src/main/java/com/elphel/imagej/orthomosaic/ItemMatch.java
View file @ 640b8df3
......@@ -4,12 +4,111 @@ import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import ij.ImagePlus;
import ij.Prefs;
import ij.gui.PointRoi;
import ij.io.FileSaver;
public class ItemMatch {
public double [] combo_pxy; // combo image {pixel_x, pixel_y}
public double [] lla; // add later
public HashMap<String, ItemPatternMatch> pattern_matches = new HashMap<String, ItemPatternMatch>();
public double abs_contrast = 0;
public ItemMatch (double [] combo_pxy) {
public CorrelationPeakStats[][] filter_data; // [scene][half/full]
public double [][] extracted_objects; // [scene][pix]
public double [][] extracted_nodc; // [scene][pix]
public double [][] extracted_masked; // [scene][pix]
public double [][][] extr_corr; // [scene][pix][patt]
public double [][] extr_corr_half; // [scene][pix][h/f]
public boolean removed = false;
public ItemMatch (
int num_scenes,
double [] combo_pxy) {
setComboPXY(combo_pxy);
this.combo_pxy = combo_pxy;
filter_data = new CorrelationPeakStats[num_scenes][2];
extracted_objects = new double [num_scenes][];
extracted_nodc = new double [num_scenes][];
extracted_masked = new double [num_scenes][];
extr_corr = new double [num_scenes][][];
extr_corr_half = new double [num_scenes][];
removed = false;
}
public void setRemoved (boolean removed) {
this.removed = removed;
}
public void remove () {
this.removed = true;
}
public boolean isRemoved() {
return removed;
}
public void setCorrFull(
int scene_num,
int num_patt,
double [] corr_data) {
extr_corr[scene_num] = new double [num_patt][];
extr_corr[scene_num][0] = corr_data;
}
public void setCorr(
int scene_num,
int patt_index,
double [] corr_data) {
extr_corr[scene_num][patt_index] = corr_data;
}
public double [] getCorr(
int scene_num,
int patt_index) {
return extr_corr[scene_num][patt_index];
}
public double [][] getCorrs(
int scene_num) {
return extr_corr[scene_num];
}
public void setCorrHalf(
int scene_num,
int best_patt, // 0-full
boolean combine_full) {
extr_corr_half[scene_num] = extr_corr[scene_num][best_patt].clone();
if (combine_full) {
double full_max_value = filter_data[scene_num][0].best_d; // [0];
if (!Double.isNaN(full_max_value)) {
// multiply by full correlation, delete by full correlation maximum
for (int i = 0; i < extr_corr_half[scene_num].length; i++){
extr_corr_half[scene_num][i] *=
Math.max(0, extr_corr[scene_num][0][i]/full_max_value);
}
}
}
}
public void setCorrHalf(
int scene_num,
double [] corr_data) {
extr_corr_half[scene_num] = corr_data;
}
public double [] getCorrHalf(
int scene_num) {
return extr_corr_half[scene_num];
}
public void setComboPXY(double [] combo_pxy) {
this.combo_pxy = combo_pxy;
}
......@@ -27,7 +126,7 @@ public class ItemMatch {
return new int [] {(int) Math.round(combo_pxy[0]), (int) Math.round(combo_pxy[1])};
}
public void addPatternMatches(
public void addPatternMatches( // updates existing too
GroundObjectPattern groundObjectPattern,
double [] matches,
int best_sub) {
......@@ -101,6 +200,8 @@ public class ItemMatch {
return match.getMatch(indx);
}
/**
* Return correlation value as the one of the pattern correlations pointed
* by the best subpattern index.
......@@ -126,28 +227,206 @@ public class ItemMatch {
public static ArrayList<Integer> sortByMatch(
ArrayList<ItemMatch> match_list,
GroundObjectPattern groundObjectPattern,
boolean keep_removed,
int indx){
String pattern_path = groundObjectPattern.getPatternPath();
return sortByMatch(match_list, pattern_path, indx);
return sortByMatch(match_list, pattern_path, keep_removed, indx);
}
public static ArrayList<Integer> sortByMatch(
ArrayList<ItemMatch> match_list,
String pattern_path,
boolean keep_removed,
int indx){
double large_enough = 1.0;
ArrayList<Integer> result_list= new ArrayList<Integer>(match_list.size());
for (int i = 0; i < match_list.size(); i++) {
if (keep_removed || !match_list.get(i).isRemoved()) {
result_list.add(i);
}
}
Collections.sort(result_list, new Comparator<Integer>() {
@Override
public int compare(Integer lhs, Integer rhs) {
double rhsd = match_list.get(rhs).getMatchValue(pattern_path, indx);
double lhsd = match_list.get(lhs).getMatchValue(pattern_path, indx);
return (rhsd > lhsd) ? 1 : (rhsd < lhsd) ? -1 : 0; // recreasing
ItemMatch rhsm = match_list.get(rhs);
ItemMatch lhsm = match_list.get(lhs);
double rhsd = rhsm.getMatchValue(pattern_path, indx) - (rhsm.isRemoved() ? large_enough:0 );
double lhsd = lhsm.getMatchValue(pattern_path, indx) - (lhsm.isRemoved() ? large_enough:0 );
return (rhsd > lhsd) ? 1 : (rhsd < lhsd) ? -1 : 0; // decreasing
}
});
return result_list;
}
public static void setExtractedObjects(
int scene_num,
int extr_size,
GroundObjectPattern gop,
ArrayList <ItemMatch> matches_list,
ArrayList<Integer> match_sort,
double[] data,
int width,
int debugLevel) {
if (match_sort.isEmpty()) {
System.out.println("List of detected objects is empty, nothing to process");
return;
}
for (int mn=0; mn < match_sort.size(); mn++) {
int indx =match_sort.get(mn);
ItemMatch match = matches_list.get(indx);
match.extractObject(
data, // double [] data,
width, // int width,
scene_num, // int scene_num,
extr_size); // int extr_size);
}
}
public void extractObject(
double [] data,
int width,
int scene_num,
int extr_size) {
double [] center_xy = getXY();
double [] extr_data = ObjectLocation.extractObjectImage(
center_xy, // double [] center_xy,
extr_size, // int size,
data, // double [] src_img,
width); // int width)
setExtractedObject(extr_data, scene_num);
}
public void extractCorrs(
double [][] corr_data,
int width,
int scene_num,
int extr_size) {
double [] center_xy = getXY();
extr_corr[scene_num] = new double [corr_data.length][];
for (int i = 0; i < corr_data.length; i++) {
extr_corr[scene_num][i] = ObjectLocation.extractObjectImage(
center_xy, // double [] center_xy,
extr_size, // int size,
corr_data[i], // double [] src_img,
width); // int width)
}
}
public void setExtractedObject(
double [] data,
int scene_num) {
extracted_objects[scene_num] = data;
extracted_nodc[scene_num] = data.clone();
OrthoMap.removeDC(extracted_nodc[scene_num]);
}
public void setMaskedObject(
double [] data,
int scene_num) {
extracted_masked[scene_num] = data;
}
public void setCorrData(
double [] data,
int scene_num,
int half_full) { // 1 - half, 0 - full
extr_corr[scene_num][half_full] = data;
}
public static ImagePlus getImageExtracts(
String prefix, // include PC, filter here
String suffix, // w/o .tiff
boolean show,
boolean save,
boolean show_removed,
String save_dir,
int scene_num,
boolean nodc,
int mode, // 0 - extract, 1 - masked, 2 - corr full, 3 - corr half
boolean show_centers,
int extr_size,
GroundObjectPattern gop,
ArrayList <ItemMatch> matches_list,
ArrayList<Integer> match_sort,
int debugLevel) {
if (match_sort.isEmpty()) {
System.out.println("List of detected objects is empty, nothing to display");
return null;
}
PointRoi roi = new PointRoi();
roi.setOptions("nolabel"); // label");
for (int mn=0; mn < match_sort.size(); mn++) {
roi.addPoint(extr_size/2,extr_size/2,mn+1); // ,1);
}
int num_slices = match_sort.size();
if (!show_removed) {
num_slices=0;
for (int mn=0; mn < match_sort.size(); mn++) {
if (!matches_list.get(match_sort.get(mn)).isRemoved()) {
num_slices++;
}
}
}
String [] extr_titles = new String [num_slices];
double [][] extr_data = new double [num_slices][];
String title = prefix;
switch (mode) {
case 0: title += "_objects_" + scene_num + (nodc ? "_no-DC":"_with-DC"); break;
case 1: title += "_contrast-masks"; break;
case 2: title += "_objects_"+scene_num+"_corr_full"; break;
case 3: title += "_objects_"+scene_num+"_corr_halves"; break;
}
if (suffix != null) {
title += suffix;
}
title += ".tiff";
int indx = 0;
for (int mn=0; mn < match_sort.size(); mn++) {
ItemMatch match = matches_list.get(match_sort.get(mn));
if (show_removed || !matches_list.get(match_sort.get(mn)).isRemoved()) {
int [] ixy = match.getIntXY();
extr_titles[indx] = ixy[0]+"/"+ ixy[1]+":"+String.format("%7.5f", match.getMatchValue(gop));
switch (mode) {
case 0: extr_data[indx] = nodc ? match.extracted_nodc[scene_num] : match.extracted_objects[scene_num]; break;
case 1: extr_data[indx] = match.extracted_masked[scene_num]; break;
case 2: extr_data[indx] = match.extr_corr[scene_num][0]; break;
case 3: extr_data[indx] = match.extr_corr_half[scene_num]; break;
}
indx++;
}
}
ImagePlus imp = ShowDoubleFloatArrays.makeArrays(
extr_data,
extr_size,
extr_size,
title, // +"-objects_"+scene_num,
extr_titles); // test_titles,
if (imp == null) {
return null;
}
if (show_centers) {
imp.setRoi(roi); // null
}
if (show) {
imp.show();
}
if (save && (save_dir != null)) {
if (!save_dir.endsWith(Prefs.getFileSeparator())) {
save_dir+=Prefs.getFileSeparator();
}
String save_path = save_dir+title;
FileSaver imp_fs = new FileSaver(imp);
imp_fs.saveAsTiff(save_path);
if (debugLevel > -4) {
System.out.println("Saved "+save_path);
}
}
return imp;
}
}
src/main/java/com/elphel/imagej/orthomosaic/ItemPatternMatch.java
View file @ 640b8df3
......@@ -13,7 +13,8 @@ public class ItemPatternMatch {
public void setMatches(double [] matches) {
sub_matches = matches;
}
public void setMatch(int indx, double match_value) {
public void setMatch(int indx, double match_value) { // not used
if (sub_matches == null) {
sub_matches = new double [indx+1];
Arrays.fill(sub_matches, Double.NaN);
......
src/main/java/com/elphel/imagej/orthomosaic/ObjectLocation.java
View file @ 640b8df3
......@@ -119,239 +119,6 @@ public class ObjectLocation {
}
return dcrop;
}
/**
* Search for maximum within specified radius from the center or specified offset from the center
* @param data square data
* @param cent_xy offset from the square center to search
* @param radius search for maximum within this radius
* @param frac_max fraction of maximum to measure area
* @param other_radius measure maximal disconnected pixel fraction of the maximum within this radius
* @return {maximum, area, dist, x, y}. No interpolation yet, each returned value is integer.
* Returns null if no local maximum within area. x,y are offsets from the (provided) center
*/
public static double [] getMaxLocAreaNaN(
double [] data, // square data
double [] cent_xy, // if null, use center of the square
double radius, // search for maximum within this radius
double frac_max,
double other_radius,
int debugLevel) {
double [] rslt = getMaxLocArea(
data, // square data
cent_xy, // if null, use center of the square
radius, // search for maximum within this radius
frac_max,
other_radius,
debugLevel);
if (rslt != null) {
return rslt;
} else {
return new double[] {Double.NaN,Double.NaN,Double.NaN,Double.NaN,Double.NaN,Double.NaN,Double.NaN};
}
}
public static double [] getMaxLocArea(
double [] data, // square data
double [] cent_xy, // if null, use center of the square
double radius, // search for maximum within this radius
double frac_max,
double other_radius,
int debugLevel) {
boolean debug = debugLevel > 1;
int size = (int) Math.sqrt(data.length);
if (cent_xy == null) {
cent_xy = new double [] {size/2, size/2};
}
double maxr2 = radius * radius;
double best_d = 0;
int best_indx = -1;
int min_x = Math.max(1, (int) Math.floor(cent_xy[0]-radius));
int min_y = Math.max(1, (int) Math.floor(cent_xy[1]-radius));
int max_x = Math.min(size-2, (int) Math.ceil(cent_xy[0]+radius));
int max_y = Math.min(size-2, (int) Math.ceil(cent_xy[1]+radius));
for (int y = min_y; y <= max_y; y++) { // do not search on very edges
double dy = (y-cent_xy[1]);
double y2 = dy*dy;
if (y2 < maxr2) {
for (int x = min_x; x <= max_x; x++) {
double dx = x - cent_xy[0];
double r2 = y2 + dx*dx;
if (r2 < maxr2) {
int indx = y * size + x;
double d = data[indx];
if (d > best_d) {
best_indx = indx;
best_d = d;
}
}
}
}
}
// is it local max?
if ((data[best_indx - 1] > best_d) || (data[best_indx + 1] > best_d) ||
(data[best_indx - size] > best_d) || (data[best_indx + size] > best_d)) {
return null; // on the edge, not a local max
}
boolean [] above_thresh = new boolean [data.length];
double thresh = best_d * frac_max;
for (int i = 0; i < data.length; i++) {
above_thresh[i] = data[i] > thresh;
}
int [] clusters = (new TileNeibs(size,size)).enumerateClusters(
above_thresh, // boolean [] tiles,
null, // int [] num_clusters,
false); // boolean ordered)
int center_cluster = clusters[best_indx];
int best_x = best_indx % size;
int best_y = best_indx / size;
double xc = best_x - cent_xy[0];
double yc = best_y - cent_xy[1];
double s0=0, sx=0, sy=0, sx2 = 0, sy2=0, sxy = 0;
for (int i = 0; i < clusters.length; i++) {
if (clusters[i] == center_cluster) {
double y = i / size - (yc + cent_xy[1]); //(yc + cent_xy[1]) - absolute, from (0,0)
double x = i % size - (xc + cent_xy[0]);
double w = data[i]-thresh;
s0 += w;
sx += w * x;
sy += w * y;
sx2 += w * x * x;
sy2 += w * y * y;
sxy += w * x * y;
}
}
double cxx = sx2 - sx * sx / s0, cyy= sy2 - sy * sy / s0, cxy = sxy - sx * sy / s0;
/*
* sum(Mi*(Xi-avg(X))^2) = SX2 - SX^2/S0
* sum(Mi*(Yi-avg(Y))^2) = SY2 - SY^2/S0
* sum(Mi*(Xi-avg(X)*(Yi-avg(Y))) = SXY - SX*SY / S0
*/
Matrix covar = new Matrix(new double[][] {{cxx, cxy},{cxy,cyy}});
double [] cent_offs = {sx/s0 + xc,sy/s0 + yc};
EigenvalueDecomposition eig = covar.eig();
double [] eigval = {eig.getD().get(0, 0),eig.getD().get(1, 1)};
Arrays.sort(eigval); // ascending
double elong = Math.sqrt(eigval[1]/eigval[0]);
double eff_rad = Math.sqrt(Math.sqrt(eigval[1]*eigval[0]));
double dist = Math.sqrt(cent_offs[0]*cent_offs[0] + cent_offs[1]*cent_offs[1]);
double max_other = 0.0;
if (other_radius > 0) { // only calculate if asked for
double or2 = other_radius*other_radius;
min_x = Math.max(1, (int) Math.floor(best_x-other_radius));
min_y = Math.max(1, (int) Math.floor(best_y-other_radius));
max_x = Math.min(size-2, (int) Math.ceil(best_x+other_radius));
max_y = Math.min(size-2, (int) Math.ceil(best_y+other_radius));
for (int y = min_y; y <= max_y; y++) { // do not search on very edges
double dy = (y-best_y);
double y2 = dy*dy;
if (y2 < or2) {
for (int x = min_x; x <= max_x; x++) {
double dx = x - best_x;
double r2 = y2 + dx*dx;
if (r2 < or2) {
int indx = y * size + x;
if (clusters[indx] != center_cluster) { // only disconnected from the central area
double d = data[indx];
if (d > max_other) {
max_other = d;
}
}
}
}
}
}
max_other /= best_d;
}
if (debug) {
System.out.println("\ncenter offset ["+(sx/s0)+","+(sy/s0)+"] , from center: ["+cent_offs[0]+","+cent_offs[1]+"]");
System.out.println("Covariance matrix:");
covar.print(8, 6);
System.out.println("eig.getV()");
eig.getV().print(8, 6);
System.out.println("eig.getD()");
eig.getD().print(8, 6);
System.out.println(String.format("best_d=%7.5f, rad=%6.3f, elong=%6.3f, dist=%6.3f, dx=%6.3f, dy=%6.3f, mo=%5.3f",
best_d, eff_rad, elong, dist, cent_offs[0], cent_offs[1], max_other));
}
return new double [] {best_d, eff_rad, elong, dist, cent_offs[0], cent_offs[1], max_other} ;
}
public static double [] getMaxLocAreaInt(
double [] data, // square data
double [] cent_xy, // if null, use center of the square
double radius, // search for maximum within this radius
double frac_max) {
int size = (int) Math.sqrt(data.length);
if (cent_xy == null) {
cent_xy = new double [] {size/2, size/2};
}
double maxr2 = radius * radius;
double best_d = 0;
int best_indx = -1;
int min_x = Math.max(1, (int) Math.floor(cent_xy[0]-radius));
int min_y = Math.max(1, (int) Math.floor(cent_xy[1]-radius));
int max_x = Math.min(size-2, (int) Math.ceil(cent_xy[0]+radius));
int max_y = Math.min(size-2, (int) Math.ceil(cent_xy[1]+radius));
for (int y = min_y; y <= max_y; y++) { // do not search on very edges
double dy = (y-cent_xy[1]);
double y2 = dy*dy;
if (y2 < maxr2) {
for (int x = min_x; x <= max_x; x++) {
double dx = x - cent_xy[0];
double r2 = y2 + dx*dx;
if (r2 < maxr2) {
int indx = y * size + x;
double d = data[indx];
if (d > best_d) {
best_indx = indx;
best_d = d;
}
}
}
}
}
// is it local max?
if ((data[best_indx - 1] > best_d) || (data[best_indx + 1] > best_d) ||
(data[best_indx - size] > best_d) || (data[best_indx + size] > best_d)) {
return null; // on the edge, not a local max
}
boolean [] above_thresh = new boolean [data.length];
double thresh = best_d * frac_max;
for (int i = 0; i < data.length; i++) {
above_thresh[i] = data[i] > thresh;
}
int [] clusters = (new TileNeibs(size,size)).enumerateClusters(
above_thresh, // boolean [] tiles,
null, // int [] num_clusters,
false); // boolean ordered)
int center_cluster = clusters[best_indx];
double area_max = 0.0;
for (int i = 0; i < clusters.length; i++) {
if (clusters[i] == center_cluster) {
area_max += 1;
}
}
double x_max = best_indx % size - cent_xy[0];
double y_max = best_indx / size - cent_xy[1];
double dist = Math.sqrt(y_max * y_max + x_max * x_max);
return new double [] {best_d, area_max, dist, x_max, y_max} ;
}
public static double [] getPatternCenter(
double [] data_in,
......@@ -429,6 +196,7 @@ public class ObjectLocation {
corr_pattern, // final double [] pattern, // [psize*psize]
false, // final boolean convolve, // convolve, not correlate
phaseCoeff, // final double phaseCoeff,
0, // final double lpf_sigma, // 0 - do not filter
debugLevel); // final int debugLevel)
if (corr_ret != null) {
System.arraycopy(corr_out, 0, corr_ret, 0, corr_out.length);
......
src/main/java/com/elphel/imagej/orthomosaic/OrthoMap.java
View file @ 640b8df3
......@@ -3009,6 +3009,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
* @param num_refine number of the plane generation refine passes
* @param frac_above remove high fraction of all tiles
* @param frac_below remove low fraction of all tiles
* @param tile_plane_elev null or double[2][] will return per-tile elevations
* @param debug_title Name of the debug image or null to skip
*
* @return plane parameters - all metric, relative to the vertical point:
......@@ -3023,6 +3024,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
int num_refine,
double frac_above,
double frac_below,
double [] tile_plane_elev, // null or double[2][] will return per-tile elevations
String debug_title) {
final int num_bins = 1024;
int diff_zoom_lev = zoom_lev - orig_zoom_level;
......@@ -3062,6 +3064,16 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
width, // final int width,
pix_vertical); // final double [] xy0)
}
if (tile_plane_elev != null) {
for (int pix = 0; pix < elev.length; pix++) {
int px = pix % width;
int py = pix / width;
double x = px - plane_pix[3];
double y = py - plane_pix[4];
double plane_elev = plane_pix[2]+ x *plane_pix[0] +y * plane_pix[1];
tile_plane_elev[pix] = elev [pix] - plane_elev;
}
}
if (debug_title != null) {
String [] dbg_titles = {"elev", "plane", "diff"};
double [][] dbg_img = new double [dbg_titles.length][elev.length];
......@@ -3082,7 +3094,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
debug_title,
dbg_titles);
}
return new double [] {plane_pix[0]*pix_size, plane_pix[1]*pix_size, plane_pix[2]}; // no xy0
return new double [] {plane_pix[0]/pix_size, plane_pix[1]/pix_size, plane_pix[2]}; // no xy0
}
......@@ -3259,7 +3271,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
double [] pix_vertical = {vert_meters[0]/pix_size,vert_meters[1]/pix_size}; // image center in pixels
double max_r_diff = metric_error * camera_height / pix_size;
double [] plane_pix =
{plane_metric[0]/pix_size,plane_metric[1]/pix_size,plane_metric[2],pix_vertical[0],pix_vertical[1]};
{plane_metric[0]*pix_size,plane_metric[1]*pix_size,plane_metric[2],pix_vertical[0],pix_vertical[1]};
boolean [] mask = removeHighProjection (
elev, // final double [] data,
null, // final boolean [] mask_in, // will be modified if provided
......@@ -3526,6 +3538,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
final double [] pattern, // [psize*psize]
final boolean convolve, // convolve, not correlate
final double phaseCoeff,
final double lpf_sigma, // 0 - do not filter
final int debugLevel) {
final int height = data.length/width;
final double [] dout = new double [data.length];
......@@ -3548,15 +3561,15 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
final DoubleFHT doubleFHT0 = new DoubleFHT();
final double [] patternFD = pattern.clone();
doubleFHT0.transformPattern(patternFD);
final double [] filter = null; // probably not needed
if ((width== psize) && (height == psize)) {
return testPhaseCorr(
return correlateWithPatternSquare(
data, // double [] data,
pattern, // double [] pattern,
convolve, //final boolean convolve, // convolve, not correlate
window, // double [] wnd,
phaseCoeff); // double phaseCoeff)
phaseCoeff, // double phaseCoeff)
lpf_sigma); // double lpf_sigma)
}
final int tilesX = (int) Math.ceil(width/(psize/2)) + 1;
......@@ -3581,6 +3594,10 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
double [] dtile = new double [pattern.length];
TileNeibs tn = new TileNeibs(psize,psize);
DoubleFHT doubleFHT = new DoubleFHT();
double [] filter = doubleFHT.getFrequencyFilter (
dtile, // double [] data, // just to find out size
0, // double highPassSigma,
lpf_sigma); // double lowPassSigma
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;
......@@ -3675,12 +3692,13 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
return dout;
}
public static double [] testPhaseCorr(
public static double [] correlateWithPatternSquare(
double [] data_in,
double [] pattern_in,
boolean convolve, // convolve, not correlate
double [] wnd_in,
double phaseCoeff) {
double phaseCoeff,
double lpf_sigma) {
// phaseCoeff = 0.5;
boolean dbg=false;
double [] data = data_in.clone();
......@@ -3703,6 +3721,10 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
rslt_titles);
}
DoubleFHT doubleFHT = new DoubleFHT();
double [] filter = doubleFHT. getFrequencyFilter (
data, // double [] data, // just to find out size
0, // double highPassSigma,
lpf_sigma); // double lowPassSigma
double [] data_orig2 = data.clone();
double [] pattern_orig = pattern.clone();
boolean use_new = true; // false;
......@@ -3714,14 +3736,14 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
corr_out = doubleFHT.convolvePattern (
data, // double [] first,
patternFD, // double [] secondFD,
null, // double [] filter, // high/low pass filtering
filter, // 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
filter, // double [] filter, // high/low pass filtering
null); // double [] first_save ) //null-OK
}
} else {
......@@ -3729,13 +3751,13 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
corr_out = doubleFHT. convolve ( // new
data, // double [] first,
pattern, // double [] second,
null);// double [] filter, // high/low pass filtering
filter); // 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
filter, // double [] filter, // high/low pass filtering
null, // double [] first_save,
null); // double [] second_save ) //null-OK
}
......@@ -3964,6 +3986,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
corr_patterns[n], //final double [] pattern, // [psize*psize]
false, // final boolean convolve, // convolve, not correlate
phaseCoeff, // final double phaseCoeff,
0, // final double lpf_sigma, // 0 - do not filter
debugLevel); // final int debugLevel) {
}
String [] patt_titles = new String[corrs_out.length];
......@@ -4056,6 +4079,7 @@ public class OrthoMap implements Comparable <OrthoMap>, Serializable{
corr_patterns[n], // final double [] pattern, // [psize*psize]
true, // final boolean convolve, // convolve, not correlate
phaseCoeff, // final double phaseCoeff,
0, // final double lpf_sigma, // 0 - do not filter
debugLevel); // final int debugLevel) {
}
double [] convolved_merged = mergePattenConvolutions(
......
src/main/java/com/elphel/imagej/orthomosaic/OrthoMapsCollection.java
View file @ 640b8df3
This source diff could not be displayed because it is too large. You can view the blob instead.
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