List of usage examples for org.apache.commons.math3.analysis.interpolation LoessInterpolator LoessInterpolator
public LoessInterpolator(double bandwidth, int robustnessIters)
From source file:gdsc.smlm.ij.plugins.PSFCreator.java
public void run(ImageProcessor ip) { loadConfiguration();/* ww w.j a v a 2s. c o m*/ BasePoint[] spots = getSpots(); if (spots.length == 0) { IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px"); return; } ImageStack stack = getImageStack(); final int width = imp.getWidth(); final int height = imp.getHeight(); final int currentSlice = imp.getSlice(); // Adjust settings for a single maxima config.setIncludeNeighbours(false); fitConfig.setDuplicateDistance(0); ArrayList<double[]> centres = new ArrayList<double[]>(spots.length); int iterations = 1; LoessInterpolator loess = new LoessInterpolator(smoothing, iterations); // TODO - The fitting routine may not produce many points. In this instance the LOESS interpolator // fails to smooth the data very well. A higher bandwidth helps this but perhaps // try a different smoothing method. // For each spot Utils.log(TITLE + ": " + imp.getTitle()); Utils.log("Finding spot locations..."); Utils.log(" %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2)); StoredDataStatistics averageSd = new StoredDataStatistics(); StoredDataStatistics averageA = new StoredDataStatistics(); Statistics averageRange = new Statistics(); for (int n = 1; n <= spots.length; n++) { BasePoint spot = spots[n - 1]; final int x = (int) spot.getX(); final int y = (int) spot.getY(); MemoryPeakResults results = fitSpot(stack, width, height, x, y); if (results.size() < 5) { Utils.log(" Spot %d: Not enough fit results %d", n, results.size()); continue; } // Get the results for the spot centre and width double[] z = new double[results.size()]; double[] xCoord = new double[z.length]; double[] yCoord = new double[z.length]; double[] sd = new double[z.length]; double[] a = new double[z.length]; int i = 0; for (PeakResult peak : results.getResults()) { z[i] = peak.peak; xCoord[i] = peak.getXPosition() - x; yCoord[i] = peak.getYPosition() - y; sd[i] = FastMath.max(peak.getXSD(), peak.getYSD()); a[i] = peak.getAmplitude(); i++; } // Smooth the amplitude plot double[] smoothA = loess.smooth(z, a); // Find the maximum amplitude int maximumIndex = findMaximumIndex(smoothA); // Find the range at a fraction of the max. This is smoothed to find the X/Y centre int start = 0, stop = smoothA.length - 1; double limit = smoothA[maximumIndex] * amplitudeFraction; for (int j = 0; j < smoothA.length; j++) { if (smoothA[j] > limit) { start = j; break; } } for (int j = smoothA.length; j-- > 0;) { if (smoothA[j] > limit) { stop = j; break; } } averageRange.add(stop - start + 1); // Extract xy centre coords and smooth double[] smoothX = new double[stop - start + 1]; double[] smoothY = new double[smoothX.length]; double[] smoothSd = new double[smoothX.length]; double[] newZ = new double[smoothX.length]; for (int j = start, k = 0; j <= stop; j++, k++) { smoothX[k] = xCoord[j]; smoothY[k] = yCoord[j]; smoothSd[k] = sd[j]; newZ[k] = z[j]; } smoothX = loess.smooth(newZ, smoothX); smoothY = loess.smooth(newZ, smoothY); smoothSd = loess.smooth(newZ, smoothSd); // Since the amplitude is not very consistent move from this peak to the // lowest width which is the in-focus spot. maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start); // Find the centre at the amplitude peak double cx = smoothX[maximumIndex] + x; double cy = smoothY[maximumIndex] + y; int cz = (int) newZ[maximumIndex]; double csd = smoothSd[maximumIndex]; double ca = smoothA[maximumIndex + start]; // The average should weight the SD using the signal for each spot averageSd.add(smoothSd[maximumIndex]); averageA.add(ca); if (ignoreSpot(n, z, a, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cx, cy, cz, csd)) { Utils.log(" Spot %d was ignored", n); continue; } // Store result Utils.log(" Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f\n", n, cx, cy, cz, csd, ca); centres.add(new double[] { cx, cy, cz, csd }); } if (interactiveMode) { imp.setSlice(currentSlice); imp.setOverlay(null); } if (centres.isEmpty()) { String msg = "No suitable spots could be identified centres"; Utils.log(msg); IJ.error(TITLE, msg); return; } // Find the limits of the z-centre int minz = (int) centres.get(0)[2]; int maxz = minz; for (double[] centre : centres) { if (minz > centre[2]) minz = (int) centre[2]; else if (maxz < centre[2]) maxz = (int) centre[2]; } IJ.showStatus("Creating PSF image"); // Create a stack that can hold all the data. ImageStack psf = createStack(stack, minz, maxz, magnification); // For each spot Statistics stats = new Statistics(); boolean ok = true; for (int i = 0; ok && i < centres.size(); i++) { double progress = (double) i / centres.size(); final double increment = 1.0 / (stack.getSize() * centres.size()); IJ.showProgress(progress); double[] centre = centres.get(i); // Extract the spot float[][] spot = new float[stack.getSize()][]; Rectangle regionBounds = null; for (int slice = 1; slice <= stack.getSize(); slice++) { ImageExtractor ie = new ImageExtractor((float[]) stack.getPixels(slice), width, height); if (regionBounds == null) regionBounds = ie.getBoxRegionBounds((int) centre[0], (int) centre[1], boxRadius); spot[slice - 1] = ie.crop(regionBounds); } float b = getBackground(spot); stats.add(b); subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height); // This takes a long time so this should track progress ok = addToPSF(maxz, magnification, psf, centre, spot, regionBounds, progress, increment); } IJ.showProgress(1); if (threadPool != null) { threadPool.shutdownNow(); threadPool = null; } if (!ok) return; final double avSd = getAverage(averageSd, averageA, 2); Utils.log(" Average background = %.2f, Av. SD = %s px", stats.getMean(), Utils.rounded(avSd, 4)); normalise(psf, maxz, avSd * magnification); IJ.showProgress(1); psfImp = Utils.display("PSF", psf); psfImp.setSlice(maxz); psfImp.resetDisplayRange(); psfImp.updateAndDraw(); double[][] fitCom = new double[2][psf.getSize()]; Arrays.fill(fitCom[0], Double.NaN); Arrays.fill(fitCom[1], Double.NaN); double fittedSd = fitPSF(psf, loess, maxz, averageRange.getMean(), fitCom); // Compute the drift in the PSF: // - Use fitted centre if available; otherwise find CoM for each frame // - express relative to the average centre double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / magnification); double[] slice = Utils.newArray(psf.getSize(), 1, 1.0); String title = TITLE + " CoM Drift"; Plot2 plot = new Plot2(title, "Slice", "Drift (nm)"); double[] limitsX = getLimits(com[0]); double[] limitsY = getLimits(com[1]); plot.setLimits(1, psf.getSize(), Math.min(limitsX[0], limitsY[0]), Math.max(limitsX[1], limitsY[1])); plot.setColor(Color.red); plot.addPoints(slice, com[0], Plot.DOT); plot.addPoints(slice, loess.smooth(slice, com[0]), Plot.LINE); plot.setColor(Color.blue); plot.addPoints(slice, com[1], Plot.DOT); plot.addPoints(slice, loess.smooth(slice, com[1]), Plot.LINE); Utils.display(title, plot); // TODO - Redraw the PSF with drift correction applied. // This means that the final image should have no drift. // This is relevant when combining PSF images. It doesn't matter too much for simulations // unless the drift is large. // Add Image properties containing the PSF details final double fwhm = getFWHM(psf, maxz); psfImp.setProperty("Info", XmlUtils .toXML(new PSFSettings(maxz, nmPerPixel / magnification, nmPerSlice, centres.size(), fwhm))); Utils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, Utils.rounded(nmPerPixel / magnification, 3), Utils.rounded(nmPerSlice, 3), centres.size(), Utils.rounded(fittedSd * nmPerPixel, 4), Utils.rounded(fwhm)); createInteractivePlots(psf, maxz, nmPerPixel / magnification, fittedSd * nmPerPixel); IJ.showStatus(""); }
From source file:com.github.lindenb.jvarkit.tools.redon.CopyNumber01.java
private UnivariateInterpolator createInterpolator() { UnivariateInterpolator interpolator = null; interpolator = new LoessInterpolator(0.5, 4); //interpolator = new NevilleInterpolator(); return interpolator; }
From source file:MSUmpire.LCMSPeakStructure.LCMSPeakMS1.java
public void GenerateMassCalibrationRTMap() throws IOException { String pngfile = FilenameUtils.getFullPath(ScanCollectionName) + "/" + FilenameUtils.getBaseName(ScanCollectionName) + "_masscaliRT.png"; XYSeries series = new XYSeries("PSM"); XYSeriesCollection xySeriesCollection = new XYSeriesCollection(); LoessInterpolator loessInterpolator = new LoessInterpolator(0.75, //bandwidth, 2//robustnessIters );//from ww w . j a va 2 s .c o m for (PSM psm : this.IDsummary.PSMList.values()) { float ppm = InstrumentParameter.CalcSignedPPM(psm.ObserPrecursorMass, psm.NeutralPepMass); series.add(new XYDataItem(psm.RetentionTime, ppm)); } double x[] = new double[IDsummary.PSMList.size()]; double y[] = new double[x.length]; double currentmin = 0f; for (int i = 0; i < series.getItemCount(); i++) { x[i] = (double) series.getX(i); if (x[i] <= currentmin) { x[i] = currentmin + 0.0001f; } currentmin = x[i]; y[i] = (double) series.getY(i); } Masscalibrationfunction = loessInterpolator.interpolate(x, y); XYSeries smoothline = new XYSeries("Loess Regression"); double xvalue = x[0]; while (xvalue < x[x.length - 1]) { smoothline.add(xvalue, Masscalibrationfunction.value(xvalue)); xvalue += 0.05d; } xySeriesCollection.addSeries(smoothline); xySeriesCollection.addSeries(series); JFreeChart chart = ChartFactory.createScatterPlot("Mass calibration", "RT", "Mass error (ppm)", xySeriesCollection, PlotOrientation.VERTICAL, true, true, false); XYPlot xyPlot = (XYPlot) chart.getPlot(); xyPlot.setDomainCrosshairVisible(true); xyPlot.setRangeCrosshairVisible(true); XYItemRenderer renderer = xyPlot.getRenderer(); renderer.setSeriesPaint(1, Color.blue); renderer.setSeriesPaint(0, Color.BLACK); renderer.setSeriesShape(1, new Ellipse2D.Double(0, 0, 3, 3)); renderer.setSeriesStroke(1, new BasicStroke(1.0f)); xyPlot.setBackgroundPaint(Color.white); ChartUtilities.saveChartAsPNG(new File(pngfile), chart, 1000, 600); }
From source file:com.github.brandtg.stl.StlDecomposition.java
/** * Performs weighted Loess smoothing on a series. * * <p>/* w w w . j a va2 s .c om*/ * Does not assume contiguous time. * </p> * * @param times * The times. * @param series * The time series data. * @param bandwidth * The amount of neighbor points to consider for each point in Loess. * @param weights * The weights to use for smoothing, if null, equal weights are assumed. * @return * Loess-smoothed series. */ private double[] loessSmooth(double[] times, double[] series, double bandwidth, double[] weights) { if (weights == null) { return new LoessInterpolator(bandwidth, config.getLoessRobustnessIterations()).smooth(times, series); } else { return new LoessInterpolator(bandwidth, config.getLoessRobustnessIterations()).smooth(times, series, weights); } }
From source file:com.github.brandtg.stl.StlDecomposition.java
private TimesAndValues padEdges(double[] ts, double[] xs) { // Find step between times double step = Math.abs(ts[1] - ts[0]); // Times (assuming uniform double[] paddedTimes = new double[ts.length + 2]; System.arraycopy(ts, 0, paddedTimes, 1, ts.length); paddedTimes[0] = paddedTimes[1] - step; paddedTimes[paddedTimes.length - 1] = paddedTimes[paddedTimes.length - 2] + step; // Series// w ww .j ava 2s .com double[] paddedSeries = new double[xs.length + 2]; System.arraycopy(xs, 0, paddedSeries, 1, xs.length); // Use Loess at ends to pad // n.b. For some reason, this can result in NaN values - perhaps similar to // https://issues.apache.org/jira/browse/MATH-296. If we see NaN, just "extrapolate" by copying // the end points :( double left = paddedSeries[1]; double right = paddedSeries[paddedSeries.length - 2]; double bandwidth = 0.3; if (ts.length * bandwidth > 2) { PolynomialSplineFunction loess = new LoessInterpolator(bandwidth, 2).interpolate(ts, xs); double loessLeft = loess.value(ts[0]); if (!Double.isNaN(loessLeft)) { left = loessLeft; } double loessRight = loess.value(ts[ts.length - 1]); if (!Double.isNaN(loessRight)) { right = loessRight; } } paddedSeries[0] = left; paddedSeries[paddedSeries.length - 1] = right; return new TimesAndValues(paddedTimes, paddedSeries); }
From source file:gdsc.smlm.ij.plugins.DriftCalculator.java
private static boolean smooth(double[] newDx, double[] newDy, double[] originalDriftTimePoints, double smoothing, int iterations) { double[][] values = extractValues(originalDriftTimePoints, 0, newDx.length - 1, newDx, newDy); // Smooth/*from ww w .j ava2 s. c om*/ LoessInterpolator loess = new LoessInterpolator(smoothing, iterations); values[1] = loess.smooth(values[0], values[1]); values[2] = loess.smooth(values[0], values[2]); // Add back int n = 0; for (int t = 0; t < newDx.length; t++) { if (originalDriftTimePoints[t] != 0) { newDx[t] = values[1][n]; newDy[t] = values[2][n]; n++; if (Double.isNaN(newDx[t])) { Utils.log("ERROR : Loess smoothing created bad X-estimate at point %d/%d", t, newDx.length); return false; } if (Double.isNaN(newDy[t])) { Utils.log("ERROR : Loess smoothing created bad Y-estimate at point %d/%d", t, newDx.length); return false; } } } return true; }
From source file:gdsc.smlm.ij.plugins.SpotAnalysis.java
private double[] interpolate(double[] xValues2, double[] yValues) { // Smooth the values not in the current on-frames double[] newX = Arrays.copyOf(xValues, xValues.length); double[] newY = Arrays.copyOf(yValues, yValues.length); for (Spot s : onFrames) { newX[s.frame - 1] = -1;//from ww w. j av a 2 s. c om } int c = 0; for (int i = 0; i < newX.length; i++) { if (newX[i] == -1) continue; newX[c] = newX[i]; newY[c] = newY[i]; c++; } newX = Arrays.copyOf(newX, c); newY = Arrays.copyOf(newY, c); double smoothing = 0.25; try { smoothing = Double.parseDouble(smoothingTextField.getText()); if (smoothing < 0.01 || smoothing > 0.9) smoothing = 0.25; } catch (NumberFormatException e) { } LoessInterpolator loess = new LoessInterpolator(smoothing, 1); PolynomialSplineFunction f = loess.interpolate(newX, newY); // Interpolate double[] plotSmooth = new double[xValues.length]; for (int i = 0; i < xValues.length; i++) { // Cannot interpolate outside the bounds of the input data if (xValues[i] < newX[0]) plotSmooth[i] = newY[0]; else if (xValues[i] > newX[newX.length - 1]) plotSmooth[i] = newY[newX.length - 1]; else plotSmooth[i] = f.value(xValues[i]); } return plotSmooth; }
From source file:gdsc.smlm.ij.plugins.BenchmarkFilterAnalysis.java
/** * @param filter/*w w w.ja v a 2s .co m*/ */ private void depthAnalysis(Filter filter) { // TODO : This analysis ignores the partial match distance. // Use the score for each result to get a weighted histogram. if (!depthRecallAnalysis || simulationParameters.fixedDepth) return; // Build a histogram of the number of spots at different depths final double[] depths = depthStats.getValues(); final double range = simulationParameters.depth / simulationParameters.a / 2; double[] limits = { -range, range }; final int bins = Math.max(10, simulationParameters.molecules / 50); double[][] h = Utils.calcHistogram(depths, limits[0], limits[1], bins); double[][] h2 = Utils.calcHistogram(depthFitStats.getValues(), limits[0], limits[1], bins); // To get the number of TP at each depth will require that the filter is run // manually to get the results that pass. MemoryPeakResults results = filter.filter(resultsList.get(0), failCount); double[] depths2 = new double[results.size()]; int count = 0; for (PeakResult r : results.getResults()) { if (r.origValue != 0) depths2[count++] = ((DepthPeakResult) r).depth; } depths2 = Arrays.copyOf(depths2, count); // Build a histogram using the same limits double[][] h3 = Utils.calcHistogram(depths2, limits[0], limits[1], bins); // Convert pixel depth to nm for (int i = 0; i < h[0].length; i++) h[0][i] *= simulationParameters.a; limits[0] *= simulationParameters.a; limits[1] *= simulationParameters.a; // Produce a histogram of the number of spots at each depth String title = TITLE + " Depth Histogram"; Plot2 plot = new Plot2(title, "Depth (nm)", "Frequency"); plot.setLimits(limits[0], limits[1], 0, Maths.max(h[1])); plot.setColor(Color.black); plot.addPoints(h[0], h[1], Plot2.BAR); plot.setColor(Color.blue); plot.addPoints(h[0], h2[1], Plot2.BAR); plot.setColor(Color.red); plot.addPoints(h[0], h3[1], Plot2.BAR); plot.addLabel(0, 0, "Black = Spots; Blue = Fitted; Red = Filtered"); PlotWindow pw = Utils.display(title, plot); // Interpolate final double halfBinWidth = (h[0][1] - h[0][0]) * 0.5; // Remove final value of the histogram as this is at the upper limit of the range (i.e. count zero) h[0] = Arrays.copyOf(h[0], h[0].length - 1); h[1] = Arrays.copyOf(h[1], h[0].length); h2[1] = Arrays.copyOf(h2[1], h[0].length); h3[1] = Arrays.copyOf(h3[1], h[0].length); // Use minimum of 3 points for smoothing // Ensure we use at least 15% of data double bandwidth = Math.max(3.0 / h[0].length, 0.15); LoessInterpolator l = new LoessInterpolator(bandwidth, 1); PolynomialSplineFunction spline = l.interpolate(h[0], h[1]); PolynomialSplineFunction spline2 = l.interpolate(h[0], h2[1]); PolynomialSplineFunction spline3 = l.interpolate(h[0], h3[1]); // Increase the number of points to show a smooth curve double[] points = new double[bins * 5]; limits = Maths.limits(h[0]); final double interval = (limits[1] - limits[0]) / (points.length - 1); double[] v = new double[points.length]; double[] v2 = new double[points.length]; double[] v3 = new double[points.length]; for (int i = 0; i < points.length - 1; i++) { points[i] = limits[0] + i * interval; v[i] = spline.value(points[i]); v2[i] = spline2.value(points[i]); v3[i] = spline3.value(points[i]); points[i] += halfBinWidth; } // Final point on the limit of the spline range int ii = points.length - 1; v[ii] = spline.value(limits[1]); v2[ii] = spline2.value(limits[1]); v3[ii] = spline3.value(limits[1]); points[ii] = limits[1] + halfBinWidth; // Calculate recall for (int i = 0; i < v.length; i++) { v2[i] = v2[i] / v[i]; v3[i] = v3[i] / v[i]; } title = TITLE + " Depth Histogram (normalised)"; plot = new Plot2(title, "Depth (nm)", "Recall"); plot.setLimits(limits[0] + halfBinWidth, limits[1] + halfBinWidth, 0, Maths.max(v2)); plot.setColor(Color.blue); plot.addPoints(points, v2, Plot2.LINE); plot.setColor(Color.red); plot.addPoints(points, v3, Plot2.LINE); plot.addLabel(0, 0, "Blue = Fitted; Red = Filtered"); PlotWindow pw2 = Utils.display(title, plot); if (Utils.isNewWindow()) { Point p = pw.getLocation(); p.y += pw.getHeight(); pw2.setLocation(p); } }
From source file:oct.analysis.application.comp.EZWorker.java
@Override protected EZEdgeCoord doInBackground() throws Exception { int foveaCenterXPosition = analysisManager.getFoveaCenterXPosition(); /*//from w w w . jav a 2 s .c om first get a sharpened version of the OCT and use that to obtain the segmentation of the Bruch's membrane. Use a Loess interpolation algorithm to smooth out imperfetions in the segmentation line. */ UnivariateInterpolator interpolator = new LoessInterpolator(0.1, 0); ArrayList<Point> rawBrmPoints = new ArrayList<>(analysisManager .getSegmentation(new SharpenOperation(15, 0.5F)).getSegment(Segmentation.BrM_SEGMENT)); double[][] brmSeg = Util.getXYArraysFromPoints(rawBrmPoints); UnivariateFunction brmInterp = interpolator.interpolate(brmSeg[0], brmSeg[1]); BufferedImage sharpOCT = analysisManager.getSharpenedOctImage(8.5D, 1.0F); setProgress(10); /* Starting from the identified location of the fovea search northward in the image until the most northern pixels northward (in a 3x3 matrix of pixels arround the the search point (X,Y) ) are black (ie. the search matrix is has found that the search point isn't totally surrounded by white pixels). Then a recursive search algorithm determines if the black area signifies the seperation between bands or simply represents a closed (a black blob entirely surrounded by white pixels) black band. It will continue searching northward in the image until it can find an open region of all blak pixels. Once this is found it will find the contour of the edge between the black and white pixels along the width of the image. */ int searchY = (int) Math.round(brmInterp.value(foveaCenterXPosition)) + 1; do { searchY--; } while (Util.calculateGrayScaleValue(sharpOCT.getRGB(foveaCenterXPosition, searchY)) > 0 || !isContrastPoint(foveaCenterXPosition, searchY, sharpOCT)); LinkedList<Point> contour = new LinkedList<>(); Point startPoint = new Point(foveaCenterXPosition, searchY); //find contour by searching for white pixel boundary to te right of the fovea contour.add(findContourRight(startPoint, Cardinality.SOUTH, startPoint, Cardinality.SOUTH, contour, sharpOCT, 0)); //search until open black area found (ie. if the search algorithm arrives back at //the starting pixel keep moving north to next black area to search) while (contour.get(0).equals(startPoint)) { contour = new LinkedList<>(); do { searchY--; } while (Util.calculateGrayScaleValue(sharpOCT.getRGB(foveaCenterXPosition, searchY)) == 0); do { searchY--; } while (Util.calculateGrayScaleValue(sharpOCT.getRGB(foveaCenterXPosition, searchY)) > 0 || isSurroundedByWhite(foveaCenterXPosition, searchY, sharpOCT)); startPoint = new Point(foveaCenterXPosition, searchY); contour.add(findContourRight(startPoint, Cardinality.SOUTH, startPoint, Cardinality.SOUTH, contour, sharpOCT, 0)); } setProgress(20); //open balck space found, complete contour to left of fovea contour.add( findContourLeft(startPoint, Cardinality.SOUTH, startPoint, Cardinality.SOUTH, contour, sharpOCT)); analysisManager.getImgPanel().setDrawPoint(new Point(foveaCenterXPosition, searchY)); setProgress(30); /* since the contour can snake around due to aberations and low image density we need to create a single line (represented by points) from left to right to represent the countour. This is easily done by building a line of points consisting of the point with the largest Y value (furthest from the top of the image) at each X value. This eliminates overhangs from the contour line. */ Map<Double, List<Point>> grouped = contour.stream().collect(Collectors.groupingBy(Point::getX)); List<Point> refinedEZContour = grouped.values().stream().map((List<Point> points) -> { int maxY = points.stream().mapToInt((Point p) -> p.y).min().getAsInt(); return new Point(points.get(0).x, maxY); }).sorted((Point p1, Point p2) -> Integer.compare(p1.x, p2.x)).collect(Collectors.toList()); setProgress(35); /* Starting from the identified location of the fovea search southward in the image until the most southern pixels (in a 3x3 matrix of pixels arround the the search point (X,Y) ) are black (ie. the search matrix has found that the search point isn't totally surrounded by white pixels). Then a recursive search algorithm determines if the black area signifies the bottom of the Bruch's membrane or simply represents a closed (a black blob entirely surrounded by white pixels) black band. It will continue searching southward in the image until it can find an open region of all black pixels. Once this is found it will find the contour of the edge between the black and white pixels, along the width of the image, of the bottom of the Bruch's membrane. */ // sharpOCT = getSharpenedOctImage(5D, 1.0F); searchY = (int) Math.round(brmInterp.value(foveaCenterXPosition)); do { searchY++; } while (Util.calculateGrayScaleValue(sharpOCT.getRGB(foveaCenterXPosition, searchY)) > 0 || isSurroundedByWhite(foveaCenterXPosition, searchY, sharpOCT)); contour = new LinkedList<>(); startPoint = new Point(foveaCenterXPosition, searchY); /* Find contour by searching for white pixel boundary to te right of the fovea. Sometimes the crap below the Bruchs membrane causes too much interferance for the algorithm to work properly so we must tweak some of the parameters of the sharpening performed on the image until the algorithm succedes or we can no longer tweak parameters. In the case of the later event we can use the raw segmented Bruchs membrane as a substitute to keep the method from failing. */ contour.add(findContourRight(startPoint, Cardinality.NORTH, startPoint, Cardinality.NORTH, contour, sharpOCT, 0)); double filtValue = 8.5D; boolean tweakFailed = false; while (contour.contains(null)) { contour = new LinkedList<>(); filtValue -= 0.5D; System.out.println("Reducing sigma to " + filtValue); if (filtValue <= 0D) { tweakFailed = true; break; } sharpOCT = analysisManager.getSharpenedOctImage(8.5D, 1.0F); contour.add(findContourRight(startPoint, Cardinality.NORTH, startPoint, Cardinality.NORTH, contour, sharpOCT, 0)); } if (tweakFailed) { contour = new LinkedList<>(rawBrmPoints); } else { //search until open black area found (ie. if the search algorithm arrives back at //the starting pixel keep moving south to next black area to search) while (contour.get(0).equals(startPoint)) { contour = new LinkedList<>(); do { searchY++; } while (Util.calculateGrayScaleValue(sharpOCT.getRGB(foveaCenterXPosition, searchY)) == 0); do { searchY++; } while (Util.calculateGrayScaleValue(sharpOCT.getRGB(foveaCenterXPosition, searchY)) > 0 || isSurroundedByWhite(foveaCenterXPosition, searchY, sharpOCT)); startPoint = new Point(foveaCenterXPosition, searchY); contour.add(findContourRight(startPoint, Cardinality.NORTH, startPoint, Cardinality.NORTH, contour, sharpOCT, 0)); } setProgress(45); //open balck space found, complete contour to left of fovea contour.add(findContourLeft(startPoint, Cardinality.NORTH, startPoint, Cardinality.NORTH, contour, sharpOCT)); } setProgress(55); /* since the contour can snake around due to aberations and low image density we need to create a single line (represented by points) from left to right to represent the countour. This is easily done by building a line of points consisting of the point with the smallest Y value (closest to the top of the image) at each X value. This eliminates overhangs from the contour line. */ grouped = contour.stream().collect(Collectors.groupingBy(Point::getX)); List<Point> refinedBruchsMembraneContour = grouped.values().stream().map((List<Point> points) -> { int minY = points.stream().mapToInt((Point p) -> p.y).min().getAsInt(); return new Point(points.get(0).x, minY); }).sorted((Point p1, Point p2) -> Integer.compare(p1.x, p2.x)).collect(Collectors.toList()); setProgress(70); /* use a Loess interpolator again to smooth the new contours of the EZ and Bruch's Membrane */ double[][] refinedContourPoints = Util.getXYArraysFromPoints(refinedEZContour); UnivariateFunction interpEZContour = interpolator.interpolate(refinedContourPoints[0], refinedContourPoints[1]); refinedContourPoints = Util.getXYArraysFromPoints(refinedBruchsMembraneContour); UnivariateFunction interpBruchsContour = interpolator.interpolate(refinedContourPoints[0], refinedContourPoints[1]); /* find the average difference in the distance in the Y between the 10 pixels at each end of the Bruch's Membrane contour and the contour created along the top of the EZ. */ //since the lines are sorted on X position it is easy to align the lines //based on the tails of each line int minX = refinedEZContour.get(0).x; int maxX; //the interpolator can shorten the range of the X values from the original supplied //so we need to test where the end of the range occurs since it isn't directly accessible for (maxX = refinedEZContour.get(refinedEZContour.size() - 1).x; maxX > minX; maxX--) { try { double tmp = interpEZContour.value(maxX) - interpBruchsContour.value(maxX); //if this break is reached we have found the max value the interpolators will allow break; } catch (OutOfRangeException oe) { //do nothing but let loop continue } } double avgDif = Stream .concat(IntStream.range(minX + 30, minX + 50).boxed(), IntStream.range(maxX - 49, maxX - 28).boxed()) .mapToDouble(x -> interpBruchsContour.value(x) - interpEZContour.value(x)).average().getAsDouble(); int height = sharpOCT.getHeight();//make to use in lambda expression List<LinePoint> ezLine = IntStream.rangeClosed(minX, maxX) .mapToObj(x -> new LinePoint(x, height - interpEZContour.value(x) - avgDif)) .collect(Collectors.toList()); List<LinePoint> bmLine = IntStream.rangeClosed(minX, maxX) .mapToObj(x -> new LinePoint(x, height - interpBruchsContour.value(x))) .collect(Collectors.toList()); List<LinePoint> bmUnfiltLine = refinedBruchsMembraneContour.stream() .map((Point p) -> new LinePoint(p.x, height - p.getY())).collect(Collectors.toList()); Util.graphPoints(ezLine, bmLine, bmUnfiltLine); analysisManager.getImgPanel().setDrawnLines( IntStream.rangeClosed(minX, maxX).mapToObj(x -> new LinePoint(x, interpEZContour.value(x))) .collect(Collectors.toList()), IntStream.rangeClosed(minX, maxX).mapToObj(x -> new LinePoint(x, interpBruchsContour.value(x))) .collect(Collectors.toList())); /* Find the difference between the two contours (Bruch's membrane and the EZ + Bruch's membrane) and use this to determine where the edge of the EZ is */ List<LinePoint> diffLine = findDiffWithAdjustment(interpBruchsContour, 0D, interpEZContour, avgDif, minX, maxX); setProgress(90); // List<LinePoint> peaks = Util.findPeaksAndVallies(diffLine); // Util.graphPoints(diffLine, peaks); /* Find the first zero crossings of the difference line on both sides of the fovea. If a zero crossing can't be found then search for the first crossing of a value of 1, then 2, then 3, etc. until an X coordinate of a crossing is found on each side of the fovea. */ OptionalInt ezLeftEdge; double crossingThreshold = 0.25D; do { double filtThresh = crossingThreshold; System.out.println("Crossing threshold = " + crossingThreshold); ezLeftEdge = diffLine.stream().filter(lp -> lp.getY() <= filtThresh && lp.getX() < foveaCenterXPosition) .mapToInt(LinePoint::getX).max(); crossingThreshold += 0.25D; } while (!ezLeftEdge.isPresent()); OptionalInt ezRightEdge; crossingThreshold = 0.25D; do { double filtThresh = crossingThreshold; System.out.println("Crossing threshold = " + crossingThreshold); ezRightEdge = diffLine.stream() .filter(lp -> lp.getY() <= filtThresh && lp.getX() > foveaCenterXPosition) .mapToInt(LinePoint::getX).min(); crossingThreshold += 0.25D; } while (!ezRightEdge.isPresent()); //return findings return new EZEdgeCoord(ezLeftEdge.getAsInt(), ezRightEdge.getAsInt()); }
From source file:org.apache.solr.client.solrj.io.eval.LoessEvaluator.java
@Override public Object doWork(Object... objects) throws IOException { Object first = objects[0];//from ww w .ja v a2 s . c o m double[] x = null; double[] y = null; if (objects.length == 1) { //Only the y values passed y = ((List) first).stream().mapToDouble(value -> ((BigDecimal) value).doubleValue()).toArray(); x = new double[y.length]; for (int i = 0; i < y.length; i++) { x[i] = i; } } else if (objects.length == 2) { Object second = objects[1]; x = ((List) first).stream().mapToDouble(value -> ((BigDecimal) value).doubleValue()).toArray(); y = ((List) second).stream().mapToDouble(value -> ((BigDecimal) value).doubleValue()).toArray(); } LoessInterpolator interpolator = new LoessInterpolator(bandwidth, robustIterations); double[] smooth = interpolator.smooth(x, y); List list = new ArrayList(); for (double yvalue : smooth) { list.add(yvalue); } return list; }