1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math3.optim.nonlinear.vector.jacobian; 18 19 import java.util.ArrayList; 20 21 import org.apache.commons.math3.analysis.MultivariateVectorFunction; 22 import org.apache.commons.math3.analysis.MultivariateMatrixFunction; 23 import org.apache.commons.math3.util.MathUtils; 24 import org.apache.commons.math3.util.FastMath; 25 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction; 26 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian; 27 28 /** 29 * Class that models a circle. 30 * The parameters of problem are: 31 * <ul> 32 * <li>the x-coordinate of the circle center,</li> 33 * <li>the y-coordinate of the circle center,</li> 34 * <li>the radius of the circle.</li> 35 * </ul> 36 * The model functions are: 37 * <ul> 38 * <li>for each triplet (cx, cy, r), the (x, y) coordinates of a point on the 39 * corresponding circle.</li> 40 * </ul> 41 */ 42 @Deprecated 43 class CircleProblem { 44 /** Cloud of points assumed to be fitted by a circle. */ 45 private final ArrayList<double[]> points; 46 /** Error on the x-coordinate of the points. */ 47 private final double xSigma; 48 /** Error on the y-coordinate of the points. */ 49 private final double ySigma; 50 /** Number of points on the circumference (when searching which 51 model point is closest to a given "observation". */ 52 private final int resolution; 53 54 /** 55 * @param xError Assumed error for the x-coordinate of the circle points. 56 * @param yError Assumed error for the y-coordinate of the circle points. 57 * @param searchResolution Number of points to try when searching the one 58 * that is closest to a given "observed" point. 59 */ CircleProblem(double xError, double yError, int searchResolution)60 public CircleProblem(double xError, 61 double yError, 62 int searchResolution) { 63 points = new ArrayList<double[]>(); 64 xSigma = xError; 65 ySigma = yError; 66 resolution = searchResolution; 67 } 68 69 /** 70 * @param xError Assumed error for the x-coordinate of the circle points. 71 * @param yError Assumed error for the y-coordinate of the circle points. 72 */ CircleProblem(double xError, double yError)73 public CircleProblem(double xError, 74 double yError) { 75 this(xError, yError, 500); 76 } 77 addPoint(double px, double py)78 public void addPoint(double px, double py) { 79 points.add(new double[] { px, py }); 80 } 81 target()82 public double[] target() { 83 final double[] t = new double[points.size() * 2]; 84 for (int i = 0; i < points.size(); i++) { 85 final double[] p = points.get(i); 86 final int index = i * 2; 87 t[index] = p[0]; 88 t[index + 1] = p[1]; 89 } 90 91 return t; 92 } 93 weight()94 public double[] weight() { 95 final double wX = 1 / (xSigma * xSigma); 96 final double wY = 1 / (ySigma * ySigma); 97 final double[] w = new double[points.size() * 2]; 98 for (int i = 0; i < points.size(); i++) { 99 final int index = i * 2; 100 w[index] = wX; 101 w[index + 1] = wY; 102 } 103 104 return w; 105 } 106 getModelFunction()107 public ModelFunction getModelFunction() { 108 return new ModelFunction(new MultivariateVectorFunction() { 109 public double[] value(double[] params) { 110 final double cx = params[0]; 111 final double cy = params[1]; 112 final double r = params[2]; 113 114 final double[] model = new double[points.size() * 2]; 115 116 final double deltaTheta = MathUtils.TWO_PI / resolution; 117 for (int i = 0; i < points.size(); i++) { 118 final double[] p = points.get(i); 119 final double px = p[0]; 120 final double py = p[1]; 121 122 double bestX = 0; 123 double bestY = 0; 124 double dMin = Double.POSITIVE_INFINITY; 125 126 // Find the angle for which the circle passes closest to the 127 // current point (using a resolution of 100 points along the 128 // circumference). 129 for (double theta = 0; theta <= MathUtils.TWO_PI; theta += deltaTheta) { 130 final double currentX = cx + r * FastMath.cos(theta); 131 final double currentY = cy + r * FastMath.sin(theta); 132 final double dX = currentX - px; 133 final double dY = currentY - py; 134 final double d = dX * dX + dY * dY; 135 if (d < dMin) { 136 dMin = d; 137 bestX = currentX; 138 bestY = currentY; 139 } 140 } 141 142 final int index = i * 2; 143 model[index] = bestX; 144 model[index + 1] = bestY; 145 } 146 147 return model; 148 } 149 }); 150 } 151 152 public ModelFunctionJacobian getModelFunctionJacobian() { 153 return new ModelFunctionJacobian(new MultivariateMatrixFunction() { 154 public double[][] value(double[] point) { 155 return jacobian(point); 156 } 157 }); 158 } 159 160 private double[][] jacobian(double[] params) { 161 final double[][] jacobian = new double[points.size() * 2][3]; 162 163 for (int i = 0; i < points.size(); i++) { 164 final int index = i * 2; 165 // Partial derivative wrt x-coordinate of center. 166 jacobian[index][0] = 1; 167 jacobian[index + 1][0] = 0; 168 // Partial derivative wrt y-coordinate of center. 169 jacobian[index][1] = 0; 170 jacobian[index + 1][1] = 1; 171 // Partial derivative wrt radius. 172 final double[] p = points.get(i); 173 jacobian[index][2] = (p[0] - params[0]) / params[2]; 174 jacobian[index + 1][2] = (p[1] - params[1]) / params[2]; 175 } 176 177 return jacobian; 178 } 179 } 180