Skip to content

Commit 1cc96ff

Browse files
gselzerctrueden
gselzer
authored and
ctrueden
committed
Add Quadric Op/Test
1 parent 67dc938 commit 1cc96ff

File tree

2 files changed

+281
-0
lines changed

2 files changed

+281
-0
lines changed

‎src/main/java/net/imagej/ops/stats/regression/leastSquares/Quadric.java‎

Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,166 @@
1+
/*-
2+
* #%L
3+
* ImageJ software for multidimensional image processing and analysis.
4+
* %%
5+
* Copyright (C) 2014 - 2018 ImageJ developers.
6+
* %%
7+
* Redistribution and use in source and binary forms, with or without
8+
* modification, are permitted provided that the following conditions are met:
9+
*
10+
* 1. Redistributions of source code must retain the above copyright notice,
11+
* this list of conditions and the following disclaimer.
12+
* 2. Redistributions in binary form must reproduce the above copyright notice,
13+
* this list of conditions and the following disclaimer in the documentation
14+
* and/or other materials provided with the distribution.
15+
*
16+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
20+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26+
* POSSIBILITY OF SUCH DAMAGE.
27+
* #L%
28+
*/
29+
30+
package net.imagej.ops.stats.regression.leastSquares;
31+
32+
import java.util.Collection;
33+
import java.util.Iterator;
34+
import java.util.function.Function;
35+
36+
import org.joml.Matrix4d;
37+
import org.joml.Vector3d;
38+
import org.ojalgo.matrix.BasicMatrix;
39+
import org.ojalgo.matrix.PrimitiveMatrix;
40+
import org.ojalgo.random.Deterministic;
41+
import org.scijava.ops.core.Op;
42+
import org.scijava.param.Parameter;
43+
import org.scijava.plugin.Plugin;
44+
import org.scijava.struct.ItemIO;
45+
46+
/**
47+
* An op that fits a quadratic surface (quadric) into a set of points.
48+
* <p>
49+
* The op first solves the quadric that best fits the point cloud by minimising
50+
* the distance by least squares fitting. It's found by solving a polynomial -
51+
* the general equation of a quadric. There are no guarantees about the type of
52+
* the quadric solved, and it can be real or imaginary. The method is sensitive
53+
* to outlier points.
54+
* </p>
55+
* <p>
56+
* The op is based on the the implementations of Yury Petrov &amp; KalebKE.
57+
* </p>
58+
*
59+
* @author Richard Domander (Royal Veterinary College, London)
60+
*/
61+
@Plugin(type = Op.class, name = "stats.leastSquares")
62+
@Parameter(key = "vectorCollection")
63+
@Parameter(key = "outputMatrix", type = ItemIO.OUTPUT)
64+
public class Quadric implements
65+
Function<Collection<Vector3d>, Matrix4d>
66+
{
67+
68+
/**
69+
* Minimum number of points in the input collection needed to fit a quadric
70+
* equation.
71+
*/
72+
public static final int MIN_DATA = 9;
73+
74+
@Override
75+
public Matrix4d apply(final Collection<Vector3d> points) {
76+
if (points.size() < MIN_DATA)
77+
throw new IllegalArgumentException("Must pass more points in order to fit a quadric equation!");
78+
final double[] vector = solveVector(points);
79+
return toQuadricMatrix(vector);
80+
}
81+
82+
/**
83+
* Creates a design matrix used for least squares fitting from a collection of
84+
* points.
85+
*
86+
* @see #solveVector(Collection)
87+
* @param points points in a 3D space.
88+
* @return a [points.size()][9] matrix of real values.
89+
*/
90+
private static PrimitiveMatrix createDesignMatrix(
91+
final Collection<Vector3d> points)
92+
{
93+
final BasicMatrix.Builder<PrimitiveMatrix> builder = PrimitiveMatrix.FACTORY
94+
.getBuilder(points.size(), MIN_DATA);
95+
final Iterator<Vector3d> iterator = points.iterator();
96+
for (int i = 0; i < points.size(); i++) {
97+
final Vector3d p = iterator.next();
98+
builder.set(i, 0, p.x * p.x);
99+
builder.set(i, 1, p.y * p.y);
100+
builder.set(i, 2, p.z * p.z);
101+
builder.set(i, 3, 2 * p.x * p.y);
102+
builder.set(i, 4, 2 * p.x * p.z);
103+
builder.set(i, 5, 2 * p.y * p.z);
104+
builder.set(i, 6, 2 * p.x);
105+
builder.set(i, 7, 2 * p.y);
106+
builder.set(i, 8, 2 * p.z);
107+
}
108+
return builder.build();
109+
}
110+
111+
/**
112+
* Solves the equation for the quadratic surface that best fits the given
113+
* points.
114+
* <p>
115+
* The vector solved is the polynomial Ax<sup>2</sup> + By<sup>2</sup> +
116+
* Cz<sup>2</sup> + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz, i.e. the general
117+
* equation of a quadric. The fitting is done with least squares.
118+
* </p>
119+
*
120+
* @param points A collection of points in a 3D space.
121+
* @return the solution vector of the surface.
122+
*/
123+
private static double[] solveVector(final Collection<Vector3d> points) {
124+
final int n = points.size();
125+
// Find (dT * d)^-1
126+
final PrimitiveMatrix d = createDesignMatrix(points);
127+
final PrimitiveMatrix dT = d.transpose();
128+
final PrimitiveMatrix dTDInv = dT.multiply(d).invert();
129+
// Multiply dT * O, where O = [1, 1, ... 1] (n x 1) matrix
130+
final PrimitiveMatrix o = PrimitiveMatrix.FACTORY.makeFilled(n, 1,
131+
new Deterministic(1.0));
132+
final PrimitiveMatrix dTO = dT.multiply(o);
133+
// Find solution A = (dT * d)^-1 * (dT * O)
134+
return dTDInv.multiply(dTO).toRawCopy1D();
135+
}
136+
137+
/**
138+
* Creates a matrix out of a quadric surface solution vector in homogeneous
139+
* coordinates.
140+
*
141+
* @see #solveVector(Collection)
142+
* @return a matrix representing the polynomial solution vector in an
143+
* algebraic form.
144+
*/
145+
private Matrix4d toQuadricMatrix(final double[] solution) {
146+
// I'm not a clever man, so I'm using local variables to
147+
// better follow the matrix assignment.
148+
final double a = solution[0];
149+
final double b = solution[1];
150+
final double c = solution[2];
151+
final double d = solution[3];
152+
final double e = solution[4];
153+
final double f = solution[5];
154+
final double g = solution[6];
155+
final double h = solution[7];
156+
final double i = solution[8];
157+
// @formatter:off
158+
return new Matrix4d(
159+
a, d, e, g,
160+
d, b, f, h,
161+
e, f, c, i,
162+
g, h, i, -1
163+
);
164+
// @formatter:on
165+
}
166+
}

‎src/test/java/net/imagej/ops/stats/regression/leastSquares/QuadricTest.java‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
/*-
2+
* #%L
3+
* ImageJ software for multidimensional image processing and analysis.
4+
* %%
5+
* Copyright (C) 2014 - 2018 ImageJ developers.
6+
* %%
7+
* Redistribution and use in source and binary forms, with or without
8+
* modification, are permitted provided that the following conditions are met:
9+
*
10+
* 1. Redistributions of source code must retain the above copyright notice,
11+
* this list of conditions and the following disclaimer.
12+
* 2. Redistributions in binary form must reproduce the above copyright notice,
13+
* this list of conditions and the following disclaimer in the documentation
14+
* and/or other materials provided with the distribution.
15+
*
16+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
20+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26+
* POSSIBILITY OF SUCH DAMAGE.
27+
* #L%
28+
*/
29+
30+
package net.imagej.ops.stats.regression.leastSquares;
31+
32+
import static java.util.stream.Collectors.toList;
33+
import static org.junit.Assert.assertEquals;
34+
35+
import java.util.List;
36+
import java.util.stream.Stream;
37+
38+
import net.imagej.ops.AbstractOpTest;
39+
40+
import org.joml.Matrix4d;
41+
import org.joml.Matrix4dc;
42+
import org.joml.Vector3d;
43+
import org.junit.Test;
44+
45+
/**
46+
* Tests for {@link Quadric}.
47+
*
48+
* @author Richard Domander (Royal Veterinary College, London)
49+
*/
50+
public class QuadricTest extends AbstractOpTest {
51+
52+
private static final double alpha = Math.cos(Math.PI / 4.0);
53+
private static final List<Vector3d> unitSpherePoints = Stream.of(new Vector3d(
54+
1, 0, 0), new Vector3d(-1, 0, 0), new Vector3d(0, 1, 0), new Vector3d(0, -1,
55+
0), new Vector3d(0, 0, 1), new Vector3d(0, 0, -1), new Vector3d(alpha,
56+
alpha, 0), new Vector3d(-alpha, alpha, 0), new Vector3d(alpha, -alpha,
57+
0), new Vector3d(-alpha, -alpha, 0), new Vector3d(0, alpha, alpha),
58+
new Vector3d(0, -alpha, alpha), new Vector3d(0, alpha, -alpha),
59+
new Vector3d(0, -alpha, -alpha), new Vector3d(alpha, 0, alpha),
60+
new Vector3d(alpha, 0, -alpha), new Vector3d(-alpha, 0, alpha),
61+
new Vector3d(-alpha, 0, -alpha)).collect(toList());
62+
63+
@Test
64+
public void testEquation() {
65+
final Matrix4dc solution = (Matrix4dc) ops.run("stats.leastSquares",
66+
unitSpherePoints);
67+
final double a = solution.m00();
68+
final double b = solution.m11();
69+
final double c = solution.m22();
70+
final double d = solution.m01();
71+
final double e = solution.m02();
72+
final double f = solution.m12();
73+
final double g = solution.m03();
74+
final double h = solution.m13();
75+
final double i = solution.m23();
76+
77+
for (final Vector3d p : unitSpherePoints) {
78+
final double polynomial = a * p.x * p.x + b * p.y * p.y + c * p.z * p.z +
79+
2 * d * p.x * p.y + 2 * e * p.x * p.z + 2 * f * p.y * p.z + 2 * g *
80+
p.x + 2 * h * p.y + 2 * i * p.z;
81+
assertEquals("The matrix does not solve the polynomial equation", 1.0,
82+
polynomial, 1e-12);
83+
}
84+
}
85+
86+
@Test(expected = IllegalArgumentException.class)
87+
public void testMatchingFailsIfTooFewPoints() {
88+
final int nPoints = Math.max(0, Quadric.MIN_DATA - 1);
89+
final List<Vector3d> points = Stream.generate(Vector3d::new).limit(nPoints)
90+
.collect(toList());
91+
92+
ops.run("stats.leastSquares", points);
93+
}
94+
95+
@Test
96+
public void testMatrixElements() {
97+
final Matrix4dc solution = (Matrix4dc) ops.run("stats.leastSquares",
98+
unitSpherePoints);
99+
100+
assertEquals("The matrix element is incorrect", 1.0, solution.m00(), 1e-12);
101+
assertEquals("The matrix element is incorrect", 1.0, solution.m11(), 1e-12);
102+
assertEquals("The matrix element is incorrect", 1.0, solution.m22(), 1e-12);
103+
assertEquals("The matrix element is incorrect", 0.0, solution.m01(), 1e-12);
104+
assertEquals("The matrix element is incorrect", 0.0, solution.m02(), 1e-12);
105+
assertEquals("The matrix element is incorrect", 0.0, solution.m03(), 1e-12);
106+
assertEquals("The matrix element is incorrect", 0.0, solution.m12(), 1e-12);
107+
assertEquals("The matrix element is incorrect", 0.0, solution.m13(), 1e-12);
108+
assertEquals("The matrix element is incorrect", 0.0, solution.m23(), 1e-12);
109+
assertEquals("The matrix element is incorrect", -1.0, solution.m33(),
110+
1e-12);
111+
final Matrix4d transposed = new Matrix4d();
112+
solution.transpose(transposed);
113+
assertEquals("Matrix is not symmetric", solution, transposed);
114+
}
115+
}

0 commit comments

Comments
 (0)