1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.fitting.leastsquares;
18
19 import java.util.ArrayList;
20
21 import org.apache.commons.geometry.euclidean.twod.Vector2D;
22 import org.apache.commons.math4.legacy.analysis.MultivariateMatrixFunction;
23 import org.apache.commons.math4.legacy.analysis.MultivariateVectorFunction;
24
25
26
27
28 class CircleVectorial {
29 private ArrayList<Vector2D> points;
30
31 CircleVectorial() {
32 points = new ArrayList<>();
33 }
34
35 public void addPoint(double px, double py) {
36 points.add(Vector2D.of(px, py));
37 }
38
39 public int getN() {
40 return points.size();
41 }
42
43 public double getRadius(Vector2D center) {
44 double r = 0;
45 for (Vector2D point : points) {
46 r += point.distance(center);
47 }
48 return r / points.size();
49 }
50
51 public MultivariateVectorFunction getModelFunction() {
52 return new MultivariateVectorFunction() {
53 @Override
54 public double[] value(double[] params) {
55 Vector2D center = Vector2D.of(params[0], params[1]);
56 double radius = getRadius(center);
57 double[] residuals = new double[points.size()];
58 for (int i = 0; i < residuals.length; i++) {
59 residuals[i] = points.get(i).distance(center) - radius;
60 }
61
62 return residuals;
63 }
64 };
65 }
66
67 public MultivariateMatrixFunction getModelFunctionJacobian() {
68 return new MultivariateMatrixFunction() {
69 @Override
70 public double[][] value(double[] params) {
71 final int n = points.size();
72 final Vector2D center = Vector2D.of(params[0], params[1]);
73
74 double dRdX = 0;
75 double dRdY = 0;
76 for (Vector2D pk : points) {
77 double dk = pk.distance(center);
78 dRdX += (center.getX() - pk.getX()) / dk;
79 dRdY += (center.getY() - pk.getY()) / dk;
80 }
81 dRdX /= n;
82 dRdY /= n;
83
84
85 double[][] jacobian = new double[n][2];
86 for (int i = 0; i < n; i++) {
87 final Vector2D pi = points.get(i);
88 final double di = pi.distance(center);
89 jacobian[i][0] = (center.getX() - pi.getX()) / di - dRdX;
90 jacobian[i][1] = (center.getY() - pi.getY()) / di - dRdY;
91 }
92
93 return jacobian;
94 }
95 };
96 }
97 }