List of usage examples for org.apache.commons.math3.analysis.interpolation LoessInterpolator smooth
public final double[] smooth(final double[] xval, final double[] yval) throws NonMonotonicSequenceException, DimensionMismatchException, NoDataException, NotFiniteNumberException, NumberIsTooSmallException
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 w w w . j ava 2 s . co m 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.PSFCreator.java
public void run(ImageProcessor ip) { loadConfiguration();//ww w . j a v a 2 s . com 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:gdsc.smlm.ij.plugins.PSFCreator.java
/** * Fit the new PSF image and show a graph of the amplitude/width * //from w w w.j a v a2s . c o m * @param psf * @param loess * @param averageRange * @param fitCom * @return The width of the PSF in the z-centre */ private double fitPSF(ImageStack psf, LoessInterpolator loess, int cz, double averageRange, double[][] fitCom) { IJ.showStatus("Fitting final PSF"); // Update the box radius since this is used in the fitSpot method. boxRadius = psf.getWidth() / 2; int x = boxRadius, y = boxRadius; FitConfiguration fitConfig = config.getFitConfiguration(); final double shift = fitConfig.getCoordinateShiftFactor(); fitConfig.setInitialPeakStdDev0(fitConfig.getInitialPeakStdDev0() * magnification); fitConfig.setInitialPeakStdDev1(fitConfig.getInitialPeakStdDev1() * magnification); // Need to be updated after the widths have been set fitConfig.setCoordinateShiftFactor(shift); fitConfig.setBackgroundFitting(false); //fitConfig.setLog(new IJLogger()); MemoryPeakResults results = fitSpot(psf, psf.getWidth(), psf.getHeight(), x, y); if (results.size() < 5) { Utils.log(" Final PSF: Not enough fit results %d", results.size()); return 0; } // 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; // Set limits for the fit final float maxWidth = (float) (FastMath.max(fitConfig.getInitialPeakStdDev0(), fitConfig.getInitialPeakStdDev1()) * magnification * 4); final float maxSignal = 2; // PSF is normalised to 1 for (PeakResult peak : results.getResults()) { // Remove bad fits where the width/signal is above the expected final float w = FastMath.max(peak.getXSD(), peak.getYSD()); if (peak.getSignal() > maxSignal || w > maxWidth) continue; z[i] = peak.peak; fitCom[0][peak.peak - 1] = xCoord[i] = peak.getXPosition() - x; fitCom[1][peak.peak - 1] = yCoord[i] = peak.getYPosition() - y; sd[i] = w; a[i] = peak.getAmplitude(); i++; } // Truncate z = Arrays.copyOf(z, i); xCoord = Arrays.copyOf(xCoord, i); yCoord = Arrays.copyOf(yCoord, i); sd = Arrays.copyOf(sd, i); a = Arrays.copyOf(a, i); // Extract the average smoothed range from the individual fits int r = (int) Math.ceil(averageRange / 2); int start = 0, stop = z.length - 1; for (int j = 0; j < z.length; j++) { if (z[j] > cz - r) { start = j; break; } } for (int j = z.length; j-- > 0;) { if (z[j] < cz + r) { stop = j; break; } } // 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[] smoothA = new double[smoothX.length]; double[] newZ = new double[smoothX.length]; int smoothCzIndex = 0; for (int j = start, k = 0; j <= stop; j++, k++) { smoothX[k] = xCoord[j]; smoothY[k] = yCoord[j]; smoothSd[k] = sd[j]; smoothA[k] = a[j]; newZ[k] = z[j]; if (newZ[k] == cz) smoothCzIndex = k; } smoothX = loess.smooth(newZ, smoothX); smoothY = loess.smooth(newZ, smoothY); smoothSd = loess.smooth(newZ, smoothSd); smoothA = loess.smooth(newZ, smoothA); // Update the widths and positions using the magnification final double scale = 1.0 / magnification; for (int j = 0; j < xCoord.length; j++) { xCoord[j] *= scale; yCoord[j] *= scale; sd[j] *= scale; } for (int j = 0; j < smoothX.length; j++) { smoothX[j] *= scale; smoothY[j] *= scale; smoothSd[j] *= scale; } showPlots(z, a, newZ, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cz); // Store the data for replotting this.z = z; this.a = a; this.smoothAz = newZ; this.smoothA = smoothA; this.xCoord = xCoord; this.yCoord = yCoord; this.sd = sd; this.newZ = newZ; this.smoothX = smoothX; this.smoothY = smoothY; this.smoothSd = smoothSd; //maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start); return smoothSd[smoothCzIndex]; }
From source file:org.apache.solr.client.solrj.io.eval.LoessEvaluator.java
@Override public Object doWork(Object... objects) throws IOException { Object first = objects[0];/*www. jav a 2 s .com*/ 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; }
From source file:ubic.basecode.math.linearmodels.MeanVarianceEstimator.java
/** * First ensures that x values are strictly increasing and performs a loess fit afterwards. The loess fit are * determined by <code>BANDWIDTH</code> and <code>ROBUSTNESS_ITERS</code>. * /* ww w .j av a 2 s .com*/ * @param xy * @return loessFit or null if there are less than 3 data points */ private DoubleMatrix2D loessFit(DoubleMatrix2D xy) { assert xy != null; DoubleMatrix1D sx = xy.viewColumn(0); DoubleMatrix1D sy = xy.viewColumn(1); Map<Double, Double> map = new TreeMap<>(); for (int i = 0; i < sx.size(); i++) { if (Double.isNaN(sx.get(i)) || Double.isInfinite(sx.get(i)) || Double.isNaN(sy.get(i)) || Double.isInfinite(sy.get(i))) { continue; } map.put(sx.get(i), sy.get(i)); } DoubleMatrix2D xyChecked = new DenseDoubleMatrix2D(map.size(), 2); xyChecked.viewColumn(0).assign(ArrayUtils.toPrimitive(map.keySet().toArray(new Double[0]))); xyChecked.viewColumn(1).assign(ArrayUtils.toPrimitive(map.values().toArray(new Double[0]))); // in R: // loess(c(1:5),c(1:5)^2,f=0.5,iter=3) // Note: we start to loose some precision here in comparison with R's loess FIXME why? does it matter? DoubleMatrix2D loessFit = new DenseDoubleMatrix2D(xyChecked.rows(), xyChecked.columns()); // try { // fit a loess curve LoessInterpolator loessInterpolator = new LoessInterpolator(MeanVarianceEstimator.BANDWIDTH, MeanVarianceEstimator.ROBUSTNESS_ITERS); double[] loessY = loessInterpolator.smooth(xyChecked.viewColumn(0).toArray(), xyChecked.viewColumn(1).toArray()); loessFit.viewColumn(0).assign(xyChecked.viewColumn(0)); loessFit.viewColumn(1).assign(loessY); return loessFit; }