1 /**
2 * Copyright (c) 2011, The University of Southampton and the individual contributors.
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without modification,
6 * are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * * Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * * Neither the name of the University of Southampton nor the names of its
16 * contributors may be used to endorse or promote products derived from this
17 * software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
23 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
26 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30 package org.openimaj.image.analysis.watershed.feature;
31
32 import org.openimaj.image.FImage;
33 import org.openimaj.image.pixel.IntValuePixel;
34 import org.openimaj.math.geometry.point.Point2d;
35 import org.openimaj.math.geometry.shape.Circle;
36 import org.openimaj.math.geometry.shape.Ellipse;
37 import org.openimaj.math.geometry.shape.EllipseUtilities;
38 import org.openimaj.math.geometry.shape.Polygon;
39 import org.openimaj.math.util.QuadraticEquation;
40
41 import Jama.Matrix;
42
43 /**
44 * Accumulate the values of u11, u20 and u02 required to fit an
45 * ellipse to the feature.
46 *
47 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
48 *
49 */
50 public class MomentFeature implements ComponentFeature
51 {
52 int n = 0;
53 double mx = 0;
54 double my = 0;
55 double Mx2 = 0;
56 double My2 = 0;
57 double sxy = 0;
58 double sx = 0;
59 double sy = 0;
60
61 @Override
62 public void merge(ComponentFeature f) {
63 MomentFeature mf = (MomentFeature) f;
64
65 double dx = mf.mx - mx;
66 double dy = mf.my - my;
67
68 mx = (n*mx + mf.n*mf.mx) / (n + mf.n);
69 my = (n*my + mf.n*mf.my) / (n + mf.n);
70
71 Mx2 += mf.Mx2 + dx*dx*n*mf.n / (n + mf.n);
72 My2 += mf.My2 + dy*dy*n*mf.n / (n + mf.n);
73
74 n += mf.n;
75 sxy += mf.sxy;
76 sx += mf.sx;
77 sy += mf.sy;
78 }
79
80 @Override
81 public void addSample(IntValuePixel p) {
82 n++;
83 double dx = p.x - mx;
84 double dy = p.y - my;
85
86 mx += dx / n;
87 my += dy / n;
88
89 Mx2 += dx * (p.x - mx);
90 My2 += dy * (p.y - my);
91
92 sx += p.x;
93 sy += p.y;
94 sxy += p.x*p.y;
95 }
96
97 /**
98 * Get the number of pixels accumulated
99 * @return the number of pixels
100 */
101 public double n() {
102 return n;
103 }
104
105 /**
106 * Get the value of u11
107 * @return the value of u11
108 */
109 public double u11() {
110 return (sxy - sx*sy / n) / n;
111 }
112
113 /**
114 * Get the value of u20
115 * @return the value of u20
116 */
117 public double u20() {
118 return Mx2 / n;
119 }
120
121 /**
122 * Get the value of u02
123 * @return the value of u02
124 */
125 public double u02() {
126 return My2 / n;
127 }
128
129 /**
130 * Get the value of m10
131 * @return the value of m10
132 */
133 public double m10() {
134 return mx;
135 }
136
137 /**
138 * Get the value of m01
139 * @return the value of m01
140 */
141 public double m01() {
142 return my;
143 }
144
145 @Override
146 public MomentFeature clone() {
147 try {
148 return (MomentFeature) super.clone();
149 } catch (CloneNotSupportedException e) {
150 throw new AssertionError(e);
151 }
152 }
153
154 /**
155 * Create an ellipse based on the features parameters
156 * @return an ellipse
157 */
158 public Ellipse getEllipse() {
159 return getEllipse(1);
160 }
161 /**
162 * Create an ellipse based on the features parameters. Scale
163 * the ellipse about its centre with the given scale factor.
164 * @param sf the scale factor
165 * @return an ellipse
166 */
167 public Ellipse getEllipse(float sf) {
168 float u = (float) m10();
169 float v = (float) m01();
170 Matrix sm = new Matrix(new double[][]{
171 {u20(),u11()},
172 {u11(),u02()}
173 });
174 return EllipseUtilities.ellipseFromCovariance(u, v, sm, sf);
175 }
176
177 /**
178 * Create an ellipse based on the features parameters. Scale
179 * the ellipse about its centre with the given scale factor.
180 * @param sf the scale factor
181 * @return an ellipse
182 */
183 public Circle getCircle(float sf) {
184 Ellipse e = getEllipse(sf);
185 Point2d p = e.calculateCentroid();
186 return new Circle(p.getX(),p.getY(),(float)(e.getMajor() + e.getMinor())/2);
187 }
188
189 /**
190 * Create a rotated rectangle that fits around an ellipse
191 * based on the features parameters. Scale the ellipse about
192 * its centre with the given scale factor.
193 * @param sf the scale factor
194 * @return a polygon representing a rotated rectangle
195 */
196 public Polygon getEllipseBoundingBox(float sf) {
197 // double xx = u20();
198 // double xy = u11();
199 // double yy = u02();
200 //
201 // double theta = 0.5 * Math.atan2(2*xy, xx-yy);
202 //
203 // double trace = xx + yy;
204 // double det = (xx*yy) - (xy*xy);
205 // double [] eigval = QuadraticEquation.solveGeneralQuadratic(1, -1*trace, det);
206 //
207 // double a = Math.sqrt(eigval[1]) * sf * 2;
208 // double b = Math.sqrt(eigval[0]) * sf * 2;
209 //
210 // float Ax = (float) (a*Math.cos(theta) + b*Math.cos((Math.PI/2) + theta) + mx);
211 // float Ay = (float) (a*Math.sin(theta) + b*Math.sin((Math.PI/2) + theta) + my);
212 // float Bx = (float) (a*Math.cos(theta) - b*Math.cos((Math.PI/2) + theta) + mx);
213 // float By = (float) (a*Math.sin(theta) - b*Math.sin((Math.PI/2) + theta) + my);
214 // float Cx = (float) (-a*Math.cos(theta) - b*Math.cos((Math.PI/2) + theta) + mx);
215 // float Cy = (float) (-a*Math.sin(theta) - b*Math.sin((Math.PI/2) + theta) + my);
216 // float Dx = (float) (-a*Math.cos(theta) + b*Math.cos((Math.PI/2) + theta) + mx);
217 // float Dy = (float) (-a*Math.sin(theta) + b*Math.sin((Math.PI/2) + theta) + my);
218 //
219 // return new Polygon(new Point2dImpl(Ax,Ay), new Point2dImpl(Bx,By), new Point2dImpl(Cx,Cy), new Point2dImpl(Dx,Dy));
220 return this.getEllipse(sf).calculateOrientedBoundingBox();
221 }
222
223 /**
224 * Get the primary orientation of the feature.
225 * @return the orientation
226 */
227 public double getOrientation() {
228 double xx = u20();
229 double xy = u11();
230 double yy = u02();
231
232 return 0.5 * Math.atan2(2*xy, xx-yy);
233 }
234
235 protected double[] getEllipseBoundingRectsData(double sf) {
236 double xx = u20();
237 double xy = u11();
238 double yy = u02();
239
240 double theta = 0.5 * Math.atan2(2*xy, xx-yy);
241
242 double trace = xx + yy;
243 double det = (xx*yy) - (xy*xy);
244 double [] eigval = QuadraticEquation.solveGeneralQuadratic(1, -1*trace, det);
245
246 double a = 1.0 + Math.sqrt(Math.abs(eigval[1])) * sf * 4;
247 double b = 1.0 + Math.sqrt(Math.abs(eigval[0])) * sf * 4;
248
249 double [] data = {Math.max(a, b), Math.min(a, b), theta};
250
251 return data;
252 }
253
254 /**
255 * Extract a rectangular image patch centered on the feature
256 * with the same primary orientation and a given scale factor.
257 * @param image the image to extract from
258 * @param sf the scale factor
259 * @return a rectangular image patch
260 */
261 public FImage extractEllipsePatch(FImage image, double sf) {
262 double [] data = getEllipseBoundingRectsData(sf);
263 double height = data[1], width = data[0], ori = data[2];
264
265 int sx = (int) Math.rint(width);
266 int sy = (int) Math.rint(height);
267
268 FImage patch = new FImage(sx, sy);
269
270 //extract pixels
271 for (int y=0; y<sy; y++) {
272 for (int x=0; x<sx; x++) {
273 double xbar = x - sx / 2.0;
274 double ybar = y - sy / 2.0;
275
276 double xx = (xbar * Math.cos(ori) - ybar * Math.sin(ori)) + mx;
277 double yy = (xbar * Math.sin(ori) + ybar * Math.cos(ori)) + my;
278
279 patch.setPixel(x, y, image.getPixelInterp(xx, yy));
280 }
281 }
282
283 return patch;
284 }
285 }