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.video.processing.motion;
31
32 import org.apache.commons.math.complex.Complex;
33 import org.apache.commons.math.linear.Array2DRowFieldMatrix;
34 import org.apache.commons.math.linear.FieldLUDecompositionImpl;
35 import org.openimaj.image.FImage;
36 import org.openimaj.image.analysis.algorithm.TemplateMatcher;
37 import org.openimaj.image.analysis.algorithm.TemplateMatcher.Mode;
38 import org.openimaj.image.pixel.FValuePixel;
39 import org.openimaj.image.processing.algorithm.FourierTransform;
40 import org.openimaj.math.geometry.point.Point2d;
41 import org.openimaj.math.geometry.point.Point2dImpl;
42 import org.openimaj.math.geometry.shape.Rectangle;
43 import org.openimaj.video.VideoFrame;
44 import org.openimaj.video.VideoSubFrame;
45
46 import edu.emory.mathcs.jtransforms.fft.FloatFFT_2D;
47
48 /**
49 * A set of algorithms for the motion estimator.
50 *
51 * @author David Dupplaw (dpd@ecs.soton.ac.uk)
52 * @created 1 Mar 2012
53 *
54 */
55 public abstract class MotionEstimatorAlgorithm
56 {
57
58 /**
59 * Within a search window around the subimages detect most likely match and
60 * thus motion.
61 *
62 * @author Sina Samangooei (ss@ecs.soton.ac.uk)
63 *
64 */
65 public static class TEMPLATE_MATCH extends MotionEstimatorAlgorithm {
66 private float searchProp;
67 private Mode mode;
68
69 /**
70 * Defaults to allowing a maximum of templatesize/2 movement using the
71 * {@link Mode#CORRELATION}
72 */
73 public TEMPLATE_MATCH() {
74 this.searchProp = .5f;
75 this.mode = TemplateMatcher.Mode.NORM_SUM_SQUARED_DIFFERENCE;
76 }
77
78 /**
79 * Given the template's size, search around a border of size
80 * template*searchWindowBorderProp
81 *
82 * @param searchWindowBorderProp
83 * @param mode
84 * the matching mode
85 */
86 public TEMPLATE_MATCH(float searchWindowBorderProp, Mode mode) {
87 this.searchProp = searchWindowBorderProp;
88 this.mode = mode;
89 }
90
91 @Override
92 Point2d estimateMotion(VideoSubFrame<FImage> img1sub,
93 VideoSubFrame<FImage>... imagesSub)
94 {
95
96 final VideoFrame<FImage> current = img1sub.extract();
97 final VideoFrame<FImage> prev = imagesSub[0];
98 final Rectangle prevSearchRect = imagesSub[0].roi;
99
100 final int sw = (int) (prevSearchRect.width * this.searchProp);
101 final int sh = (int) (prevSearchRect.height * this.searchProp);
102 final int sx = (int) (prevSearchRect.x + ((prevSearchRect.width - sw) / 2.f));
103 final int sy = (int) (prevSearchRect.y + ((prevSearchRect.height - sh) / 2.f));
104
105 final Rectangle searchRect = new Rectangle(sx, sy, sw, sh);
106 // System.out.println("Search window: " + searchRect);
107 // MBFImage searchRectDraw = new
108 // MBFImage(img1sub.frame.clone(),img1sub.frame.clone(),img1sub.frame.clone());
109 // searchRectDraw.drawShape(searchRect, RGBColour.RED);
110 // searchRectDraw.drawPoint(img1sub.roi.getCOG(), RGBColour.GREEN,
111 // 3);
112 final TemplateMatcher matcher = new TemplateMatcher(current.frame, mode, searchRect);
113 matcher.analyseImage(prev.frame);
114 final FValuePixel[] responses = matcher.getBestResponses(1);
115 final FValuePixel firstBest = responses[0];
116 // for (FValuePixel bestRespose : responses) {
117 // if(bestRespose == null) continue;
118 // if(firstBest == null) firstBest = bestRespose;
119 // bestRespose.translate(current.frame.width/2,
120 // current.frame.height/2);
121 // // searchRectDraw.drawPoint(bestRespose, RGBColour.BLUE, 3);
122 // }
123 final Point2d centerOfGrid = img1sub.roi.calculateCentroid();
124 // System.out.println("First reponse: " + firstBest );
125 // System.out.println("Center of template: " + centerOfGrid);
126
127 // DisplayUtilities.displayName(searchRectDraw, "searchWindow");
128 if (firstBest == null || Float.isNaN(firstBest.value))
129 return new Point2dImpl(0, 0);
130 // firstBest.translate(current.frame.width/2,
131 // current.frame.height/2);
132 // System.out.println("First reponse (corrected): " + firstBest );
133 // System.out.println("Diff: " + centerOfGrid.minus(firstBest));
134 return centerOfGrid.minus(firstBest);
135 }
136 }
137
138 /**
139 * Basic phase correlation algorithm that finds peaks in the cross-power
140 * spectrum between two images. This is the basic implementation without
141 * sub-pixel accuracy.
142 */
143 public static class PHASE_CORRELATION extends MotionEstimatorAlgorithm
144 {
145 /**
146 * Calculate the estimated motion vector between <code>images</code>
147 * which [0] is first in the sequence and <code>img2</code> which is
148 * second in the sequence. This method uses phase correlation - the fact
149 * that translations in space can be seen as shifts in phase in the
150 * frequency domain. The returned vector will have a maximum horizontal
151 * displacement of <code>img2.getWidth()/2</code> and a minimum
152 * displacement of <code>-img2.getWidth()/2</code> and similarly for the
153 * vertical displacement and height.
154 *
155 * @param img2sub
156 * The second image in the sequence
157 * @param imagesSub
158 * The previous image in the sequence
159 * @return the estimated motion vector as a {@link Point2d} in absolute
160 * x and y coordinates.
161 */
162 @Override
163 public Point2d estimateMotion(VideoSubFrame<FImage> img2sub,
164 VideoSubFrame<FImage>... imagesSub)
165 {
166 // The previous image will be the first in the images array
167 final FImage img1 = imagesSub[0].extract().frame;
168 final VideoFrame<FImage> img2 = img2sub.extract();
169
170 // No previous frame?
171 if (img1 == null)
172 return new Point2dImpl(0, 0);
173
174 // The images must have comparable shapes and must be square
175 if (img1.getRows() != img2.frame.getRows() ||
176 img1.getCols() != img2.frame.getCols() ||
177 img1.getCols() != img2.frame.getRows())
178 return new Point2dImpl(0, 0);
179
180 // Prepare and perform an FFT for each of the incoming images.
181 final int h = img1.getRows();
182 final int w = img1.getCols();
183
184 try
185 {
186 final FloatFFT_2D fft1 = new FloatFFT_2D(h, w);
187 final FloatFFT_2D fft2 = new FloatFFT_2D(h, w);
188 final float[][] data1 = FourierTransform.prepareData(img1, h, w, false);
189 final float[][] data2 = FourierTransform.prepareData(img2.frame, h, w, false);
190 fft1.complexForward(data1);
191 fft2.complexForward(data2);
192
193 // Multiply (element-wise) the fft and the conjugate of the fft.
194 Complex[][] cfft = new Complex[h][w];
195 for (int y = 0; y < h; y++)
196 {
197 for (int x = 0; x < w; x++)
198 {
199 final float re1 = data1[y][x * 2];
200 final float im1 = data1[y][1 + x * 2];
201 final float re2 = data2[y][x * 2];
202 final float im2 = data2[y][1 + x * 2];
203
204 final Complex c1 = new Complex(re1, im1);
205 final Complex c2 = new Complex(re2, -im2);
206 cfft[y][x] = c1.multiply(c2);
207 }
208 }
209
210 // ----------------------------------------
211 // Normalise by the determinant
212 // ----------------------------------------
213 // First calculate the determinant
214 final Array2DRowFieldMatrix<Complex> cmat =
215 new Array2DRowFieldMatrix<Complex>(cfft);
216 final FieldLUDecompositionImpl<Complex> luDecomp =
217 new FieldLUDecompositionImpl<Complex>(cmat);
218 final Complex det = luDecomp.getDeterminant();
219 cmat.scalarMultiply(new Complex(1d / det.abs(), 0));
220
221 // Convert back to an array for doing the inverse FFTing
222 cfft = cmat.getData();
223 for (int y = 0; y < h; y++)
224 {
225 for (int x = 0; x < w; x++)
226 {
227 data1[y][x * 2] = (float) cfft[y][x].getReal();
228 data1[y][1 + x * 2] = (float) cfft[y][x].getImaginary();
229 }
230 }
231
232 // Perform the inverse FFT
233 fft1.complexInverse(data1, false);
234
235 // Get the data back out
236 FourierTransform.unprepareData(data1, img1, false);
237
238 // Get the estimated motion vector from the peak in the space
239 final FValuePixel p = img1.maxPixel();
240 return new Point2dImpl(
241 -(p.x > w / 2 ? p.x - w : p.x),
242 -(p.y > w / 2 ? p.y - w : p.y));
243 } catch (final Exception e)
244 {
245 return new Point2dImpl(0, 0);
246 }
247 }
248 };
249
250 /**
251 * Estimate the motion to the given subimage, <code>img1sub</code> from the
252 * previous frames. The previous frames will be given in reverse order so
253 * that imagesSub[0] will be the previous frame, imagesSub[1] the frame
254 * before that, etc. The number of frames given will be at most that given
255 * by {@link #requiredNumberOfFrames()}. It could be less if at the
256 * beginning of the video. If you require more frames, return an empty
257 * motion vector - that is (0,0).
258 *
259 * @param img1sub
260 * The image to which we want to estimate the motion.
261 * @param imagesSub
262 * The previous frames in reverse order
263 * @return The estimated motion vector.
264 */
265 abstract Point2d estimateMotion(VideoSubFrame<FImage> img1sub,
266 VideoSubFrame<FImage>... imagesSub);
267
268 /**
269 * The required number of frames required for the given motion estimation
270 * algorithm to work. The default is 1 which means the algorithm will only
271 * receive the previous frame. If more are required, override this method
272 * and return the required number.
273 *
274 * @return The required number of frames to pass to the algorithm.
275 */
276 protected int requiredNumberOfFrames()
277 {
278 return 1;
279 }
280 }