List of usage examples for org.apache.commons.math3.analysis.interpolation UnivariateInterpolator interpolate
UnivariateFunction interpolate(double xval[], double yval[]);
From source file:oct.analysis.application.comp.EZWorker.java
@Override protected EZEdgeCoord doInBackground() throws Exception { int foveaCenterXPosition = analysisManager.getFoveaCenterXPosition(); /*//from w w w . j a va 2s . 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.meteoinfo.math.interpolate.InterpUtil.java
/** * Make interpolation function//from www. j a v a2s . c om * @param x X data * @param y Y data * @param kind Specifies the kind of interpolation as a string (linear, 'spline'). * @return Interpolation function */ public static UnivariateFunction getInterpFunc(Array x, Array y, String kind) { double[] xd = (double[]) ArrayUtil.copyToNDJavaArray(x); double[] yd = (double[]) ArrayUtil.copyToNDJavaArray(y); UnivariateInterpolator li; switch (kind) { case "spline": case "cubic": li = new SplineInterpolator(); break; case "akima": li = new AkimaSplineInterpolator(); break; case "divided": li = new DividedDifferenceInterpolator(); break; case "loess": li = new LoessInterpolator(); break; case "neville": li = new NevilleInterpolator(); break; default: li = new LinearInterpolator(); break; } UnivariateFunction psf = li.interpolate(xd, yd); return psf; }
From source file:uk.ac.diamond.scisoft.analysis.dataset.function.Interpolation1D.java
public static Dataset interpolate(IDataset oldx, IDataset oldy, IDataset newx, UnivariateInterpolator interpolator) { //TODO more sanity checks on inputs DoubleDataset dx = (DoubleDataset) DatasetUtils.cast(oldx, Dataset.FLOAT64); DoubleDataset dy = (DoubleDataset) DatasetUtils.cast(oldy, Dataset.FLOAT64); boolean sorted = true; double maxtest = dx.getDouble(0); int count = dx.getSize(); for (int i = 1; i < count; i++) { if (maxtest > dx.getDouble(i)) { sorted = false;// w ww . jav a2 s . c o m break; } maxtest = dx.getDouble(i); } double[] sortedx = null; double[] sortedy = null; if (!sorted) { IntegerDataset sIdx = getIndiciesOfSorted(dx); sortedx = new double[dx.getData().length]; sortedy = new double[dy.getData().length]; for (int i = 0; i < sIdx.getSize(); i++) { sortedx[i] = dx.getDouble(sIdx.get(i)); sortedy[i] = dy.getDouble(sIdx.get(i)); } } else { sortedx = dx.getData(); sortedy = dy.getData(); } UnivariateFunction func = interpolator.interpolate(sortedx, sortedy); Dataset newy = DatasetFactory.zeros(newx.getShape(), Dataset.FLOAT64); newy.setName(oldy.getName() + "_interpolated"); count = newy.getSize(); double val = 0; for (int i = 0; i < count; i++) { try { val = func.value(newx.getDouble(i)); } catch (OutOfRangeException e) { val = 0; } newy.set(val, i); } return newy; }
From source file:uk.ac.diamond.scisoft.ncd.calibration.NCDAbsoluteCalibration.java
public void setAbsoluteData(List<Amount<ScatteringVector>> lstAbsQ, Dataset absI, Unit<ScatteringVector> unit) { absQ = new DoubleDataset(lstAbsQ.size()); for (int idx = 0; idx < lstAbsQ.size(); idx++) { Amount<ScatteringVector> vec = lstAbsQ.get(idx); absQ.set(vec.doubleValue(unit), idx); }/*from www .ja va 2 s . c o m*/ this.absI = absI.clone(); UnivariateInterpolator interpolator = new SplineInterpolator(); absInterpolate = interpolator.interpolate((double[]) absQ.getBuffer(), (double[]) absI.getBuffer()); }
From source file:uk.ac.diamond.scisoft.ncd.core.DegreeOfOrientation.java
public Object[] process(Serializable buffer, Serializable axis, final int[] dimensions) { double[] parentaxis = (double[]) ConvertUtils.convert(axis, double[].class); float[] parentdata = (float[]) ConvertUtils.convert(buffer, float[].class); int size = dimensions[dimensions.length - 1]; double[] myaxis = new double[size]; double[] mydata = new double[size]; double[] cos2data = new double[size]; double[] sin2data = new double[size]; double[] sincosdata = new double[size]; for (int i = 0; i < parentaxis.length; i++) { myaxis[i] = Math.toRadians(parentaxis[i]); mydata[i] = parentdata[i];//from ww w. j av a 2 s . c o m float cos2alpha = (float) Math.cos(2.0 * myaxis[i]); float sin2alpha = (float) Math.sin(2.0 * myaxis[i]); cos2data[i] = (1.0f + cos2alpha) * parentdata[i] / 2.0; sin2data[i] = (1.0f - cos2alpha) * parentdata[i] / 2.0; sincosdata[i] = sin2alpha * parentdata[i] / 2.0; } UnivariateInterpolator interpolator = new SplineInterpolator(); UnivariateFunction function = interpolator.interpolate(myaxis, mydata); UnivariateFunction cos2Function = interpolator.interpolate(myaxis, cos2data); UnivariateFunction sin2Function = interpolator.interpolate(myaxis, sin2data); UnivariateFunction sincosFunction = interpolator.interpolate(myaxis, sincosdata); UnivariateIntegrator integrator = new IterativeLegendreGaussIntegrator(15, BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY, BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY); try { float cos2mean = (float) integrator.integrate(INTEGRATION_POINTS, cos2Function, myaxis[0], myaxis[myaxis.length - 1]); float sin2mean = (float) integrator.integrate(INTEGRATION_POINTS, sin2Function, myaxis[0], myaxis[myaxis.length - 1]); float sincosmean = (float) integrator.integrate(INTEGRATION_POINTS, sincosFunction, myaxis[0], myaxis[myaxis.length - 1]); float norm = (float) integrator.integrate(INTEGRATION_POINTS, function, myaxis[0], myaxis[myaxis.length - 1]); cos2mean /= norm; sin2mean /= norm; sincosmean /= norm; float result = (float) Math.sqrt(Math.pow(cos2mean - sin2mean, 2) - 4.0 * sincosmean * sincosmean); double angle = MathUtils.normalizeAngle(Math.atan2(2.0 * sincosmean, cos2mean - sin2mean) / 2.0, Math.PI); Object[] output = new Object[] { new float[] { result }, new float[] { (float) Math.toDegrees(angle) }, new float[] { (float) (result * Math.cos(angle)), (float) (result * Math.sin(angle)) }, }; return output; } catch (TooManyEvaluationsException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } catch (MaxCountExceededException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } }
From source file:uk.ac.diamond.scisoft.ncd.core.SaxsInvariant.java
public Object[] process(Serializable buffer, Serializable errors, Serializable axis, final int[] dimensions) { double[] parentaxis = (double[]) ConvertUtils.convert(axis, double[].class); float[] parentdata = (float[]) ConvertUtils.convert(buffer, float[].class); double[] parenterrors = (double[]) ConvertUtils.convert(errors, double[].class); int shift = (parentaxis[0] > 0 ? 1 : 0); int size = dimensions[dimensions.length - 1] + shift; double[] myaxis = new double[size]; double[] mydata = new double[size]; double[] myerrors = new double[size]; if (shift > 0) { myaxis[0] = 0.0;/*from w w w . ja v a2 s . c o m*/ mydata[0] = 0.0; myerrors[0] = 0.0; } for (int i = 0; i < parentaxis.length; i++) { myaxis[i + shift] = parentaxis[i]; mydata[i + shift] = parentdata[i] * parentaxis[i] * parentaxis[i]; myerrors[i + shift] = parenterrors[i] * Math.pow(parentaxis[i], 4); } UnivariateInterpolator interpolator = new SplineInterpolator(); UnivariateFunction function = interpolator.interpolate(myaxis, mydata); UnivariateIntegrator integrator = new IterativeLegendreGaussIntegrator(15, BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY, BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY); try { float result = (float) integrator.integrate(INTEGRATION_POINTS, function, 0.0, myaxis[myaxis.length - 1]); IDataset data = new FloatDataset(parentdata, dimensions); IDataset qaxis = new DoubleDataset(parentaxis, dimensions); PorodPlotData porodPlotData = (PorodPlotData) SaxsAnalysisPlotType.POROD_PLOT.getSaxsPlotDataObject(); SimpleRegression regression = porodPlotData.getPorodPlotParameters(data.squeeze(), qaxis.squeeze()); Amount<Dimensionless> c4 = porodPlotData.getC4(regression); result += (float) (c4.getEstimatedValue() / myaxis[myaxis.length - 1]); double error = 0.0; for (int i = 0; i < myaxis.length; i++) { int idx1 = Math.max(0, i - 1); int idx2 = Math.min(myaxis.length - 1, i + 1); error += Math.pow((myaxis[idx2] - myaxis[idx1]), 2) * myerrors[i] / 4.0; } error += Math.pow(c4.getAbsoluteError() / myaxis[myaxis.length - 1], 2); return new Object[] { new float[] { result }, new double[] { error } }; } catch (TooManyEvaluationsException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } catch (MaxCountExceededException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } }