Example usage for org.apache.commons.math3.analysis.interpolation SplineInterpolator SplineInterpolator

List of usage examples for org.apache.commons.math3.analysis.interpolation SplineInterpolator SplineInterpolator

Introduction

In this page you can find the example usage for org.apache.commons.math3.analysis.interpolation SplineInterpolator SplineInterpolator.

Prototype

SplineInterpolator

Source Link

Usage

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/*w w  w  . java2  s  .  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:cz.cuni.lf1.lge.ThunderSTORM.results.ModifiedLoess.java

/**
 * Compute an interpolating function by performing a loess fit
 * on the data at the original abscissae and then building a cubic spline
 * with a//from  w ww  .ja  v  a  2s .  c om
 * {@link org.apache.commons.math3.analysis.interpolation.SplineInterpolator}
 * on the resulting fit.
 *
 * @param xval the arguments for the interpolation points
 * @param yval the values for the interpolation points
 * @return A cubic spline built upon a loess fit to the data at the original abscissae
 * @throws NonMonotonicSequenceException if {@code xval} not sorted in
 * strictly increasing order.
 * @throws DimensionMismatchException if {@code xval} and {@code yval} have
 * different sizes.
 * @throws NoDataException if {@code xval} or {@code yval} has zero size.
 * @throws NotFiniteNumberException if any of the arguments and values are
 * not finite real numbers.
 * @throws NumberIsTooSmallException if the bandwidth is too small to
 * accomodate the size of the input data (i.e. the bandwidth must be
 * larger than 2/n).
 */
public final PolynomialSplineFunction interpolate(double[] xval, double[] yval)
        throws NonMonotonicSequenceException, DimensionMismatchException, NoDataException,
        NotFiniteNumberException, NumberIsTooSmallException {
    double[] smoothed = smooth(xval, yval);
    DoubleList newX = new ArrayDoubleList();
    DoubleList newSmoothed = new ArrayDoubleList();
    newX.add(xval[0]);
    newSmoothed.add(smoothed[0]);
    for (int i = 1; i < xval.length; i++) {
        if (xval[i] != xval[i - 1]) {
            newX.add(xval[i]);
            newSmoothed.add(smoothed[i]);
        }
    }
    xval = newX.toArray();
    smoothed = newSmoothed.toArray();

    return new SplineInterpolator().interpolate(xval, smoothed);
}

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;/*from  w  w w.j  a v  a 2  s. c o m*/
    }
    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:com.bwc.ora.models.Lrp.java

/**
 * Get the local maximums from a collection of Points.
 *
 * @param lrpSeries//from   w  w w . j a va  2  s .  c  o  m
 * @param seriesTitle
 * @return
 */
public static XYSeries getMaximumsWithHiddenPeaks(XYSeries lrpSeries, String seriesTitle) {
    XYSeries maxPoints = new XYSeries(seriesTitle);

    //convert to x and y coordinate arrays
    double[][] xyline = lrpSeries.toArray();

    //use a spline interpolator to converts points into an equation
    UnivariateInterpolator interpolator = new SplineInterpolator();
    UnivariateFunction function = interpolator.interpolate(xyline[0], xyline[1]);

    // create a differentiator using 5 points and 0.01 step
    FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(5, 0.01);

    // create a new function that computes both the value and the derivatives
    // using DerivativeStructure
    UnivariateDifferentiableFunction completeF = differentiator.differentiate(function);

    // now we can compute the value and its derivatives
    // here we decided to display up to second order derivatives,
    // because we feed completeF with order 2 DerivativeStructure instances
    //find local minima in second derivative, these indicate the peaks (and hidden peaks)
    //of the input
    for (double x = xyline[0][0] + 1; x < xyline[0][xyline[0].length - 1] - 1; x += 0.5) {
        DerivativeStructure xDSc = new DerivativeStructure(1, 2, 0, x);
        DerivativeStructure xDSl = new DerivativeStructure(1, 2, 0, x - 0.5);
        DerivativeStructure xDSr = new DerivativeStructure(1, 2, 0, x + 0.5);
        DerivativeStructure yDSc = completeF.value(xDSc);
        DerivativeStructure yDSl = completeF.value(xDSl);
        DerivativeStructure yDSr = completeF.value(xDSr);
        double c2d = yDSc.getPartialDerivative(2);
        if (c2d < yDSl.getPartialDerivative(2) && c2d < yDSr.getPartialDerivative(2)) {
            maxPoints.add((int) Math.round(x), yDSc.getValue());
        }
    }

    return maxPoints;
}

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;/* ww w .  j  a  v  a2 s .  co m*/
    }
    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: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 ww. j  a  v a  2  s.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: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/*from  ww  w .ja v  a 2s .  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:gdsc.smlm.ij.plugins.DriftCalculator.java

private void interpolate(double[] dx, double[] dy, double[] originalDriftTimePoints) {
    // Interpolator can only create missing values within the range provided by the input values.
    // The two ends have to be extrapolated.
    // TODO: Perform extrapolation. Currently the end values are used.

    // Find end points
    int startT = 0;
    while (originalDriftTimePoints[startT] == 0)
        startT++;// w w  w.  ja  v a  2s  .  c  om
    int endT = originalDriftTimePoints.length - 1;
    while (originalDriftTimePoints[endT] == 0)
        endT--;

    // Extrapolate using a constant value
    for (int t = startT; t-- > 0;) {
        dx[t] = dx[startT];
        dy[t] = dy[startT];
    }
    for (int t = endT; ++t < dx.length;) {
        dx[t] = dx[endT];
        dy[t] = dy[endT];
    }

    double[][] values = extractValues(originalDriftTimePoints, startT, endT, dx, dy);

    PolynomialSplineFunction fx, fy;
    if (values[0].length < 3) {
        fx = new LinearInterpolator().interpolate(values[0], values[1]);
        fy = new LinearInterpolator().interpolate(values[0], values[2]);
    } else {
        fx = new SplineInterpolator().interpolate(values[0], values[1]);
        fy = new SplineInterpolator().interpolate(values[0], values[2]);
    }

    for (int t = startT; t <= endT; t++) {
        if (originalDriftTimePoints[t] == 0) {
            dx[t] = fx.value(t);
            dy[t] = fy.value(t);
        }
    }

    this.interpolationStart = startT;
    this.interpolationEnd = endT;
}

From source file:gdsc.smlm.model.AiryPSFModel.java

private static synchronized void createAiryDistribution() {
    if (spline != null)
        return;//www . j  a v  a 2  s  . c  om

    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:com.github.vatbub.tictactoe.view.Main.java

private void updateAILevelLabel(boolean forceUpdate) {
    double sliderPos = 100 * Math.round(aiLevelSlider.getValue() * 3.0 / 100.0) / 3.0;

    if (sliderPos != aiLevelLabelPositionProperty.get() || forceUpdate) {
        aiLevelLabelPositionProperty.set(sliderPos);

        // request focus of the current ai label for accessibility
        aiLevelLabelHBox.getChildren().get((int) (sliderPos * 3 / 100)).requestFocus();
        updateAccessibleTexts();//from w w  w .j  a  va  2 s .c  o m

        // get the slider position
        double[] xDouble = new double[] { 0, 100.0 / 3.0, 200.0 / 3.0, 300.0 / 3.0 };
        double[] translationYDouble = new double[4];
        double[] widthYDouble = new double[4];
        double[] trueWidthYDouble = new double[4];
        for (int i = 0; i < translationYDouble.length; i++) {
            // {-getAILevelLabelCenter(0), -getAILevelLabelCenter(1), -getAILevelLabelCenter(2), -getAILevelLabelCenter(3)};
            translationYDouble[i] = -getAILevelLabelCenter(i);
            widthYDouble[i] = Math.max(90, ((Label) aiLevelLabelHBox.getChildren().get(i)).getWidth()
                    + 8 * aiLevelLabelHBox.getSpacing());
            trueWidthYDouble[i] = ((Label) aiLevelLabelHBox.getChildren().get(i)).getWidth();
        }

        SplineInterpolator splineInterpolator = new SplineInterpolator();
        PolynomialSplineFunction translateFunction = splineInterpolator.interpolate(xDouble,
                translationYDouble);
        PolynomialSplineFunction widthFunction = splineInterpolator.interpolate(xDouble, widthYDouble);
        PolynomialSplineFunction trueWidthFunction = splineInterpolator.interpolate(xDouble, trueWidthYDouble);

        KeyValue hBoxLayoutXKeyValue1 = new KeyValue(aiLevelLabelHBox.layoutXProperty(),
                aiLevelLabelHBox.getLayoutX(), Interpolator.EASE_BOTH);
        KeyValue aiLevelLabelClipRectangleWidthKeyValue1 = new KeyValue(
                aiLevelLabelClipRectangle.widthProperty(), aiLevelLabelClipRectangle.getWidth(),
                Interpolator.EASE_BOTH);
        KeyValue aiLevelLabelClipRectangleXKeyValue1 = new KeyValue(aiLevelLabelClipRectangle.xProperty(),
                aiLevelLabelClipRectangle.getX(), Interpolator.EASE_BOTH);
        KeyValue aiLevelCenterLineStartXKeyValue1 = new KeyValue(aiLevelCenterLine.startXProperty(),
                aiLevelCenterLine.getStartX(), Interpolator.EASE_BOTH);
        KeyValue aiLevelCenterLineEndXKeyValue1 = new KeyValue(aiLevelCenterLine.endXProperty(),
                aiLevelCenterLine.getEndX(), Interpolator.EASE_BOTH);
        KeyFrame keyFrame1 = new KeyFrame(Duration.seconds(0), hBoxLayoutXKeyValue1,
                aiLevelLabelClipRectangleWidthKeyValue1, aiLevelLabelClipRectangleXKeyValue1,
                aiLevelCenterLineStartXKeyValue1, aiLevelCenterLineEndXKeyValue1);

        double interpolatedLabelWidth = trueWidthFunction.value(sliderPos);

        KeyValue hBoxLayoutXKeyValue2 = new KeyValue(aiLevelLabelHBox.layoutXProperty(),
                translateFunction.value(sliderPos), Interpolator.EASE_BOTH);
        KeyValue aiLevelLabelClipRectangleWidthKeyValue2 = new KeyValue(
                aiLevelLabelClipRectangle.widthProperty(), widthFunction.value(sliderPos),
                Interpolator.EASE_BOTH);
        KeyValue aiLevelLabelClipRectangleXKeyValue2 = new KeyValue(aiLevelLabelClipRectangle.xProperty(),
                aiLevelLabelPane.getWidth() / 2 - widthFunction.value(sliderPos) / 2, Interpolator.EASE_BOTH);
        KeyValue aiLevelCenterLineStartXKeyValue2 = new KeyValue(aiLevelCenterLine.startXProperty(),
                (aiLevelLabelPane.getWidth() - interpolatedLabelWidth) / 2, Interpolator.EASE_BOTH);
        KeyValue aiLevelCenterLineEndXKeyValue2 = new KeyValue(aiLevelCenterLine.endXProperty(),
                (aiLevelLabelPane.getWidth() + interpolatedLabelWidth) / 2, Interpolator.EASE_BOTH);
        KeyFrame keyFrame2 = new KeyFrame(Duration.seconds(animationSpeed * 0.9), hBoxLayoutXKeyValue2,
                aiLevelLabelClipRectangleWidthKeyValue2, aiLevelLabelClipRectangleXKeyValue2,
                aiLevelCenterLineStartXKeyValue2, aiLevelCenterLineEndXKeyValue2);

        Timeline timeline = new Timeline(keyFrame1, keyFrame2);
        timeline.play();
    }
}