001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.video.processing.motion;
031
032import org.apache.commons.math.complex.Complex;
033import org.apache.commons.math.linear.Array2DRowFieldMatrix;
034import org.apache.commons.math.linear.FieldLUDecompositionImpl;
035import org.openimaj.image.FImage;
036import org.openimaj.image.analysis.algorithm.TemplateMatcher;
037import org.openimaj.image.analysis.algorithm.TemplateMatcher.Mode;
038import org.openimaj.image.pixel.FValuePixel;
039import org.openimaj.image.processing.algorithm.FourierTransform;
040import org.openimaj.math.geometry.point.Point2d;
041import org.openimaj.math.geometry.point.Point2dImpl;
042import org.openimaj.math.geometry.shape.Rectangle;
043import org.openimaj.video.VideoFrame;
044import org.openimaj.video.VideoSubFrame;
045
046import edu.emory.mathcs.jtransforms.fft.FloatFFT_2D;
047
048/**
049 * A set of algorithms for the motion estimator.
050 *
051 * @author David Dupplaw (dpd@ecs.soton.ac.uk)
052 * @created 1 Mar 2012
053 *
054 */
055public abstract class MotionEstimatorAlgorithm
056{
057
058        /**
059         * Within a search window around the subimages detect most likely match and
060         * thus motion.
061         *
062         * @author Sina Samangooei (ss@ecs.soton.ac.uk)
063         *
064         */
065        public static class TEMPLATE_MATCH extends MotionEstimatorAlgorithm {
066                private float searchProp;
067                private Mode mode;
068
069                /**
070                 * Defaults to allowing a maximum of templatesize/2 movement using the
071                 * {@link Mode#CORRELATION}
072                 */
073                public TEMPLATE_MATCH() {
074                        this.searchProp = .5f;
075                        this.mode = TemplateMatcher.Mode.NORM_SUM_SQUARED_DIFFERENCE;
076                }
077
078                /**
079                 * Given the template's size, search around a border of size
080                 * template*searchWindowBorderProp
081                 *
082                 * @param searchWindowBorderProp
083                 * @param mode
084                 *            the matching mode
085                 */
086                public TEMPLATE_MATCH(float searchWindowBorderProp, Mode mode) {
087                        this.searchProp = searchWindowBorderProp;
088                        this.mode = mode;
089                }
090
091                @Override
092                Point2d estimateMotion(VideoSubFrame<FImage> img1sub,
093                                VideoSubFrame<FImage>... imagesSub)
094                {
095
096                        final VideoFrame<FImage> current = img1sub.extract();
097                        final VideoFrame<FImage> prev = imagesSub[0];
098                        final Rectangle prevSearchRect = imagesSub[0].roi;
099
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}