Author: erans
Date: Thu Apr 26 12:13:49 2012
New Revision: 1330801
URL: http://svn.apache.org/viewvc?rev=1330801&view=rev
Log:
Added unit test.
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
(with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
(with props)
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java?rev=1330801&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
Thu Apr 26 12:13:49 2012
@@ -0,0 +1,173 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.optimization.general;
+
+import java.util.Arrays;
+import java.util.ArrayList;
+import
org.apache.commons.math3.analysis.DifferentiableMultivariateVectorFunction;
+import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
+import org.apache.commons.math3.util.MathUtils;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * Class that models a circle.
+ * The parameters of problem are:
+ * <ul>
+ * <li>the x-coordinate of the circle center,</li>
+ * <li>the y-coordinate of the circle center,</li>
+ * <li>the radius of the circle.</li>
+ * </ul>
+ * The model functions are:
+ * <ul>
+ * <li>for each triplet (cx, cy, r), the (x, y) coordinates of a point on the
+ * corresponding circle.</li>
+ * </ul>
+ */
+class CircleProblem implements DifferentiableMultivariateVectorFunction {
+ /** Cloud of points assumed to be fitted by a circle. */
+ private final ArrayList<double[]> points;
+ /** Error on the x-coordinate of the points. */
+ private final double xSigma;
+ /** Error on the y-coordinate of the points. */
+ private final double ySigma;
+ /** Number of points on the circumference (when searching which
+ model point is closest to a given "observation". */
+ private final int resolution;
+
+ /**
+ * @param xError Assumed error for the x-coordinate of the circle points.
+ * @param xError Assumed error for the y-coordinate of the circle points.
+ * @param searchResolution Number of points to try when searching the one
+ * that is closest to a given "observed" point.
+ */
+ public CircleProblem(double xError,
+ double yError,
+ int searchResolution) {
+ points = new ArrayList<double[]>();
+ xSigma = xError;
+ ySigma = yError;
+ resolution = searchResolution;
+ }
+
+ /**
+ * @param xError Assumed error for the x-coordinate of the circle points.
+ * @param xError Assumed error for the y-coordinate of the circle points.
+ */
+ public CircleProblem(double xError,
+ double yError) {
+ this(xError, yError, 500);
+ }
+
+ public void addPoint(double px, double py) {
+ points.add(new double[] { px, py });
+ }
+
+ public double[] target() {
+ final double[] t = new double[points.size() * 2];
+ for (int i = 0; i < points.size(); i++) {
+ final double[] p = points.get(i);
+ final int index = i * 2;
+ t[index] = p[0];
+ t[index + 1] = p[1];
+ }
+
+ return t;
+ }
+
+ public double[] weight() {
+ final double wX = 1 / (xSigma * xSigma);
+ final double wY = 1 / (ySigma * ySigma);
+ final double[] w = new double[points.size() * 2];
+ for (int i = 0; i < points.size(); i++) {
+ final int index = i * 2;
+ w[index] = wX;
+ w[index + 1] = wY;
+ }
+
+ return w;
+ }
+
+ public double[] value(double[] params) {
+ final double cx = params[0];
+ final double cy = params[1];
+ final double r = params[2];
+
+ final double[] model = new double[points.size() * 2];
+
+ final double deltaTheta = MathUtils.TWO_PI / resolution;
+ for (int i = 0; i < points.size(); i++) {
+ final double[] p = points.get(i);
+ final double px = p[0];
+ final double py = p[1];
+
+ double bestX = 0;
+ double bestY = 0;
+ double dMin = Double.POSITIVE_INFINITY;
+
+ // Find the angle for which the circle passes closest to the
+ // current point (using a resolution of 100 points along the
+ // circumference).
+ for (double theta = 0; theta <= MathUtils.TWO_PI; theta +=
deltaTheta) {
+ final double currentX = cx + r * FastMath.cos(theta);
+ final double currentY = cy + r * FastMath.sin(theta);
+ final double dX = currentX - px;
+ final double dY = currentY - py;
+ final double d = dX * dX + dY * dY;
+ if (d < dMin) {
+ dMin = d;
+ bestX = currentX;
+ bestY = currentY;
+ }
+ }
+
+ final int index = i * 2;
+ model[index] = bestX;
+ model[index + 1] = bestY;
+ }
+
+ return model;
+ }
+
+ public MultivariateMatrixFunction jacobian() {
+ return new MultivariateMatrixFunction() {
+ public double[][] value(double[] point) {
+ return jacobian(point);
+ }
+ };
+ }
+
+ private double[][] jacobian(double[] params) {
+ final double[][] jacobian = new double[points.size() * 2][3];
+
+ for (int i = 0; i < points.size(); i++) {
+ final int index = i * 2;
+ // Partial derivative wrt x-coordinate of center.
+ jacobian[index][0] = 1;
+ jacobian[index + 1][0] = 0;
+ // Partial derivative wrt y-coordinate of center.
+ jacobian[index][1] = 0;
+ jacobian[index + 1][1] = 1;
+ // Partial derivative wrt radius.
+ final double[] p = points.get(i);
+ jacobian[index][2] = (p[0] - params[0]) / params[2];
+ jacobian[index + 1][2] = (p[1] - params[1]) / params[2];
+ }
+
+ return jacobian;
+ }
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/CircleProblem.java
------------------------------------------------------------------------------
svn:eol-style = native
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java?rev=1330801&r1=1330800&r2=1330801&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java
Thu Apr 26 12:13:49 2012
@@ -590,6 +590,55 @@ public class LevenbergMarquardtOptimizer
}
}
+ @Test
+ public void testCircleFitting2() {
+ final double xCenter = 123.456;
+ final double yCenter = 654.321;
+ final double xSigma = 10;
+ final double ySigma = 15;
+ final double radius = 111.111;
+ final RandomCirclePointGenerator factory
+ = new RandomCirclePointGenerator(xCenter, yCenter, radius,
+ xSigma, ySigma,
+ 59421063L);
+ final CircleProblem circle = new CircleProblem(xSigma, ySigma);
+
+ final int numPoints = 10;
+ for (Point2D.Double p : factory.generate(numPoints)) {
+ circle.addPoint(p.x, p.y);
+ // System.out.println(p.x + " " + p.y);
+ }
+
+ // First guess for the center's coordinates and radius.
+ final double[] init = { 90, 659, 115 };
+
+ final LevenbergMarquardtOptimizer optimizer
+ = new LevenbergMarquardtOptimizer();
+ final PointVectorValuePair optimum = optimizer.optimize(100, circle,
+
circle.target(), circle.weight(),
+ init);
+
+ final double[] paramFound = optimum.getPoint();
+
+ // Retrieve errors estimation.
+ final double[][] covMatrix = optimizer.getCovariances();
+ final double[] asymptoticStandardErrorFound =
optimizer.guessParametersErrors();
+ final double[] sigmaFound = new double[covMatrix.length];
+ for (int i = 0; i < covMatrix.length; i++) {
+ sigmaFound[i] = FastMath.sqrt(covMatrix[i][i]);
+// System.out.println("i=" + i + " value=" + paramFound[i]
+// + " sigma=" + sigmaFound[i]
+// + " ase=" + asymptoticStandardErrorFound[i]);
+ }
+
+ // System.out.println("chi2=" + optimizer.getChiSquare());
+
+ // Check that the parameters are found within the assumed error bars.
+ Assert.assertEquals(xCenter, paramFound[0],
asymptoticStandardErrorFound[0]);
+ Assert.assertEquals(yCenter, paramFound[1],
asymptoticStandardErrorFound[1]);
+ Assert.assertEquals(radius, paramFound[2],
asymptoticStandardErrorFound[2]);
+ }
+
private static class LinearProblem implements
DifferentiableMultivariateVectorFunction, Serializable {
private static final long serialVersionUID = 703247177355019415L;
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java?rev=1330801&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
Thu Apr 26 12:13:49 2012
@@ -0,0 +1,95 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.optimization.general;
+
+import java.awt.geom.Point2D;
+import org.apache.commons.math3.random.RandomData;
+import org.apache.commons.math3.random.RandomDataImpl;
+import org.apache.commons.math3.random.Well44497b;
+import org.apache.commons.math3.util.MathUtils;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * Factory for generating a cloud of points that approximate a circle.
+ */
+public class RandomCirclePointGenerator {
+ /** RNG. */
+ private final RandomData random;
+ /** Radius of the circle. */
+ private final double radius;
+ /** x-coordinate of the circle center. */
+ private final double x;
+ /** y-coordinate of the circle center. */
+ private final double y;
+ /** Error on the x-coordinate of the center. */
+ private final double xSigma;
+ /** Error on the x-coordinate of the center. */
+ private final double ySigma;
+
+ /**
+ * @param x Abscissa of the circle center.
+ * @param y Ordinate of the circle center.
+ * @param radius Radius of the circle.
+ * @param xSigma Error on the x-coordinate of the circumferenc points.
+ * @param ySigma Error on the y-coordinate of the circumferenc points.
+ * @param seed RNG seed.
+ */
+ public RandomCirclePointGenerator(double x,
+ double y,
+ double radius,
+ double xSigma,
+ double ySigma,
+ long seed) {
+ random = new RandomDataImpl(new Well44497b((seed)));
+ this.radius = radius;
+ this.x = x;
+ this.y = y;
+ this.xSigma = xSigma;
+ this.ySigma = ySigma;
+ }
+
+ /**
+ * Point generator.
+ *
+ * @param n Number of points to create.
+ * @return the cloud of {@code n} points.
+ */
+ public Point2D.Double[] generate(int n) {
+ final Point2D.Double[] cloud = new Point2D.Double[n];
+ for (int i = 0; i < n; i++) {
+ cloud[i] = create();
+ }
+ return cloud;
+ }
+
+ /**
+ * Create one point.
+ *
+ * @return a point.
+ */
+ private Point2D.Double create() {
+ final double cX = random.nextGaussian(x, xSigma);
+ final double cY = random.nextGaussian(y, ySigma);
+ final double t = random.nextUniform(0, MathUtils.TWO_PI);
+
+ final double pX = cX + radius * FastMath.cos(t);
+ final double pY = cY + radius * FastMath.sin(t);
+
+ return new Point2D.Double(pX, pY);
+ }
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/RandomCirclePointGenerator.java
------------------------------------------------------------------------------
svn:eol-style = native
