List of usage examples for org.apache.commons.math3.analysis.interpolation SplineInterpolator interpolate
public PolynomialSplineFunction interpolate(double x[], double y[]) throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException
From source file:de.mprengemann.intellij.plugin.androidicons.util.ImageUtils.java
private static BufferedImage rescaleBorder(BufferedImage image, int targetWidth, int targetHeight) throws IOException { if (targetWidth == 0) { targetWidth = 1;// w w w . j a v a2 s .c om } if (targetHeight == 0) { targetHeight = 1; } if (targetHeight > 1 && targetWidth > 1) { throw new IOException(); } int w = image.getWidth(); int h = image.getHeight(); int[] data = image.getRGB(0, 0, w, h, null, 0, w); int[] newData = new int[targetWidth * targetHeight]; for (int x = 0; x < w; x++) { for (int y = 0; y < h; y++) { newData[y * targetWidth + x] = 0x00; } } List<Integer> startPositions = new ArrayList<Integer>(); List<Integer> endPositions = new ArrayList<Integer>(); boolean inBlock = false; if (targetHeight == 1) { for (int x = 0; x < w; x++) { if ((0xff000000 & data[x]) != 0) { if (!inBlock) { inBlock = true; startPositions.add(x); } } else if (inBlock) { endPositions.add(x - 1); inBlock = false; } } if (inBlock) { endPositions.add(w - 1); } } else { for (int y = 0; y < h; y++) { if ((0xff000000 & data[y]) != 0) { if (!inBlock) { inBlock = true; startPositions.add(y); } } else if (inBlock) { endPositions.add(y - 1); inBlock = false; } } if (inBlock) { endPositions.add(h - 1); } } try { SplineInterpolator interpolator = new SplineInterpolator(); PolynomialSplineFunction function = interpolator.interpolate( new double[] { 0f, 1f, Math.max(w - 1, h - 1) }, new double[] { 0f, 1f, Math.max(targetHeight - 1, targetWidth - 1) }); for (int i = 0; i < startPositions.size(); i++) { int start = startPositions.get(i); int end = endPositions.get(i); for (int j = (int) function.value(start); j <= (int) function.value(end); j++) { newData[j] = 0xff000000; } } BufferedImage img = UIUtil.createImage(targetWidth, targetHeight, BufferedImage.TYPE_INT_ARGB); img.setRGB(0, 0, targetWidth, targetHeight, newData, 0, targetWidth); return img; } catch (Exception e) { Logger.getInstance(ImageUtils.class).error("resizeBorder", e); } return null; }
From source file:de.mprengemann.intellij.plugin.androidicons.images.ImageUtils.java
private static BufferedImage rescaleBorder(BufferedImage image, int targetWidth, int targetHeight) { if (targetWidth == 0) { targetWidth = 1;/*from ww w.j a v a2 s . c om*/ } if (targetHeight == 0) { targetHeight = 1; } if (targetHeight > 1 && targetWidth > 1) { throw new Wrong9PatchException(); } int w = image.getWidth(); int h = image.getHeight(); int[] data = image.getRGB(0, 0, w, h, null, 0, w); int[] newData = new int[targetWidth * targetHeight]; for (int x = 0; x < w; x++) { for (int y = 0; y < h; y++) { newData[y * targetWidth + x] = 0x00; } } List<Integer> startPositions = new ArrayList<Integer>(); List<Integer> endPositions = new ArrayList<Integer>(); boolean inBlock = false; if (targetHeight == 1) { for (int x = 0; x < w; x++) { if ((0xff000000 & data[x]) != 0) { if (!inBlock) { inBlock = true; startPositions.add(x); } } else if (inBlock) { endPositions.add(x - 1); inBlock = false; } } if (inBlock) { endPositions.add(w - 1); } } else { for (int y = 0; y < h; y++) { if ((0xff000000 & data[y]) != 0) { if (!inBlock) { inBlock = true; startPositions.add(y); } } else if (inBlock) { endPositions.add(y - 1); inBlock = false; } } if (inBlock) { endPositions.add(h - 1); } } try { SplineInterpolator interpolator = new SplineInterpolator(); PolynomialSplineFunction function = interpolator.interpolate( new double[] { 0f, 1f, Math.max(w - 1, h - 1) }, new double[] { 0f, 1f, Math.max(targetHeight - 1, targetWidth - 1) }); for (int i = 0; i < startPositions.size(); i++) { int start = startPositions.get(i); int end = endPositions.get(i); for (int j = (int) function.value(start); j <= (int) function.value(end); j++) { newData[j] = 0xff000000; } } BufferedImage img = UIUtil.createImage(targetWidth, targetHeight, BufferedImage.TYPE_INT_ARGB); img.setRGB(0, 0, targetWidth, targetHeight, newData, 0, targetWidth); return img; } catch (Exception e) { Logger.getInstance(ImageUtils.class).error("resizeBorder", e); } return null; }
From source file:gdsc.smlm.model.AiryPSFModel.java
private static synchronized void createAiryDistribution() { if (spline != null) return;/* w w w . j a v a2s. c o m*/ final double relativeAccuracy = 1e-4; final double absoluteAccuracy = 1e-8; final int minimalIterationCount = 3; final int maximalIterationCount = 32; UnivariateIntegrator integrator = new SimpsonIntegrator(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); UnivariateFunction f = new UnivariateFunction() { public double value(double x) { // The pattern profile is in one dimension. // Multiply by the perimeter of a circle to convert to 2D volume then normalise by 4 pi //return AiryPattern.intensity(x) * 2 * Math.PI * x / (4 * Math.PI); return AiryPattern.intensity(x) * 0.5 * x; } }; // Integrate up to a set number of dark rings int samples = 1000; final double step = RINGS[SAMPLE_RINGS] / samples; double to = 0, from = 0; r = new double[samples + 1]; sum = new double[samples + 1]; for (int i = 1; i < sum.length; i++) { from = to; r[i] = to = step * i; sum[i] = integrator.integrate(2000, f, from, to) + sum[i - 1]; } if (DoubleEquality.relativeError(sum[samples], POWER[SAMPLE_RINGS]) > 1e-3) throw new RuntimeException("Failed to create the Airy distribution"); SplineInterpolator si = new SplineInterpolator(); spline = si.interpolate(sum, r); }
From source file:imperial.modaclouds.monitoring.sda.weka.CreateArff.java
/** * Create arff file given the data/*from w ww . j a v a2s . c o m*/ * * @param timestamps_str the timestamps data * @param data the values of the metrics * @param metricName the metric name * @param fileName the file name to keep the arff file */ public static void create(ArrayList<ArrayList<String>> timestamps_str, ArrayList<ArrayList<String>> data, ArrayList<String> metricName, String fileName) { System.out.println("data: " + data.get(0)); long min_timestamp = Long.valueOf(Collections.min(timestamps_str.get(0))); long max_timestamp = Long.valueOf(Collections.max(timestamps_str.get(0))); for (int i = 1; i < timestamps_str.size(); i++) { long min_temp = Long.valueOf(Collections.min(timestamps_str.get(i))); long max_temp = Long.valueOf(Collections.max(timestamps_str.get(i))); if (max_temp < max_timestamp) { max_timestamp = max_temp; } if (min_temp > min_timestamp) { min_timestamp = min_temp; } } for (int i = 0; i < timestamps_str.size(); i++) { Iterator<String> iter_time = timestamps_str.get(i).iterator(); Iterator<String> iter_data = data.get(i).iterator(); while (iter_time.hasNext()) { long temp_timestamps = Long.valueOf(iter_time.next()); if (temp_timestamps < min_timestamp || temp_timestamps > max_timestamp) { iter_time.remove(); iter_data.next(); iter_data.remove(); } } } double[] timestamps = convertDoubles(timestamps_str.get(0)); double[] targetData = convertDoubles(data.get(0)); double[][] otherData = new double[data.size() - 1][timestamps.length]; for (int i = 0; i < data.size() - 1; i++) { double[] timestamps_temp = convertDoubles(timestamps_str.get(i)); double[] targetData_temp = convertDoubles(data.get(i)); SplineInterpolator spline = new SplineInterpolator(); Map<Double, Integer> map = new TreeMap<Double, Integer>(); for (int j = 0; j < timestamps_temp.length; j++) { map.put(timestamps_temp[j], j); } Collection<Integer> indices = map.values(); int[] indices_int = ArrayUtils.toPrimitive(indices.toArray(new Integer[indices.size()])); double[] timestamps_temp_new = new double[indices_int.length]; double[] targetData_temp_new = new double[indices_int.length]; for (int j = 0; j < indices_int.length; j++) { timestamps_temp_new[j] = timestamps_temp[indices_int[j]]; targetData_temp_new[j] = targetData_temp[indices_int[j]]; } PolynomialSplineFunction polynomical = spline.interpolate(timestamps_temp_new, targetData_temp_new); for (int j = 0; j < timestamps.length; j++) { try { otherData[i][j] = polynomical.value(timestamps[j]); } catch (Exception ex) { otherData[i][j] = targetData_temp_new[j]; } } } ArrayList<Attribute> attributes; Instances dataSet; attributes = new ArrayList<Attribute>(); for (String metric : metricName) { attributes.add(new Attribute(metric)); } dataSet = new Instances("data", attributes, 0); for (int i = 0; i < timestamps.length; i++) { double[] instanceValue1 = new double[dataSet.numAttributes()]; instanceValue1[0] = timestamps[i]; instanceValue1[1] = targetData[i]; for (int j = 0; j < data.size() - 1; j++) { instanceValue1[2 + j] = otherData[j][i]; } DenseInstance denseInstance1 = new DenseInstance(1.0, instanceValue1); dataSet.add(denseInstance1); } ArffSaver saver = new ArffSaver(); saver.setInstances(dataSet); try { String workingDir = System.getProperty("user.dir"); System.out.println("workingDir: " + workingDir); saver.setFile(new File(workingDir + "/" + fileName)); saver.writeBatch(); } catch (IOException e) { e.printStackTrace(); } }
From source file:net.djnn.addons.interpolation.CubicSplineInterpolation.java
public CubicSplineInterpolation(Element parent, String name, Element toMove, double[] t, double[] x) { super(parent, name); mu = new DoubleProperty(this, "mu", 0); SplineInterpolator sp = new SplineInterpolator(); psf = sp.interpolate(x, t); this.property = new DoubleProperty(toMove); new Binding(this, null, mu, new NativeAction(this, null, null, 1) { public void callback(Element e) { double muVal = mu.getValue(); if (muVal > 1) muVal = 1;//from ww w .jav a 2 s. c o m if (muVal < 0) muVal = 0; property.setValue(psf.value(muVal)); } }); }
From source file:de.treichels.hott.mdlviewer.swt.SwtCurveImageGenerator.java
private Image getImage(final Curve curve, final double scale, final boolean description) { // pitch curves start with 0% instead of -100% final boolean pitchCurve = curve.getPoint().get(0).getPosition() == 0; Image image;// w ww . j a v a 2 s. com GC g = null; try { // 200x250 pixel image with 5 pixel border scaled by factor scale image = new Image(Display.getDefault(), (int) (10 + 200 * scale), (int) (10 + 250 * scale)); g = new GC(image); g.setAntialias(SWT.ON); g.setTextAntialias(SWT.ON); // clear background g.setBackground(Display.getDefault().getSystemColor(SWT.COLOR_WHITE)); g.fillRectangle(0, 0, (int) (10 + 200 * scale), (int) (10 + 250 * scale)); // outer rectangle g.setForeground(Display.getDefault().getSystemColor(SWT.COLOR_BLACK)); g.setLineWidth(1); g.drawRectangle(5, 5, (int) (200 * scale), (int) (250 * scale)); g.setForeground(Display.getDefault().getSystemColor(SWT.COLOR_GRAY)); // +100% horizontal line g.drawLine(5, (int) (5 + 25 * scale), (int) (5 + 200 * scale), (int) (5 + 25 * scale)); // -100% horizontal line g.drawLine(5, (int) (5 + 225 * scale), (int) (5 + 200 * scale), (int) (5 + 225 * scale)); // 0% horizontal line g.drawLine(5, (int) (5 + 125 * scale), (int) (5 + 200 * scale), (int) (5 + 125 * scale)); if (!pitchCurve) { // 0% vertical line g.drawLine((int) (5 + 100 * scale), 5, (int) (5 + 100 * scale), (int) (5 + 250 * scale)); } g.setFont(new Font(Display.getDefault(), new FontData("Arial", 10, SWT.BOLD))); g.setForeground(Display.getDefault().getSystemColor(SWT.COLOR_BLACK)); // determine number of enabled curve points int numPoints = 0; for (final CurvePoint p : curve.getPoint()) if (p.isEnabled()) numPoints++; final double[] xVals = new double[numPoints]; final double[] yVals = new double[numPoints]; int i = 0; // store coordinates for (final CurvePoint p : curve.getPoint()) if (p.isEnabled()) { if (i == 0) // first point x coordinate is fixed to -100% (0% // for pitch curve) xVals[i] = pitchCurve ? 0 : -100; else if (i == numPoints - 1) // last point x coordinate is fixed to +100% xVals[i] = 100; else xVals[i] = p.getPosition(); yVals[i] = p.getValue(); if (description) { int x0; int y0; if (pitchCurve) { x0 = (int) (5 + xVals[i] * 2 * scale); y0 = (int) (5 + (225 - yVals[i] * 2) * scale); } else { x0 = (int) (5 + (100 + xVals[i]) * scale); y0 = (int) (5 + (125 - yVals[i]) * scale); } // draw point g.setBackground(Display.getDefault().getSystemColor(SWT.COLOR_BLACK)); g.fillOval(x0 - 3, y0 - 3, 6, 6); // draw text g.setBackground(Display.getDefault().getSystemColor(SWT.COLOR_WHITE)); g.fillRectangle(x0 - 5, y0 + 3, 9, 13); g.drawText(Integer.toString(p.getNumber() + 1), x0 - 4, y0 + 2, true); } i++; } g.setLineWidth(2); if (numPoints > 2 && curve.isSmoothing()) { // use a spline interpolator to smooth the curve final SplineInterpolator s = new SplineInterpolator(); final PolynomialSplineFunction function = s.interpolate(xVals, yVals); int x0 = 5; int y0; // starting point screen coordinates if (pitchCurve) y0 = (int) (5 + (225 - yVals[0] * 2) * scale); else y0 = (int) (5 + (125 - yVals[0]) * scale); // draw line pixel-by-pixel while (x0 < (int) (4 + 200 * scale)) { final int x1 = x0 + 1; int y1; if (pitchCurve) y1 = (int) (5 + (225 - function.value((x1 - 5) / scale / 2) * 2) * scale); else y1 = (int) (5 + (125 - function.value((x1 - 5) / scale - 100)) * scale); g.drawLine(x0, y0, x1, y1); x0 = x1; y0 = y1; } } else // draw line segments for (i = 0; i < numPoints - 1; i++) { int x0, y0, x1, y1; if (pitchCurve) { x0 = (int) (5 + xVals[i] * 2 * scale); y0 = (int) (5 + (225 - yVals[i] * 2) * scale); x1 = (int) (5 + xVals[i + 1] * 2 * scale); y1 = (int) (5 + (225 - yVals[i + 1] * 2) * scale); } else { x0 = (int) (5 + (100 + xVals[i]) * scale); y0 = (int) (5 + (125 - yVals[i]) * scale); x1 = (int) (5 + (100 + xVals[i + 1]) * scale); y1 = (int) (5 + (125 - yVals[i + 1]) * scale); } g.drawLine(x0, y0, x1, y1); } } finally { if (g != null) g.dispose(); } return image; }
From source file:de.treichels.hott.ui.android.html.AndroidCurveImageGenerator.java
private Bitmap getBitmap(final Curve curve, final float scale, final boolean description) { final boolean pitchCurve = curve.getPoint()[0].getPosition() == 0; final float scale1 = scale * 0.75f; // smaller images on the android // platform/*from ww w . j a v a 2s . c o m*/ final Bitmap image = Bitmap.createBitmap((int) (10 + 200 * scale1), (int) (10 + 250 * scale1), Bitmap.Config.ARGB_8888); final Canvas canvas = new Canvas(image); final Paint backgroundPaint = new Paint(); backgroundPaint.setColor(Color.WHITE); backgroundPaint.setStyle(Style.FILL); final Paint forgroundPaint = new Paint(); forgroundPaint.setColor(Color.BLACK); forgroundPaint.setStyle(Style.STROKE); forgroundPaint.setStrokeWidth(1.0f); forgroundPaint.setStrokeCap(Cap.BUTT); forgroundPaint.setStrokeJoin(Join.ROUND); forgroundPaint.setStrokeMiter(0.0f); final Paint curvePaint = new Paint(forgroundPaint); curvePaint.setFlags(Paint.ANTI_ALIAS_FLAG); curvePaint.setStrokeWidth(2.0f); final Paint pointPaint = new Paint(curvePaint); pointPaint.setStrokeWidth(5.0f); pointPaint.setStyle(Style.FILL_AND_STROKE); final Paint helpLinePaint = new Paint(forgroundPaint); helpLinePaint.setColor(Color.GRAY); helpLinePaint.setPathEffect(new DashPathEffect(new float[] { 5.0f, 5.0f }, 2.5f)); final Paint textPaint = new Paint(forgroundPaint); textPaint.setTypeface(Typeface.create(Typeface.SANS_SERIF, Typeface.BOLD)); textPaint.setTextSize(12.0f); textPaint.setTextAlign(Align.CENTER); textPaint.setStyle(Style.FILL); canvas.drawRect(0, 0, 10 + 200 * scale1, 10 + 250 * scale1, backgroundPaint); canvas.drawRect(5, 5, 5 + 200 * scale1, 5 + 250 * scale1, forgroundPaint); canvas.drawLine(5, 5 + 25 * scale1, 5 + 200 * scale1, 5 + 25 * scale1, helpLinePaint); canvas.drawLine(5, 5 + 225 * scale1, 5 + 200 * scale1, 5 + 225 * scale1, helpLinePaint); if (!pitchCurve) { canvas.drawLine(5, 5 + 125 * scale1, 5 + 200 * scale1, 5 + 125 * scale1, helpLinePaint); canvas.drawLine(5 + 100 * scale1, 5, 5 + 100 * scale1, 5 + 250 * scale1, helpLinePaint); } if (curve.getPoint() != null) { int numPoints = 0; for (final CurvePoint p : curve.getPoint()) { if (p.isEnabled()) { numPoints++; } } final double[] xVals = new double[numPoints]; final double[] yVals = new double[numPoints]; int i = 0; for (final CurvePoint p : curve.getPoint()) { if (p.isEnabled()) { if (i == 0) { xVals[i] = pitchCurve ? 0 : -100; } else if (i == numPoints - 1) { xVals[i] = 100; } else { xVals[i] = p.getPosition(); } yVals[i] = p.getValue(); if (description) { float x0; float y0; if (pitchCurve) { x0 = (float) (5 + xVals[i] * 2 * scale1); y0 = (float) (5 + (225 - yVals[i] * 2) * scale1); } else { x0 = (float) (5 + (100 + xVals[i]) * scale1); y0 = (float) (5 + (125 - yVals[i]) * scale1); } canvas.drawPoint(x0, y0, pointPaint); if (y0 < 5 + 125 * scale1) { canvas.drawRect(x0 - 4, y0 + 5, x0 + 3, y0 + 18, backgroundPaint); canvas.drawText(Integer.toString(p.getNumber() + 1), x0 - 1, y0 + 16, textPaint); } else { canvas.drawRect(x0 - 4, y0 - 5, x0 + 3, y0 - 18, backgroundPaint); canvas.drawText(Integer.toString(p.getNumber() + 1), x0 - 1, y0 - 7, textPaint); } } i++; } } if (numPoints > 2 && curve.isSmoothing()) { final SplineInterpolator s = new SplineInterpolator(); final PolynomialSplineFunction function = s.interpolate(xVals, yVals); float x0 = 5; float y0; if (pitchCurve) { y0 = (float) (5 + (225 - yVals[0] * 2) * scale1); } else { y0 = (float) (5 + (125 - yVals[0]) * scale1); } while (x0 < 4 + 200 * scale1) { final float x1 = x0 + 1; float y1; if (pitchCurve) { y1 = (float) (5 + (225 - function.value((x1 - 5) / scale1 / 2) * 2) * scale1); } else { y1 = (float) (5 + (125 - function.value((x1 - 5) / scale1 - 100)) * scale1); } canvas.drawLine(x0, y0, x1, y1, curvePaint); x0 = x1; y0 = y1; } } else { for (i = 0; i < numPoints - 1; i++) { float x0, y0, x1, y1; if (pitchCurve) { x0 = (float) (5 + xVals[i] * 2 * scale1); y0 = (float) (5 + (225 - yVals[i] * 2) * scale1); x1 = (float) (5 + xVals[i + 1] * 2 * scale1); y1 = (float) (5 + (225 - yVals[i + 1] * 2) * scale1); } else { x0 = (float) (5 + (100 + xVals[i]) * scale1); y0 = (float) (5 + (125 - yVals[i]) * scale1); x1 = (float) (5 + (100 + xVals[i + 1]) * scale1); y1 = (float) (5 + (125 - yVals[i + 1]) * scale1); } canvas.drawLine(x0, y0, x1, y1, curvePaint); } } } return image; }
From source file:de.biomedical_imaging.traj.math.TrajectorySplineFit.java
/** * Calculates a spline to a trajectory. Attention: The spline is fitted through a rotated version of the trajectory. * To fit the spline the trajectory is rotated into its main direction. You can access this rotated trajectory by * {@link #getRotatedTrajectory() getRotatedTrajectory}. * @return The fitted spline// w w w . j av a2s .c om */ public PolynomialSplineFunction calculateSpline() { successfull = false; /* * 1.Calculate the minimum bounding rectangle */ ArrayList<Point2D.Double> points = new ArrayList<Point2D.Double>(); for (int i = 0; i < t.size(); i++) { Point2D.Double p = new Point2D.Double(); p.setLocation(t.get(i).x, t.get(i).y); points.add(p); } /* * 1.1 Rotate that the major axis is parallel with the xaxis */ Array2DRowRealMatrix gyr = RadiusGyrationTensor2D.getRadiusOfGyrationTensor(t); EigenDecomposition eigdec = new EigenDecomposition(gyr); double inRad = -1 * Math.atan2(eigdec.getEigenvector(0).getEntry(1), eigdec.getEigenvector(0).getEntry(0)); boolean doTransform = (Math.abs(Math.abs(inRad) - Math.PI) > 0.001); if (doTransform) { angleRotated = inRad; for (int i = 0; i < t.size(); i++) { double x = t.get(i).x; double y = t.get(i).y; double newX = x * Math.cos(inRad) - y * Math.sin(inRad); double newY = x * Math.sin(inRad) + y * Math.cos(inRad); rotatedTrajectory.add(newX, newY, 0); points.get(i).setLocation(newX, newY); } //for(int i = 0; i < rect.length; i++){ // rect[i].setLocation(rect[i].x*Math.cos(inRad)-rect[i].y*Math.sin(inRad), rect[i].x*Math.sin(inRad)+rect[i].y*Math.cos(inRad)); //} } else { angleRotated = 0; rotatedTrajectory = t; } /* * 2. Divide the rectangle in n equal segments * 2.1 Calculate line in main direction * 2.2 Project the points in onto this line * 2.3 Calculate the distance between the start of the line and the projected point * 2.4 Assign point to segment according to distance of (2.3) */ List<List<Point2D.Double>> pointsInSegments = null; boolean allSegmentsContainingAtLeastTwoPoints = true; int indexSmallestX = 0; double segmentWidth = 0; do { double smallestX = Double.MAX_VALUE; double largestX = Double.MIN_VALUE; for (int i = 0; i < points.size(); i++) { if (points.get(i).x < smallestX) { smallestX = points.get(i).x; indexSmallestX = i; } if (points.get(i).x > largestX) { largestX = points.get(i).x; } } allSegmentsContainingAtLeastTwoPoints = true; segmentWidth = (largestX - smallestX) / nSegments; pointsInSegments = new ArrayList<List<Point2D.Double>>(nSegments); for (int i = 0; i < nSegments; i++) { pointsInSegments.add(new ArrayList<Point2D.Double>()); } for (int i = 0; i < points.size(); i++) { int index = (int) Math.abs((points.get(i).x / segmentWidth)); if (index > (nSegments - 1)) { index = (nSegments - 1); } pointsInSegments.get(index).add(points.get(i)); } for (int i = 0; i < pointsInSegments.size(); i++) { if (pointsInSegments.get(i).size() < 2) { if (nSegments > 2) { nSegments--; i = pointsInSegments.size(); allSegmentsContainingAtLeastTwoPoints = false; } } } } while (allSegmentsContainingAtLeastTwoPoints == false); /* * 3. Calculate the mean standard deviation over each segment: <s> */ //Point2D.Double eMajorP1 = new Point2D.Double(p1.x - (p3.x-p1.x)/2.0,p1.y - (p3.y-p1.y)/2.0); // Point2D.Double eMajorP2 = new Point2D.Double(p2.x - (p3.x-p1.x)/2.0,p2.y - (p3.y-p1.y)/2.0); double sumSDOrthogonal = 0; int Nsum = 0; CenterOfGravityFeature cogf = new CenterOfGravityFeature(rotatedTrajectory); Point2D.Double cog = new Point2D.Double(cogf.evaluate()[0], cogf.evaluate()[1]); Point2D.Double eMajorP1 = cog; Point2D.Double eMajorP2 = new Point2D.Double(cog.getX() + 1, cog.getY()); for (int i = 0; i < nSegments; i++) { StandardDeviation sd = new StandardDeviation(); double[] distances = new double[pointsInSegments.get(i).size()]; for (int j = 0; j < pointsInSegments.get(i).size(); j++) { int factor = 1; if (isLeft(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j))) { factor = -1; } distances[j] = factor * distancePointLine(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j)); } if (distances.length > 0) { sd.setData(distances); sumSDOrthogonal += sd.evaluate(); Nsum++; } } double s = sumSDOrthogonal / Nsum; if (segmentWidth < Math.pow(10, -15)) { return null; } if (s < Math.pow(10, -15)) { //If standard deviation is zero, replace it with the half of the segment width s = segmentWidth / 2; } //rotatedTrajectory.showTrajectory("rot"); /* * 4. Build a kd-tree */ KDTree<Point2D.Double> kdtree = new KDTree<Point2D.Double>(2); for (int i = 0; i < points.size(); i++) { try { //To ensure that all points have a different key, add small random number kdtree.insert(new double[] { points.get(i).x, points.get(i).y }, points.get(i)); } catch (KeySizeException e) { e.printStackTrace(); } catch (KeyDuplicateException e) { //Do nothing! It is not important } } /* * 5. Using the first point f in trajectory and calculate the center of mass * of all points around f (radius: 3*<s>)) */ List<Point2D.Double> near = null; Point2D.Double first = points.get(indexSmallestX);//minDistancePointToLine(p1, p3, points); double r1 = 3 * s; try { near = kdtree.nearestEuclidean(new double[] { first.x, first.y }, r1); } catch (KeySizeException e) { e.printStackTrace(); } double cx = 0; double cy = 0; for (int i = 0; i < near.size(); i++) { cx += near.get(i).x; cy += near.get(i).y; } cx /= near.size(); cy /= near.size(); splineSupportPoints = new ArrayList<Point2D.Double>(); splineSupportPoints.add(new Point2D.Double(cx, cy)); /* * 6. The second point is determined by finding the center-of-mass of particles in the p/2 radian * section of an annulus, r1 < r < 2r1, that is directed toward the angle with the highest number * of particles within p/2 radians. * 7. This second point is then used as the center of the annulus for choosing the third point, and the process is repeated (6. & 7.). */ /* * 6.1 Find all points in the annolous */ /* * 6.2 Write each point in a coordinate system centered at the center of the sphere, calculate direction and * check if it in the allowed bounds */ int nCircleSegments = 100; double deltaRad = 2 * Math.PI / nCircleSegments; boolean stop = false; int minN = 7; double tempr1 = r1; double allowedDeltaDirection = 0.5 * Math.PI; while (stop == false) { List<Point2D.Double> nearestr1 = null; List<Point2D.Double> nearest2xr1 = null; try { nearestr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, tempr1); nearest2xr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, 2 * tempr1); } catch (KeySizeException e) { // TODO Auto-generated catch block e.printStackTrace(); } nearest2xr1.removeAll(nearestr1); double lThreshRad = 0; double hThreshRad = Math.PI / 2; double stopThresh = 2 * Math.PI; if (splineSupportPoints.size() > 1) { double directionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; lThreshRad = directionInRad - allowedDeltaDirection / 2 - Math.PI / 4; if (lThreshRad < 0) { lThreshRad = 2 * Math.PI + lThreshRad; } if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } hThreshRad = directionInRad + allowedDeltaDirection / 2 + Math.PI / 4; if (hThreshRad < 0) { hThreshRad = 2 * Math.PI + hThreshRad; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } stopThresh = directionInRad + allowedDeltaDirection / 2 - Math.PI / 4; if (stopThresh > 2 * Math.PI) { stopThresh = stopThresh - 2 * Math.PI; } } double newCx = 0; double newCy = 0; int newCN = 0; int candN = 0; //Find center with highest density of points double lastDist = 0; double newDist = 0; do { lastDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); candN = 0; double candCx = 0; double candCy = 0; for (int i = 0; i < nearest2xr1.size(); i++) { Point2D.Double centerOfCircle = splineSupportPoints.get(splineSupportPoints.size() - 1); Vector2d relativeToCircle = new Vector2d(nearest2xr1.get(i).x - centerOfCircle.x, nearest2xr1.get(i).y - centerOfCircle.y); relativeToCircle.normalize(); double angleInRadians = Math.atan2(relativeToCircle.y, relativeToCircle.x) + Math.PI; if (lThreshRad < hThreshRad) { if (angleInRadians > lThreshRad && angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } else { if (angleInRadians > lThreshRad || angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } } if (candN > 0 && candN > newCN) { candCx /= candN; candCy /= candN; newCx = candCx; newCy = candCy; newCN = candN; } lThreshRad += deltaRad; hThreshRad += deltaRad; if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } newDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); } while ((newDist - lastDist) > 0); //Check if the new center is valid if (splineSupportPoints.size() > 1) { double currentDirectionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; double candDirectionInRad = Math.atan2( newCy - splineSupportPoints.get(splineSupportPoints.size() - 1).y, newCx - splineSupportPoints.get(splineSupportPoints.size() - 1).x) + Math.PI; double dDir = Math.max(currentDirectionInRad, candDirectionInRad) - Math.min(currentDirectionInRad, candDirectionInRad); if (dDir > 2 * Math.PI) { dDir = 2 * Math.PI - dDir; } if (dDir > allowedDeltaDirection) { stop = true; } } boolean enoughPoints = (newCN < minN); boolean isNormalRadius = Math.abs(tempr1 - r1) < Math.pow(10, -18); boolean isExtendedRadius = Math.abs(tempr1 - 3 * r1) < Math.pow(10, -18); if (enoughPoints && isNormalRadius) { //Not enough points, extend search radius tempr1 = 3 * r1; } else if (enoughPoints && isExtendedRadius) { //Despite radius extension: Not enough points! stop = true; } else if (stop == false) { splineSupportPoints.add(new Point2D.Double(newCx, newCy)); tempr1 = r1; } } //Sort Collections.sort(splineSupportPoints, new Comparator<Point2D.Double>() { public int compare(Point2D.Double o1, Point2D.Double o2) { if (o1.x < o2.x) { return -1; } if (o1.x > o2.x) { return 1; } return 0; } }); //Add endpoints if (splineSupportPoints.size() > 1) { Vector2d start = new Vector2d(splineSupportPoints.get(0).x - splineSupportPoints.get(1).x, splineSupportPoints.get(0).y - splineSupportPoints.get(1).y); start.normalize(); start.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + start.x, splineSupportPoints.get(0).y + start.y)); Vector2d end = new Vector2d( splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y); end.normalize(); end.scale(r1 * 6); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + end.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + end.y)); } else { Vector2d majordir = new Vector2d(-1, 0); majordir.normalize(); majordir.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + majordir.x, splineSupportPoints.get(0).y + majordir.y)); majordir.scale(-1); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + majordir.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + majordir.y)); } //Interpolate spline double[] supX = new double[splineSupportPoints.size()]; double[] supY = new double[splineSupportPoints.size()]; for (int i = 0; i < splineSupportPoints.size(); i++) { supX[i] = splineSupportPoints.get(i).x; supY[i] = splineSupportPoints.get(i).y; } SplineInterpolator sIinter = new SplineInterpolator(); spline = sIinter.interpolate(supX, supY); successfull = true; return spline; }
From source file:de.biomedical_imaging.traj.math.TrajectorySplineFitLegacy.java
/** * Calculates a spline to a trajectory. Attention: The spline is fitted through a rotated version of the trajectory. * To fit the spline the trajectory is rotated into its main direction. You can access this rotated trajectory by * {@link #getRotatedTrajectory() getRotatedTrajectory}. * @return The fitted spline//w w w . j ava 2s.c om */ public PolynomialSplineFunction calculateSpline() { /* * 1.Calculate the minimum bounding rectangle */ ArrayList<Point2D.Double> points = new ArrayList<Point2D.Double>(); for (int i = 0; i < t.size(); i++) { Point2D.Double p = new Point2D.Double(); p.setLocation(t.get(i).x, t.get(i).y); points.add(p); } Point2D.Double[] rect = null; try { rect = RotatingCalipers.getMinimumBoundingRectangle(points); } catch (IllegalArgumentException e) { } catch (EmptyStackException e) { } /* * 1.1 Rotate that the major axis is parallel with the xaxis */ Point2D.Double majorDirection = null; Point2D.Double p1 = rect[2]; //top left Point2D.Double p2 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[1] : rect[3]; //Point to long side Point2D.Double p3 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[3] : rect[1]; //Point to short side majorDirection = new Point2D.Double(p2.x - p1.x, p2.y - p1.y); double width = p1.distance(p2); double inRad = -1 * Math.atan2(majorDirection.y, majorDirection.x); boolean doTransform = (Math.abs(Math.abs(inRad) - Math.PI) > 0.001); if (doTransform) { angleRotated = inRad; for (int i = 0; i < t.size(); i++) { double x = t.get(i).x; double y = t.get(i).y; double newX = x * Math.cos(inRad) - y * Math.sin(inRad); double newY = x * Math.sin(inRad) + y * Math.cos(inRad); rotatedTrajectory.add(newX, newY, 0); points.get(i).setLocation(newX, newY); } for (int i = 0; i < rect.length; i++) { rect[i].setLocation(rect[i].x * Math.cos(inRad) - rect[i].y * Math.sin(inRad), rect[i].x * Math.sin(inRad) + rect[i].y * Math.cos(inRad)); } p1 = rect[2]; //top left p2 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[1] : rect[3]; //Point to long side p3 = p1.distance(rect[1]) > p1.distance(rect[3]) ? rect[3] : rect[1]; //Point to short side } else { angleRotated = 0; rotatedTrajectory = t; } /* * 2. Divide the rectangle in n equal segments * 2.1 Calculate line in main direction * 2.2 Project the points in onto this line * 2.3 Calculate the distance between the start of the line and the projected point * 2.4 Assign point to segment according to distance of (2.3) */ List<List<Point2D.Double>> pointsInSegments = null; boolean allSegmentsContainingAtLeastTwoPoints = true; do { allSegmentsContainingAtLeastTwoPoints = true; double segmentWidth = p1.distance(p2) / nSegments; pointsInSegments = new ArrayList<List<Point2D.Double>>(nSegments); for (int i = 0; i < nSegments; i++) { pointsInSegments.add(new ArrayList<Point2D.Double>()); } for (int i = 0; i < points.size(); i++) { Point2D.Double projPoint = projectPointToLine(p1, p2, points.get(i)); int index = (int) (p1.distance(projPoint) / segmentWidth); if (index > (nSegments - 1)) { index = (nSegments - 1); } pointsInSegments.get(index).add(points.get(i)); } for (int i = 0; i < pointsInSegments.size(); i++) { if (pointsInSegments.get(i).size() < 2) { if (nSegments > 2) { nSegments--; i = pointsInSegments.size(); allSegmentsContainingAtLeastTwoPoints = false; } } } } while (allSegmentsContainingAtLeastTwoPoints == false); /* * 3. Calculate the mean standard deviation over each segment: <s> */ Point2D.Double eMajorP1 = new Point2D.Double(p1.x - (p3.x - p1.x) / 2.0, p1.y - (p3.y - p1.y) / 2.0); Point2D.Double eMajorP2 = new Point2D.Double(p2.x - (p3.x - p1.x) / 2.0, p2.y - (p3.y - p1.y) / 2.0); double sumMean = 0; int Nsum = 0; for (int i = 0; i < nSegments; i++) { StandardDeviation sd = new StandardDeviation(); double[] distances = new double[pointsInSegments.get(i).size()]; for (int j = 0; j < pointsInSegments.get(i).size(); j++) { int factor = 1; if (isLeft(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j))) { factor = -1; } distances[j] = factor * distancePointLine(eMajorP1, eMajorP2, pointsInSegments.get(i).get(j)); } if (distances.length > 0) { sd.setData(distances); sumMean += sd.evaluate(); Nsum++; } } double s = sumMean / Nsum; if (s < 0.000000000001) { s = width / nSegments; } /* * 4. Build a kd-tree */ KDTree<Point2D.Double> kdtree = new KDTree<Point2D.Double>(2); for (int i = 0; i < points.size(); i++) { try { //To ensure that all points have a different key, add small random number kdtree.insert(new double[] { points.get(i).x, points.get(i).y }, points.get(i)); } catch (KeySizeException e) { e.printStackTrace(); } catch (KeyDuplicateException e) { //Do nothing! It is not important } } /* * 5. Using the first point f in trajectory and calculate the center of mass * of all points around f (radius: 3*<s>)) */ List<Point2D.Double> near = null; Point2D.Double first = minDistancePointToLine(p1, p3, points); double r1 = 3 * s; try { near = kdtree.nearestEuclidean(new double[] { first.x, first.y }, r1); } catch (KeySizeException e) { e.printStackTrace(); } double cx = 0; double cy = 0; for (int i = 0; i < near.size(); i++) { cx += near.get(i).x; cy += near.get(i).y; } cx /= near.size(); cy /= near.size(); splineSupportPoints = new ArrayList<Point2D.Double>(); splineSupportPoints.add(new Point2D.Double(cx, cy)); /* * 6. The second point is determined by finding the center-of-mass of particles in the p/2 radian * section of an annulus, r1 < r < 2r1, that is directed toward the angle with the highest number * of particles within p/2 radians. * 7. This second point is then used as the center of the annulus for choosing the third point, and the process is repeated (6. & 7.). */ /* * 6.1 Find all points in the annolous */ /* * 6.2 Write each point in a coordinate system centered at the center of the sphere, calculate direction and * check if it in the allowed bounds */ int nCircleSegments = 100; double deltaRad = 2 * Math.PI / nCircleSegments; boolean stop = false; int minN = 7; double tempr1 = r1; double allowedDeltaDirection = 0.5 * Math.PI; while (stop == false) { List<Point2D.Double> nearestr1 = null; List<Point2D.Double> nearest2xr1 = null; try { nearestr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, tempr1); nearest2xr1 = kdtree .nearestEuclidean(new double[] { splineSupportPoints.get(splineSupportPoints.size() - 1).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y }, 2 * tempr1); } catch (KeySizeException e) { // TODO Auto-generated catch block e.printStackTrace(); } nearest2xr1.removeAll(nearestr1); double lThreshRad = 0; double hThreshRad = Math.PI / 2; double stopThresh = 2 * Math.PI; if (splineSupportPoints.size() > 1) { double directionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; lThreshRad = directionInRad - allowedDeltaDirection / 2 - Math.PI / 4; if (lThreshRad < 0) { lThreshRad = 2 * Math.PI + lThreshRad; } if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } hThreshRad = directionInRad + allowedDeltaDirection / 2 + Math.PI / 4; if (hThreshRad < 0) { hThreshRad = 2 * Math.PI + hThreshRad; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } stopThresh = directionInRad + allowedDeltaDirection / 2 - Math.PI / 4; if (stopThresh > 2 * Math.PI) { stopThresh = stopThresh - 2 * Math.PI; } } double newCx = 0; double newCy = 0; int newCN = 0; int candN = 0; //Find center with highest density of points double lastDist = 0; double newDist = 0; do { lastDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); candN = 0; double candCx = 0; double candCy = 0; for (int i = 0; i < nearest2xr1.size(); i++) { Point2D.Double centerOfCircle = splineSupportPoints.get(splineSupportPoints.size() - 1); Vector2d relativeToCircle = new Vector2d(nearest2xr1.get(i).x - centerOfCircle.x, nearest2xr1.get(i).y - centerOfCircle.y); relativeToCircle.normalize(); double angleInRadians = Math.atan2(relativeToCircle.y, relativeToCircle.x) + Math.PI; if (lThreshRad < hThreshRad) { if (angleInRadians > lThreshRad && angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } else { if (angleInRadians > lThreshRad || angleInRadians < hThreshRad) { candCx += nearest2xr1.get(i).x; candCy += nearest2xr1.get(i).y; candN++; } } } if (candN > 0 && candN > newCN) { candCx /= candN; candCy /= candN; newCx = candCx; newCy = candCy; newCN = candN; } lThreshRad += deltaRad; hThreshRad += deltaRad; if (lThreshRad > 2 * Math.PI) { lThreshRad = lThreshRad - 2 * Math.PI; } if (hThreshRad > 2 * Math.PI) { hThreshRad = hThreshRad - 2 * Math.PI; } newDist = Math.min(Math.abs(lThreshRad - stopThresh), 2 * Math.PI - Math.abs(lThreshRad - stopThresh)); } while ((newDist - lastDist) > 0); //Check if the new center is valid if (splineSupportPoints.size() > 1) { double currentDirectionInRad = Math.atan2( splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y, splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x) + Math.PI; double candDirectionInRad = Math.atan2( newCy - splineSupportPoints.get(splineSupportPoints.size() - 1).y, newCx - splineSupportPoints.get(splineSupportPoints.size() - 1).x) + Math.PI; double dDir = Math.max(currentDirectionInRad, candDirectionInRad) - Math.min(currentDirectionInRad, candDirectionInRad); if (dDir > 2 * Math.PI) { dDir = 2 * Math.PI - dDir; } if (dDir > allowedDeltaDirection) { stop = true; } } boolean enoughPoints = (newCN < minN); boolean isNormalRadius = Math.abs(tempr1 - r1) < Math.pow(10, -18); boolean isExtendedRadius = Math.abs(tempr1 - 3 * r1) < Math.pow(10, -18); if (enoughPoints && isNormalRadius) { //Not enough points, extend search radius tempr1 = 3 * r1; } else if (enoughPoints && isExtendedRadius) { //Despite radius extension: Not enough points! stop = true; } else if (stop == false) { splineSupportPoints.add(new Point2D.Double(newCx, newCy)); tempr1 = r1; } } //Sort Collections.sort(splineSupportPoints, new Comparator<Point2D.Double>() { public int compare(Point2D.Double o1, Point2D.Double o2) { if (o1.x < o2.x) { return -1; } if (o1.x > o2.x) { return 1; } return 0; } }); //Add endpoints if (splineSupportPoints.size() > 1) { Vector2d start = new Vector2d(splineSupportPoints.get(0).x - splineSupportPoints.get(1).x, splineSupportPoints.get(0).y - splineSupportPoints.get(1).y); start.normalize(); start.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + start.x, splineSupportPoints.get(0).y + start.y)); Vector2d end = new Vector2d( splineSupportPoints.get(splineSupportPoints.size() - 1).x - splineSupportPoints.get(splineSupportPoints.size() - 2).x, splineSupportPoints.get(splineSupportPoints.size() - 1).y - splineSupportPoints.get(splineSupportPoints.size() - 2).y); end.normalize(); end.scale(r1 * 6); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + end.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + end.y)); } else { Vector2d majordir = new Vector2d(-1, 0); majordir.normalize(); majordir.scale(r1 * 8); splineSupportPoints.add(0, new Point2D.Double(splineSupportPoints.get(0).x + majordir.x, splineSupportPoints.get(0).y + majordir.y)); majordir.scale(-1); splineSupportPoints .add(new Point2D.Double(splineSupportPoints.get(splineSupportPoints.size() - 1).x + majordir.x, splineSupportPoints.get(splineSupportPoints.size() - 1).y + majordir.y)); } //Interpolate spline double[] supX = new double[splineSupportPoints.size()]; double[] supY = new double[splineSupportPoints.size()]; for (int i = 0; i < splineSupportPoints.size(); i++) { supX[i] = splineSupportPoints.get(i).x; supY[i] = splineSupportPoints.get(i).y; } SplineInterpolator sIinter = new SplineInterpolator(); spline = sIinter.interpolate(supX, supY); return spline; }
From source file:be.ugent.maf.cellmissy.analysis.singlecell.processing.impl.interpolation.TrackSplineInterpolator.java
@Override public InterpolatedTrack interpolateTrack(double[] time, double[] x, double[] y) { // create a new spline interpolator SplineInterpolator splineInterpolator = new SplineInterpolator(); int interpolationPoints = PropertiesConfigurationHolder.getInstance().getInt("numberOfInterpolationPoints"); // create arrays to hold the interpolant time, the interpolated X and the interpolated Y double[] interpolantTime = new double[interpolationPoints]; double[] interpolatedX = new double[interpolationPoints]; double[] interpolatedY = new double[interpolationPoints]; // the step used for the interpolation in both direction double interpolationStep = (time[time.length - 1] - time[0]) / interpolationPoints; // check for monotonicity boolean monotonic = MathArrays.isMonotonic(time, MathArrays.OrderDirection.INCREASING, false); // in case time is not monotonic, sort in place time, x and y coordinates if (!monotonic) { MathArrays.sortInPlace(time, x, y); }// w ww .j a va2 s. com // call the interpolator, and actually do the interpolation try { PolynomialSplineFunction functionX = splineInterpolator.interpolate(time, x); PolynomialSplineFunction functionY = splineInterpolator.interpolate(time, y); // get the polynomial functions in both directions PolynomialFunction polynomialFunctionX = functionX.getPolynomials()[0]; PolynomialFunction polynomialFunctionY = functionY.getPolynomials()[0]; for (int i = 0; i < interpolationPoints; i++) { interpolantTime[i] = time[0] + (i * interpolationStep); interpolatedX[i] = functionX.value(interpolantTime[i]); interpolatedY[i] = functionY.value(interpolantTime[i]); } for (int k = 0; k < interpolationPoints; k++) { if (Double.isNaN(interpolatedX[k]) | Double.isNaN(interpolatedY[k])) { return null; } } return new InterpolatedTrack(interpolantTime, interpolatedX, interpolatedY, polynomialFunctionX, polynomialFunctionY); } catch (NumberIsTooSmallException e) { LOG.error(e.getMessage()); return null; } }