detecting peak near interference

src/main/java/com/elphel/imagej/cuas/CuasMotionLMA.java

@@ -50,6 +50,7 @@ public class CuasMotionLMA {
 	
 	private int          width;
 	private double [][]  window;
+	private double       pedestal;
 	private double []    y_vector;
 	private double []    weights;
 	private double []    full_vector = new double[INDX_LEN];
@@ -67,15 +68,40 @@ public class CuasMotionLMA {
 			double     sigma,
 			double     wnd_pedestal) {
 		this.width = width;
+		this.pedestal = wnd_pedestal;
 		window = new double [width][width];
 		double k = -0.5/(sigma*sigma);
 		for (int i = 0; i < width; i++) {
 			for (int j = 0; j < width; j++) {
-				window[i][j] = Math.exp(k*(i*i+j*j)) + wnd_pedestal; // will be normalized when used
+				window[i][j] = Math.exp(k*(i*i+j*j));//  + wnd_pedestal; // will be normalized when used
 			}
 		}
 	}
 	
+	/**
+	 * Multiply (in place) data by window. Used only when stray data (local maximum is not the absolute maximum) is present
+	 * Center data will be multiplied by almost 1.0
+	 * @param tile_data data array, will be modified
+	 * @param xc relative to center =width/2
+	 * @param yc relative to center =width/2
+	 */
+	public void applyWindowToData(
+			double []  tile_data,
+			double     xc, // relative to center =width/2
+			double     yc) { // relative to center =width/2
+		double x0 = Math.min(Math.max(xc + width/2, 0), width-1);
+		double y0 = Math.min(Math.max(yc + width/2, 0), width-1);
+		int ix0 = (int) Math.round(x0);
+		int iy0 = (int) Math.round(y0);
+		for (int y = 0; y < width; y++) {
+			int ay = Math.abs(y-iy0);
+			for (int x = 0; x < width; x++) {
+				int ax = Math.abs(x-ix0);
+				double w = window[ay][ax]; // window to the nearest integer x,y
+				tile_data[x + y*width] *= w;
+			}
+		}
+	}
 	
 	public int prepareLMA(
 			boolean [] param_select,
@@ -103,7 +129,7 @@ public class CuasMotionLMA {
 			int ay = Math.abs(y-iy0);
 			for (int x = 0; x < width; x++) {
 				int ax = Math.abs(x-ix0);
-				double w = window[ay][ax]; // window to the nearest integer x,y
+				double w = window[ay][ax]+pedestal; // window to the nearest integer x,y
 				weights[x + y*width] = w;
 				sw += w;
 			}

src/main/java/com/elphel/imagej/tileprocessor/Correlation2d.java

@@ -2412,7 +2412,6 @@ public class Correlation2d {
 		return rslt;
 	}
 
-	
 	public static double [] getMaxXYCm(
 			double []  data, // will be modified if fpn_mask != null;
 			int        data_width,      //  = 2 * transform_size - 1; // negative - calculate fraction
@@ -2420,28 +2419,64 @@ public class Correlation2d {
 			int        refine, //  re-center window around new maximum. 0 -no refines (single-pass)
 			boolean [] fpn_mask,
 			boolean    ignore_border, // only if fpn_mask != null - ignore tile if maximum touches fpn_mask
-			boolean    debug)
-	{
+			boolean    debug) {
+		
 		boolean exclude_margins = false;
 		if ((fpn_mask!=null) && (fpn_mask.length==0)) {
 			exclude_margins = true;
 			fpn_mask = null;
 		}
+		
 		boolean calc_fraction = data_width < 0;
 		if (calc_fraction) {
 			data_width = -data_width;
 		}
+		int imax= 0;
+		for (int i= 1; i < data.length;i++) {
+			if (data[i] > data[imax]) {
+				imax = i;
+			}
+		}
+		return getMaxXYCm(
+				data,            // double []  data, // will be modified if fpn_mask != null;
+				data_width,      // int        data_width,      //  = 2 * transform_size - 1; // negative - calculate fraction
+				 radius,         // double     radius, // 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
+				refine,          // int        refine, //  re-center window around new maximum. 0 -no refines (single-pass)
+				fpn_mask,        // boolean [] fpn_mask,
+				ignore_border,   // boolean    ignore_border, // only if fpn_mask != null - ignore tile if maximum touches fpn_mask
+				exclude_margins, // boolean    exclude_margins,
+				calc_fraction,   // boolean    calc_fraction,
+				imax,            // int        pix_max,
+				debug);          // boolean    debug);
+	}
+
+	
+	public static double [] getMaxXYCm(
+			double []  data, // will be modified if fpn_mask != null;
+			int        data_width,      //  = 2 * transform_size - 1; // negative - calculate fraction
+			double     radius, // 0 - all same weight, > 0 cosine(PI/2*sqrt(dx^2+dy^2)/rad)
+			int        refine, //  re-center window around new maximum. 0 -no refines (single-pass)
+			boolean [] fpn_mask,
+			boolean    ignore_border, // only if fpn_mask != null - ignore tile if maximum touches fpn_mask
+			boolean    exclude_margins,
+			boolean    calc_fraction,
+			int        imax,  // index of the maximum in data[]
+			boolean    debug)
+	{
 		int data_height = data.length/data_width;
-//		int center_xy = (data_width - 1)/2; //  = transform_size - 1;
-		int center_xy = data_width / 2; //  = transform_size - 1;
+		int center_xy =   data_width / 2; //  = transform_size - 1;
 		double x0 = center_xy, y0 = center_xy;
+		/*
 		int imax= 0;
 		for (int i= 1; i < data.length;i++) {
 			if (data[i] > data[imax]) {
 				imax = i;
 			}
 		}
+		*/
+		
 		double mx = data[imax];
+		
 		if (mx < 0) {
 			return null; //negative maximum
 		}

src/main/java/com/elphel/imagej/tileprocessor/IntersceneMatchParameters.java

@@ -723,6 +723,12 @@ min_str_neib_fpn	0.35
     public double   cuas_speed_pref =    0.0;     // preferable speed (boost weights for faster targets)
     public double   cuas_speed_boost =   1.0;     // speed boost limit
 	// target filtering after constant velocity accumulation
+    
+    public double   cuas_lmax_fraction = 0.7;     // Check if local maximum is separated from tye surrounding by this fraction of the maximum value 
+    public double   cuas_lmax_radius =   3.5;     // look inside ((int)cuas_lmax_radius) * 2 + 1 square for the local maximum isolation 
+	public boolean  cuas_lmax_zero =     true;    // zero all data outside this radius from the maximum
+    
+    
 	public double   cuas_target_radius = 3.0;     // target centroids center radius
 	public double   cuas_target_strength =0.8;    // target centroids center radius
 	public double [][] cuas_target_frac = {{0,0.15},{2.5,0.18},{5,0.3}};
@@ -2242,6 +2248,14 @@ min_str_neib_fpn	0.35
 				"Boost effective strength when speed is above this.");
 		gd.addNumericField("Maximal speed boost",                    this.cuas_speed_boost, 5,8,"",
 				"Maximal speed-caused effective strength boost.");
+		
+		gd.addNumericField("Maximums separation fraction",           this.cuas_lmax_fraction, 5,8,"",
+				"Check if local maximum is separated from tye surrounding by this fraction of the maximum value.");
+		gd.addNumericField("Maximums separation radius",             this.cuas_lmax_radius, 5,8,"pix",
+				"Look inside ((int)cuas_lmax_radius) * 2 + 1 square for the local maximum isolation .");
+		gd.addCheckbox    ("Zero far from maximums pixels",          this.cuas_lmax_zero,
+				"Zero all data outside this radius from the maximum.");
+		
 		gd.addNumericField("Target radius",                          this.cuas_target_radius, 5,8,"pix",
 				"Target radius, also used to calculate fraction of totals inside (windowed) to all positive values.");
 		gd.addNumericField("Minimal target strength",                this.cuas_target_strength, 5,8,"",
@@ -3355,6 +3369,10 @@ min_str_neib_fpn	0.35
 		this.cuas_speed_pref =          gd.getNextNumber();
 		this.cuas_speed_boost =         gd.getNextNumber();
 		
+		this.cuas_lmax_fraction =       gd.getNextNumber();
+		this.cuas_lmax_radius =         gd.getNextNumber();
+		this.cuas_lmax_zero =           gd.getNextBoolean();   
+		
 		this.cuas_target_radius =       gd.getNextNumber();
 		this.cuas_target_strength =     gd.getNextNumber();
 		this.cuas_target_frac = stringToDouble2d(gd.getNextString());
@@ -4324,6 +4342,10 @@ min_str_neib_fpn	0.35
 		properties.setProperty(prefix+"cuas_speed_pref",      this.cuas_speed_pref+"");     // double
 		properties.setProperty(prefix+"cuas_speed_boost",     this.cuas_speed_boost+"");    // double
 		
+		properties.setProperty(prefix+"cuas_lmax_fraction",   this.cuas_lmax_fraction+"");  // double
+		properties.setProperty(prefix+"cuas_lmax_radius",     this.cuas_lmax_radius+"");    // double
+        properties.setProperty(prefix+"cuas_lmax_zero",       this.cuas_lmax_zero+"");      // boolean
+		
 		properties.setProperty(prefix+"cuas_target_radius",   this.cuas_target_radius+"");  // double
 		properties.setProperty(prefix+"cuas_target_strength", this.cuas_target_strength+"");// double
 		properties.setProperty(prefix+"cuas_target_frac", double2dToString(cuas_target_frac)+"");  // double[][]
@@ -5260,6 +5282,10 @@ min_str_neib_fpn	0.35
 		if (properties.getProperty(prefix+"cuas_speed_pref")!=null)      this.cuas_speed_pref=Double.parseDouble(properties.getProperty(prefix+"cuas_speed_pref"));
 		if (properties.getProperty(prefix+"cuas_speed_boost")!=null)     this.cuas_speed_boost=Double.parseDouble(properties.getProperty(prefix+"cuas_speed_boost"));
 		
+		if (properties.getProperty(prefix+"cuas_lmax_fraction")!=null)   this.cuas_lmax_fraction=Double.parseDouble(properties.getProperty(prefix+"cuas_lmax_fraction"));
+		if (properties.getProperty(prefix+"cuas_lmax_radius")!=null)     this.cuas_lmax_radius=Double.parseDouble(properties.getProperty(prefix+"cuas_lmax_radius"));
+		if (properties.getProperty(prefix+"cuas_lmax_zero")!=null)       this.cuas_lmax_zero=Boolean.parseBoolean(properties.getProperty(prefix+"cuas_lmax_zero"));
+		
 		if (properties.getProperty(prefix+"cuas_target_radius")!=null)   this.cuas_target_radius=Double.parseDouble(properties.getProperty(prefix+"cuas_target_radius"));
 		if (properties.getProperty(prefix+"cuas_target_strength")!=null) this.cuas_target_strength=Double.parseDouble(properties.getProperty(prefix+"cuas_target_strength"));
 		if (properties.getProperty(prefix+"cuas_target_frac")!=          null) cuas_target_frac=stringToDouble2d((String) properties.getProperty(prefix+"cuas_target_frac"));
@@ -6196,6 +6222,11 @@ min_str_neib_fpn	0.35
 		imp.cuas_speed_min =        this.cuas_speed_min;
 		imp.cuas_speed_pref =       this.cuas_speed_pref;
 		imp.cuas_speed_boost =      this.cuas_speed_boost;
+		
+		imp.cuas_lmax_fraction =    this.cuas_lmax_fraction;
+		imp.cuas_lmax_radius =      this.cuas_lmax_radius;
+		imp.cuas_lmax_zero =        this.cuas_lmax_zero;
+		
 		imp.cuas_target_frac =      new double [this.cuas_target_frac.length][];
 		for (int i = 0; i < this.cuas_target_frac.length; i++) {
 			imp.cuas_target_frac[i] = this.cuas_target_frac[i].clone();

src/main/java/com/elphel/imagej/tileprocessor/TileNeibs.java

 package com.elphel.imagej.tileprocessor;
+import java.awt.Point;
 import java.awt.Rectangle;
 import java.util.ArrayList;
 import java.util.Arrays;
@@ -45,26 +46,47 @@ public class TileNeibs{
 	final public static int DIRS =        DIR_XY.length; //8; // total dirs
 	final static int THREADS_MAX =        100;
 	
+	final int              sizeX;
+	final int              sizeY;
+	public Rectangle       inner = null;
+	public Rectangle       full = null;
+
+	public Rectangle getFullRectangle() {
+		if (full == null) {
+			full = new Rectangle(0,0,sizeX,sizeY);
+		}
+		return full;
+	}
+	
+	public Rectangle getInnerRectangle() {
+		if (inner == null) {
+			inner = new Rectangle(1,1,sizeX-2, sizeY-2);
+		}
+		return inner;
+	}
+	
+	int last_grown; // actual grown until died
+	public int dirs = DIRS;
+
 	public static int reverseDir(int dir) {
 		if ((dir < 0) || (dir >= DIRS)) {
 			return dir;
 		}
 		return (dir+DIRS/2) % DIRS;
 	}
-	int sizeX;
-	int sizeY;
-	int last_grown; // actual grown until died
-	public int dirs = DIRS;
+	
 	public int getLastGrown() {
 		return last_grown;
 	}
 	public TileNeibs(int size){
 		this.sizeX = size;
 		this.sizeY = size;
+		this.inner = new Rectangle(1,1,size-2, size-2);
 	}
 	public TileNeibs(int sizeX, int sizeY){
 		this.sizeX = sizeX;
 		this.sizeY = sizeY;
+		this.inner = new Rectangle(1,1,sizeX-2, sizeY-2);
 	}
 	public int numNeibs() // TODO: make configurable to
 	{
@@ -100,6 +122,10 @@ public class TileNeibs{
 	
 	int getY(int indx) {return indx / sizeX;};
 	
+	Point getPoint(int indx) {
+		return new Point (indx % sizeX, indx / sizeX);
+	}
+	
 	public int getSizeX() {
 		return sizeX;
 	}
@@ -1412,4 +1438,116 @@ public class TileNeibs{
 		return filled;
 	}
 	
+	/**
+	 * Get array of indices for all local maximums (pixels that are not less than their 8 neighbors). NaNs are considered smaller than normal values.
+	 * @param data data in line-scan order, should be [sizeX*sizeY] (not verified)
+	 * @param exclude_margins do not allow maximums on the outer edge of the defined rectangle 
+	 * @return variable-length interger array of local maxima
+	 */
+	public int [] getLocalMaxes(
+			double [] data,
+			boolean   exclude_margins) {
+		ArrayList<Integer> max_list = new ArrayList<Integer>();
+		int x0 = exclude_margins ? 1 : 0;
+		int y0 = exclude_margins ? 1 : 0;
+		int x1 = sizeX - x0 - 1;
+		int y1 = sizeY - y0 - 1;
+		for (int y = y0; y <= y1; y++) {
+			for (int x = x0; x <= x1; x++) {
+				int ntile = y * sizeX + x;
+				double d = data[ntile];
+				if (!Double.isNaN(d)) {
+					check_ismax: {
+						for (int dir = 0; dir < DIRS; dir++) {
+							int ntile1 = getNeibIndex(ntile,dir);
+							if ((ntile >=0) && (data[ntile1] > d)){
+								break check_ismax;
+							}
+						}
+						max_list.add(ntile);
+					}
+				}
+			}
+		}
+        return max_list.stream().mapToInt(Integer::intValue).toArray();
+	}
+
+	/**
+	 * Calculate index of an absolute maximum in array, ignore NaNs. If all elements are NaN, will return -1
+	 * @param data Data array, may contain NaN values
+	 * @return index of a largest element (NaN is considered smaller than any normal value) or -1 if all are NaN
+	 */
+	public static int getAmaxTile(
+			double [] data) {
+		int ntile = - 1;
+		double maxv = Double.NaN;
+		for (int i = 0; i < data.length; i++) {
+			double d = data[i];
+			if (!(d <= maxv)) {
+				ntile = i;
+				maxv = d;
+			}
+		}
+		return ntile;
+	}
+
+	/**
+	 * Verify that this local maximum is disconnected from any tile larger it by the tiles having
+	 * value of the specified fraction of the maximal value. 
+	 * @param data data array
+	 * @param ntile tile index
+	 * @param fraction fraction of the maximal value that separates the maximum
+	 * @param radius if >0, limits allowable expansion (and returns false if reached)
+	 * @return true if it is an isolated maximum, false if not.
+	 */
+	public boolean isMaxIsolated(
+			double [] data,
+			int       ntile,
+			double    fraction,
+			int       radius) {
+		if ((ntile <0) || (ntile >= data.length)) {
+			return false; // outside
+		}
+		int xc = getX(ntile);
+		int yc = getY(ntile);
+		Rectangle full = getFullRectangle();
+		Rectangle limits = (radius <=0) ?  null : (full.intersection(new Rectangle(xc-radius, yc-radius, 2*radius+1,2*radius+1)));
+		boolean [] tried = new boolean [data.length];
+		double local_max = data[ntile];
+		double threshold = local_max * fraction;
+		ArrayList<Integer> wave = new ArrayList<Integer>();
+		wave.add(ntile);
+		tried[ntile] = true;
+		while (!wave.isEmpty()) {
+			ntile = wave.remove(0);
+			for (int dir = 0; dir < DIRS; dir++) {
+				int ntile1 = getNeibIndex(ntile, dir);
+				if ((ntile1 >= 0) && !tried[ntile1]) {
+					tried[ntile1] = true;
+					double d = data[ntile1]; // may be NaN
+					if (d > local_max) {
+						return false;
+					} else if (d > threshold) { // NaN will result in false
+						wave.add(ntile1);
+						if ((limits != null) && !limits.contains(getPoint(ntile1))) {
+							return false;
+						}
+					}
+				}
+			}
+			
+		}
+		return true;
+	}
+	
+	/**
+	 * Check if the point indx is on the edge of the defined rectangle
+	 * @param indx
+	 * @return
+	 */
+	public boolean isEdge(int indx) {
+		return !getFullRectangle().contains(getPoint(indx));
+	}
+	
+	
 }