HatSnakeFilter.java
package com.github.celldynamics.quimp.geom.filters;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.nio.file.Path;
import java.nio.file.Paths;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.TreeSet;
import java.util.stream.Collectors;
import java.util.stream.IntStream;
import org.apache.commons.lang3.tuple.MutablePair;
import org.apache.commons.lang3.tuple.Pair;
import org.scijava.vecmath.Point2d;
import org.scijava.vecmath.Tuple2d;
import org.scijava.vecmath.Vector2d;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import com.github.celldynamics.quimp.Outline;
import com.github.celldynamics.quimp.QuimP;
import com.github.celldynamics.quimp.geom.BasicPolygons;
import com.github.celldynamics.quimp.plugin.QuimpPluginException;
import com.github.celldynamics.quimp.plugin.ecmm.ODEsolver;
import com.github.celldynamics.quimp.plugin.utils.IPadArray;
import com.github.celldynamics.quimp.plugin.utils.QuimpDataConverter;
import ij.IJ;
import ij.process.ImageProcessor;
/**
* Implementation of HatFilter for removing convexes from polygon.
*
* <h3>List of user parameters:</h3>
* <ol>
* <li><i>window</i> - Size of window in pixels. It is responsible for sensitivity to protrusions of
* given size. Larger window can eliminate small and large protrusions whereas smaller window is
* sensitive only to small protrusions. Window size should be from 3 to number of outline points.
* <li><i>pnum</i> - Number of protrusions/cavities that will be removed from outline. If not
* limited by <i>alev</i> parameter the algorithm will eliminate <i>pnum</i> features
* (convex/concave parts) from outline. Feature is removed if its rank is bigger than <i>alev</i>
* threshold. If <i>pnum</i> is set to 0 algorithm will remove all candidates with rank higher than
* <i>alev</i> regardless their number. <i>pnum</i> should be from 0 to any value. Algorithm stops
* searching when there is no candidates to remove.
* <li><i>alev</i> - Threshold value, if circularity computed for given window position is lower
* than threshold this window is not eliminated regarding <i>pnum</i> or its rank in circularities.
* <i>alev</i> should be in range form 0 to ?, where 0 stands for accepting every candidate
* </ol>
*
* <p><h3>General description of algorithm:</h3> The window slides over the wrapped contour. Points
* inside window for its position <i>p</i> are considered as candidates to removal from contour if
* they meet the following criterion:
* <ol>
* <li>The window has achieved for position p circularity parameter <i>c</i> larger than
* <i>alev</i>
* <li>The window on position <i>p</i> does not touch any other previously found window.
* <li>Points of window <i>p</i> are convex/concave (can be configured by {@link #setMode(int)}).
* </ol>
*
* <p>Every window <i>p</i> has assigned a <i>rank</i>. Bigger <i>rank</i> stands for better
* candidate
* to remove. Algorithm tries to remove first <i>pnum</i> windows (those with biggest ranks) that
* meet above rules. If <i>pnum</i> is set to 0 all protrusions with rank > <i>alev</i> are deleted.
*
* <p><H3>Detailed description of algorithm</H3> The algorithm comprises of three main steps:
* <ol>
* <li>Preparing <i>rank</i> table of candidates to remove
* <li>Iterating over <i>rank</i> table to find <i>pnum</i> such candidates who meet rules and store
* their coordinates in <i>ind2rem</i> array. By candidates it is understood sets of polygon indexes
* that is covered by window on given position. For simplification those vertexes are identified by
* lover and upper index of window in outline array (input). <i>pnum</i> can be 0, see note above.
* <li>Forming output table without protrusions.
* </ol>
*
* <p><H2>First step</H2> The window of size <i>window</i> slides over looped data. Looping is
* performed by {@link Collections#rotate(List, int)} method that shift data left copying
* falling out indexes
* to end of the set (finally the window is settled in constant position between indexes
* <0;window-1>). For each its position <i>r</i> the candidate points are deleted from original
* contour and circularity is computed (see {@link #getCircularity(List)}). Then candidate points
* are passed to {@link #getWeighting(List)} method where weight is evaluated. The role of weight is
* to promote in <i>rank</i> candidate points that are cumulated in small area over distributed
* sets. Thus weight should give larger values for that latter distribution than for cumulated one.
* Currently weights are calculated as squared standard deviation of distances of all candidate
* points to
* center of mass of these points (or mean point if polygon is invalid). There is also mean
* intensity calculated by {@link #getIntensity(List, List, ImageProcessor)} (not obligatory, it can
* be switched off using {@link #runPlugin(List)} method.). It is sampled along
* shrunk outline produced from input points. Finally
* circularity(r) is divided by weight (<i>r</i>) and intensity and stored in <i>circ</i> array.
* Additionally in this step the convex/concave is checked. All candidate points are tested whether
* <b>all</b> they are inside the contour made without them (for concave option) or whether
* <b>any</b> of
* them is inside the modified contour (for convex option). Convex/concave sensitivity option can be
* modified by {@link #setMode(int)}
* This information is stored in <i>convex</i> array. Finally rank array <i>circ</i> is NOT
* normalised.
*
* <p><H2>Second step</H2> In second step array of ranks <i>circ</i> is sorted in descending order.
* For every rank in sorted table the real position of window is retrieved (that gave this rank).
* The window position is defined here by two numbers - lover and upper range of indexes
* covered by it. The candidate points from this window are validated for criterion:
* <ol>
* <li><i>rank</i> must be greater than <i>alev</i>
* <li>lower and upper index of window (index means here number of polygon vertex in array) must not
* be included in any previously found window. This checking is done by deriving own class
* WindowIndRange with overwritten {@link HatSnakeFilter.WindowIndRange#compareTo(Object)} method
* that defines
* rules of equality and non relations between ranges. Basically any overlapping range or included
* is considered as equal and rejected from storing in <i>ind2rem</i> array.
* <li>candidate points must be convex or concave or any.
* </ol>
* If all above criterion are meet the window (l;u) is stored in <i>ind2rem</i>. Windows on the end
* of data are wrapped by dividing them for two sub-windows: (w;end) and (0;c) otherwise they may
* cover the whole range (e.g. <10;3> does not stand for window from 10 wrapped to 3 but window from
* 3 to 10).
*
* <p>The second step is repeated until <i>pnum</i> object will be found or end of candidates will
* be reached. If <tt>pnum</tt> is 0 all objects that meet above criteria will be found.
*
* <p><H2>Third step</H2> In third step every point from original contour is tested for including in
* array <i>ind2rem</i> that contains ranges of indexes to remove. Points on index that is not
* included in any of ranges stored in <i>ind2rem</i> are copied to output.
*
* @author p.baniukiewicz
*/
public class HatSnakeFilter implements IPadArray {
static final Logger LOGGER = LoggerFactory.getLogger(HatSnakeFilter.class.getName());
/**
* Algorithm will look for cavities.
*
* @see #setMode(int)
*/
public static final int CAVITIES = 1;
/**
* Algorithm will look for protrusions.
*
* @see #setMode(int)
*/
public static final int PROTRUSIONS = 2;
/**
* Algorithm will look for cavities and protrusions.
*
* @see #setMode(int)
*/
public static final int ALL = 3;
/**
* Number of steps used for outline shrinking for intensity sampling.
*/
public int shrinkAmount = 3;
/**
* Set it to 1 to look for inclusions, 2 for protrusions, 3 for all.
*/
private int lookFor = PROTRUSIONS; // default for compatibility with old code
private int window; // filter's window size
private int pnum; // how many protrusions to remove
private double alev; // minimal acceptance level
private int debugCounter = 0; // for naming of debug files
/**
* Construct HatFilter Input array with data is virtually circularly padded.
*/
public HatSnakeFilter() {
this.window = 15;
this.pnum = 1;
this.alev = 0;
LOGGER.debug("Set default parameter: window=" + window + " pnum=" + pnum + " alev=" + alev);
}
/**
* Constructor passing algorithm parameters.
*
* @param window size of the window
* @param pnum number of protrusions to find
* @param alev threshold
* @see HatSnakeFilter
*/
public HatSnakeFilter(int window, int pnum, double alev) {
super();
this.window = window;
this.pnum = pnum;
this.alev = alev;
}
/**
* In contrary to {@link #runPlugin(List, ImageProcessor, Pair)} this method does not use
* intensity weighting nor externally calculated ranks.
*
* @param data contour to process
* @return Processed input list, size of output list may be different than input. Empty output
* is also allowed.
* @throws QuimpPluginException on wrong input data
* @see #runPlugin(List, ImageProcessor, Pair)
* @see #runPlugin(List, ImageProcessor)
*/
public List<Point2d> runPlugin(List<Point2d> data) throws QuimpPluginException {
return runPlugin(data, null);
}
/**
* In contrary to {@link #runPlugin(List, ImageProcessor, Pair)} this method does not use
* externally calculated ranks.
*
* @param data contour to process
* @param orgIp Original image used for sampling intensity, can be <tt>null</tt>
* @return Processed input list, size of output list may be different than input. Empty output
* is also allowed.
* @throws QuimpPluginException on wrong input data
* @see #runPlugin(List)
* @see #runPlugin(List, ImageProcessor, Pair)
*/
public List<Point2d> runPlugin(List<Point2d> data, ImageProcessor orgIp)
throws QuimpPluginException {
return runPlugin(data, orgIp, null);
}
/**
* Main filter runner.
*
* <p>This version assumes that user clicked Apply button to populate data from UI to plugin or
* any other ui element. User can expect that points will be always valid but they optionally may
* have 0 length.
*
* @param points contour to process
* @param orgIp Original image used for sampling intensity, can be <tt>null</tt>
* @param ranks ranks calculated by {@link #calculateRank(List, ImageProcessor)} or <tt>null</tt>
* to let them be evaluated by this method. Ranks values must comply with <tt>alev</tt>
* value
*
* @return Processed input list, size of output list may be different than input. Empty output
* is also allowed.
* @throws QuimpPluginException on problem with input data
* @see #calculateRank(List, ImageProcessor)
*/
public List<Point2d> runPlugin(List<Point2d> points, ImageProcessor orgIp,
Pair<ArrayList<Double>, ArrayList<Boolean>> ranks) throws QuimpPluginException {
// internal parameters are not updated here but when user click apply
LOGGER.debug(String.format("Run plugin with params: window %d, pnum %d, alev %f", window, pnum,
alev));
// check input conditions
if (window % 2 == 0 || window < 0) {
throw new QuimpPluginException("Window must be uneven, positive and larger than 0");
}
if (window >= points.size()) {
throw new QuimpPluginException("Processing window to long");
}
if (window < 3) {
throw new QuimpPluginException("Window should be larger than 2");
}
if (pnum < 0) {
throw new QuimpPluginException("Number of protrusions should be larger than 0");
}
if (alev < 0) {
throw new QuimpPluginException("Acceptacne level should be positive");
}
// Step 1 - Calculate ranks and normalise it - Give up - see tag:LocalAveraging
if (ranks == null) {
ranks = calculateRank(points, orgIp);
}
ArrayList<Double> circ = ranks.getLeft(); // out ranks
ArrayList<Boolean> cavprot = ranks.getRight(); // shape flag
// Step 2 - Check criterion for all windows
if (circ == null || circ.isEmpty()) {
LOGGER.info("No candidates found for selected mode: " + lookFor);
return points; // just return non-modified data;
}
TreeSet<WindowIndRange> ind2rem = new TreeSet<>(); // <l;u> range of indexes to remove
// need sorted but the old one as well to identify windows positions
// we do not touch circ itself but rater get list of indexes that refer to sorted element in
// original array. This is sort in ascending order
IntStream tmpStr = IntStream.range(0, circ.size()).boxed()
.sorted((i, j) -> circ.get(i).compareTo(circ.get(j))).mapToInt(el -> el);
List<Integer> sortedIndexes = tmpStr.boxed().collect(Collectors.toList());
Collections.reverse(sortedIndexes); // to have in descending orer - first index in this array
// point to index in sorted aray where largest element sits
LOGGER.trace("cirI: " + sortedIndexes.toString());
LOGGER.trace("circ: " + circ.toString());
LOGGER.trace("circM: " + circ.get(sortedIndexes.get(0)).toString() + " (among ALL)");
for (int i = 0; i < sortedIndexes.size(); i++) { // for debug only
if (cavprot.get(sortedIndexes.get(i)) == true) {
LOGGER.trace(
"circMacc: " + circ.get(sortedIndexes.get(i)).toString() + " (max requested type)");
break;
}
}
if (circ.get(sortedIndexes.get(0)) < alev) {
LOGGER.info("Skipped this frame due to rank[0]=" + circ.get(sortedIndexes.get(0)));
return points; // just return non-modified data; TODO does it matter if circ is normalised?
}
int found = 0; // how many protrusions we have found already
// current index in circsorted - number of window to analyze. Real position of
// window on data can be retrieved by finding value from sorted array circularity
// in non sorted, where order of data is related to window position starting
// from 0 for most left point of window
int i = 0; // position where window begins (linear index of sortedIndexes)
boolean contains; // temporary result of test if current window is included in any prev
// temporary variable for keeping window currently tested for containing in ind2rem
WindowIndRange indexTest = new WindowIndRange();
// do as long as we can find candidates with rank higher than level and we found less than
// defined number of candidates. Latter criterion can be switched off by setting pnum to 0
// note that rank list contains all candidates and those with unwanted curvature type are
// filtered out in loop
while (i < sortedIndexes.size() && circ.get(sortedIndexes.get(i)) > alev
&& (found < pnum || pnum == 0)) {
// find where it was before sorting and store in window positions
int startpos = sortedIndexes.get(i);
if (cavprot.get(startpos) == false || lookFor == ALL) { // not valid, skip unless ALL mode
i++;
continue;
}
// check if we already have this index in list indexes to remove
if (startpos + window - 1 >= points.size()) { // if at end, we must turn to begin
indexTest.setRange(startpos, points.size() - 1); // to end
contains = ind2rem.contains(indexTest); // beginning of window at the end of dat
indexTest.setRange(0, window - (points.size() - startpos) - 1); // turn to start
contains &= ind2rem.contains(indexTest); // check rotated part at beginning
} else {
indexTest.setRange(startpos, startpos + window - 1);
contains = ind2rem.contains(indexTest);
}
if (contains == true) { // if contains go to next candidate
i++;
continue;
}
// if (cavprot.get(startpos) == true || lookFor == ALL) {
// store range of indexes that belongs to window
if (startpos + window - 1 >= points.size()) { // as prev split to two windows
// if we are on the end of data
ind2rem.add(new WindowIndRange(startpos, points.size() - 1));
// turn window to beginning
ind2rem.add(new WindowIndRange(0, window - (points.size() - startpos) - 1));
} else {
ind2rem.add(new WindowIndRange(startpos, startpos + window - 1));
}
LOGGER.info("added win for i=" + i + " startpos=" + startpos + " coord:"
+ points.get(startpos).toString() + "alev=" + circ.get(sortedIndexes.get(i)));
found++;
// }
i++;
}
if (i >= sortedIndexes.size()) { // no more data to check, probably we have less prot. pnum
LOGGER.info("Can find next candidate. You can use smaller window");
}
if (found == 0) {
LOGGER.info("No acceptable candidates found, having desired rank AND being requested type");
}
LOGGER.trace("[indexRange;Point]: " + ind2rem.toString());
LOGGER.trace("Found :" + found + " accepted windows");
// Step 3 - remove selected windows from input data
// array will be copied to new one skipping points to remove
List<Point2d> out = new ArrayList<Point2d>(); // output table for plotting temporary results
for (i = 0; i < points.size(); i++) {
// set upper and lower index to the same value - allows to test particular index for its
// presence in any defined range
indexTest.setSame(i);
if (!ind2rem.contains(indexTest)) { // check if any window position (l and u bound)
out.add(new Point2d(points.get(i)));
} else { // include tested point. Copy it to new array if not
LOGGER.trace("winpos: " + ind2rem.toString() + " " + points.get(i));
}
}
return out;
}
/**
* Evaluate rank for given outline. Return raw unnormalised values.
*
* <p>Rank is evaluated for each position of window (defined in
* {@link #HatSnakeFilter(int, int, double)} and it takes under account shape and local image
* intensity. If image is <tt>null</tt> only shape features are considered.
*
* @param points Outline as list of points
* @param orgIp Image for sampling intensity along constricted outline. Can be null
* @return Ranks (left) and cavprot flag (right). Indexes in these arrays correlate with indexes
* in input array. E.g. rank[i] stand for rank evaluated for window at position input[i]
* (points input[i] - input[i+window-1]). Shape flag is always true for selected type
* (CAVITIES or PROTRUSIONS). For ALL it is ignored.
*/
public Pair<ArrayList<Double>, ArrayList<Boolean>> calculateRank(List<Point2d> points,
ImageProcessor orgIp) {
Path tmpDebug;
PrintWriter pw = null;
List<Point2d> shCont = new ArrayList<>();
// create shrunk outline to sample intensity - one of the parameters used for candidate rank
// this is unnecessary if there is no image provided but kept here for code simplicity
Outline outline = new QuimpDataConverter(points).getOutline(); // FIXME What if more contours?
// shrink original to sample intensities close cortex - used for detecting vesicles that are
// holes in cortex area.
outline.scaleOutline(shrinkAmount, -0.3, 0.1, 0.01);
outline.unfreezeAll();
outline.correctDensity(1, 0.5); // dense shape, shrank has different number of verts than org
shCont = outline.asList();
if (QuimP.SUPER_DEBUG) {
tmpDebug = Paths.get(System.getProperty("java.io.tmpdir"),
"HatSnakeFilter_debug" + debugCounter + ".csv");
try {
pw = new PrintWriter(new FileWriter(tmpDebug.toFile()), true);
} catch (IOException e) {
e.printStackTrace();
}
Path shContDebug =
Paths.get(System.getProperty("java.io.tmpdir"), "shCont_debug" + debugCounter);
Path pointsDebug =
Paths.get(System.getProperty("java.io.tmpdir"), "points_debug" + debugCounter);
debugSaveList(shContDebug, shCont);
debugSaveList(pointsDebug, points);
debugCounter++;
}
BasicPolygons bp = new BasicPolygons(); // provides geometry processing
// Step 1 - Build circularity table
// array to store circularity for window positions. Index is related to window position
// (negative shift in rotate)
ArrayList<Double> circ = new ArrayList<Double>();
// store information if points for window at r position are convex compared to shape without
// these points
ArrayList<Boolean> cavprot = new ArrayList<Boolean>();
double tmpCirc;
double tmpInt; // mean intensity along contour in window
double tmpWei; // weighting based on local shape disturbation
for (int r = 0; r < points.size(); r++) {
// get all points except window. Window has constant position 0 - (window-1)
List<Point2d> pointsnowindow = points.subList(window, points.size());
tmpCirc = getCircularity(pointsnowindow);
// calculate weighting for circularity
List<Point2d> pointswindow = points.subList(0, window); // get points for window only
// will return 1.0 if there is no image provided
tmpInt = getIntensity(shCont, pointswindow, orgIp); // mean intensity for window points
tmpInt = tmpInt == 0.0 ? 1.0 : tmpInt; // remove 0 as we divide weight later
tmpWei = getWeighting(pointswindow);
double rank = tmpCirc / (tmpWei * tmpInt); // calculate rank for window content
circ.add(rank); // store weighted circularity for shape without window
// check if points of window are convex according to shape without these points
boolean tmpCon;
switch (lookFor) {
case CAVITIES:
tmpCon = bp.areAllPointsInside(pointsnowindow, pointswindow); // true if concave
break;
case PROTRUSIONS:
case ALL:
default:
tmpCon = bp.areAllPointOutside(pointsnowindow, pointswindow); // true if protrusions
}
// true values are those that interest us (either concave or protrusions)
// for ALL yhis list is ignored
cavprot.add(tmpCon);
// move window to next position rotates by -1 what means that on first n positions
// of points there are different values simulate window
// first iter 0 1 2 3 4 5... (w=[0 1 2])
// second itr 1 2 3 4 5 0... (w=[1 2 3])
// last itera 5 0 1 2 3 4... (w=[5 0 1])
Collections.rotate(points, -1);
Collections.rotate(shCont, -1);
// dump to file
if (QuimP.SUPER_DEBUG) {
pw.print(r + ",");
pw.print(IJ.d2s(tmpCirc, 15) + ","); // circularity only
pw.print(IJ.d2s(tmpInt, 15) + ","); // intensity
pw.print(IJ.d2s(tmpWei, 15) + ","); // shape distribution weight
pw.print(IJ.d2s(rank, 15) + ","); // final rank unscaled
pw.print(tmpCon + ","); // boolean concavity
pw.print(lookFor + ","); // boolean type of algorithm used
pw.print(pointswindow); // coordinates of points in window
pw.println();
}
}
if (QuimP.SUPER_DEBUG) {
pw.close();
}
Pair<ArrayList<Double>, ArrayList<Boolean>> ret = new MutablePair<>(circ, cavprot);
return ret;
}
/**
* Debug - saves list as tab separated coordinates.
*
* @param name name and path
* @param list list to save
*/
private void debugSaveList(Path name, List<Point2d> list) {
try {
PrintWriter pw = new PrintWriter(new FileWriter(name.toFile()), true);
int l = 0;
pw.print(list.size() + "\n");
for (Point2d p : list) {
pw.print(l + "\t" + IJ.d2s(p.getX(), 2) + "\t" + IJ.d2s(p.getY(), 2) + "\n");
l++;
}
pw.close();
} catch (IOException e) {
e.printStackTrace();
}
}
/**
* Set detection mode.
*
* @param mode can be CAVITIES, PROTRUSIONS, BOTH
*/
public void setMode(int mode) {
lookFor = mode;
}
/**
* Compute mean intensity for list of points. 3x3 stencil is used for each point.
*
* <p>Shrank contour and original contour have different number of points so direct mapping is not
* possible. For each point in window the method finds closest point in shrank contour and sample
* 3x3 stencil around. All stencils for all window points are averaged then.
*
* @param shpoints shrink contour used for sampling intensities
* @param pointswindow coordinates of window in original outline
* @param orgIp image
* @return 1.0 if input image is null. mean otherwise
*/
double getIntensity(List<Point2d> shpoints, List<Point2d> pointswindow, ImageProcessor orgIp) {
double meanI = 0.0;
if (orgIp == null) {
return 1.0; // case where intensity is not used
}
for (Point2d p : pointswindow) {
Point2d closest = findClosest(shpoints, p);
meanI += ODEsolver.sampleFluo(orgIp, (int) Math.round(closest.getX()),
(int) Math.round(closest.getY()));
}
meanI /= pointswindow.size();
return meanI;
}
/**
* Find point in list that is closest to given reference point.
*
* @param shpoints List of points to search in
* @param p Reference point
* @return Point from list that is closes to reference point
*/
private Point2d findClosest(List<Point2d> shpoints, Point2d p) {
double dist = Double.MAX_VALUE;
int minDistIndex = 0;
for (int i = 0; i < shpoints.size(); i++) {
Point2d loc = shpoints.get(i);
double d = Math.sqrt((loc.x - p.x) * (loc.x - p.x) + (loc.y - p.y) * (loc.y - p.y));
if (d < dist) {
dist = d;
minDistIndex = i;
}
}
return shpoints.get(minDistIndex);
}
/**
* Calculate circularity of polygon.
*
* <p>Circularity is computed as: \f[ circ=\frac{4*\pi*A}{P^2} \f] where \f$A\f$ is polygon area
* and \f$P\f$ is its perimeter
*
* @param p Polygon vertices
* @return circularity
*/
double getCircularity(final List<? extends Tuple2d> p) {
double area;
double perim;
BasicPolygons b = new BasicPolygons();
area = b.getPolyArea(p);
perim = b.getPolyPerim(p);
return (4 * Math.PI * area) / (perim * perim);
}
/**
* Calculates weighting based on distribution of window points.
*
* <p>Calculates center of mass of window points and then standard deviations of lengths between
* this point and every other point. Cumulated distributions like protrusions give smaller
* values than elongated ones.
* If input polygon <i>p</i> (which is only part of whole cell shape) is defective, i.e its
* edges cross, the weight is calculated using middle vector defined as mean of coordinates.
*
* @param p Polygon vertices
* @return Weight
*/
double getWeighting(final List<Point2d> p) {
double[] len = new double[p.size()];
BasicPolygons bp = new BasicPolygons();
Vector2d middle;
try { // check if input polygon is correct
middle = new Vector2d(bp.polygonCenterOfMass(p));
} catch (IllegalArgumentException e) { // if not get middle point as mean
double mx = 0;
double my = 0;
for (Point2d v : p) {
mx += v.x;
my += v.y;
}
middle = new Vector2d(mx / p.size(), my / p.size());
}
int i = 0;
// get lengths
for (Point2d v : p) {
Vector2d vec = new Vector2d(middle); // vector between px and middle
vec.sub(v);
len[i++] = vec.length();
}
// get mean
double mean = 0;
for (double d : len) {
mean += d;
}
mean /= p.size();
// get std
double std = 0;
for (double d : len) {
std += Math.pow(d - mean, 2.0);
}
std /= p.size();
std = Math.sqrt(std);
// max of len
// double maxLen = QuimPArrayUtils.arrayMax(len);
LOGGER.debug("getWeighting= " + std);
return std * std;
}
/**
* Class holding lower and upper index of window. Supports comparisons.
*
* <p>Two ranges (l;u) and (l1;u1) are equal if any of these conditions is met:
* <ol>
* <li>they overlap
* <li>they are the same
* <li>one is included in second
* </ol>
*
* @author p.baniukiewicz
* @see WindowIndRange#compareTo(Object)
*/
class WindowIndRange implements Comparable<Object> {
public int lower;
public int upper;
public WindowIndRange() {
upper = 0;
lower = 0;
}
/**
* Create pair of indexes that define window.
*
* @param l lower index
* @param u upper index
*/
WindowIndRange(int l, int u) {
setRange(l, u);
}
@Override
public String toString() {
return "{" + lower + "," + upper + "}";
}
public int hashCode() {
int result = 1;
result = 31 * result + lower;
result = 31 * result + upper;
return result;
}
/**
* Compare two WindowIndRange objects.
*
* @param obj object to compare
* @return true only if ranges does not overlap
*/
@Override
public boolean equals(Object obj) {
if (this == obj) {
return true;
}
if (obj == null) {
return false;
}
if (getClass() != obj.getClass()) {
return false;
}
final WindowIndRange other = (WindowIndRange) obj;
if (upper < other.lower) {
return true;
} else if (lower > other.upper) {
return true;
} else {
return false;
}
}
/**
* Compare two WindowIndRange objects.
*
* <p>The following rules of comparison are used:
* <ol>
* <li>If range1 is below range2 they are not equal
* <li>If range1 is above range2 they are not equal
* </ol>
*
* <p>They are equal in all other cases:
* <ol>
* <li>They are sticked
* <li>One includes other
* <li>They overlap
* </ol>
*
* @param obj Object to compare to this
* @return -1,0,1 expressing relations in windows positions
*/
@Override
public int compareTo(Object obj) {
final WindowIndRange other = (WindowIndRange) obj;
if (this == obj) {
return 0;
}
if (upper < other.lower) {
return -1;
} else if (lower > other.upper) {
return 1;
} else {
return 0;
}
}
/**
* Sets upper and lower indexes to the same value.
*
* @param i Value to set for u and l
*/
public void setSame(int i) {
lower = i;
upper = i;
}
/**
* Set pair of indexes that define window assuring that l < u.
*
* @param l lower index, always smaller
* @param u upper index
*/
public void setRange(int l, int u) {
if (l > u) {
this.lower = u;
this.upper = l;
} else {
this.lower = l;
this.upper = u;
}
}
}
}