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.image.processing.convolution;
031
032import org.openimaj.image.FImage;
033import org.openimaj.image.analyser.ImageAnalyser;
034
035import net.jafama.FastMath;
036
037/**
038 * Image processor for calculating gradients and orientations using
039 * finite-differences. Both signed (+/- PI) orientations and unsigned (+/- PI/2)
040 * are computable.
041 *
042 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
043 *
044 */
045public class FImageGradients implements ImageAnalyser<FImage> {
046        /**
047         * Modes of operation for signed and unsigned orientations
048         *
049         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
050         *
051         */
052        public enum Mode {
053                /**
054                 * Unsigned orientations in +/- PI/2 computed using <code>atan</code>.
055                 */
056                Unsigned(-PI_OVER_TWO_FLOAT, PI_OVER_TWO_FLOAT) {
057                        @Override
058                        void gradientMagnitudesAndOrientations(FImage image, FImage magnitudes, FImage orientations) {
059                                FImageGradients.gradientMagnitudesAndUnsignedOrientations(image, magnitudes, orientations);
060                        }
061                },
062                /**
063                 * Signed orientations +/- PI computed using <code>atan2</code>.
064                 */
065                Signed(-PI_FLOAT, PI_FLOAT) {
066                        @Override
067                        void gradientMagnitudesAndOrientations(FImage image, FImage magnitudes, FImage orientations) {
068                                FImageGradients.gradientMagnitudesAndOrientations(image, magnitudes, orientations);
069                        }
070                };
071
072                private float min;
073                private float max;
074
075                private Mode(float min, float max) {
076                        this.min = min;
077                        this.max = max;
078                }
079
080                abstract void gradientMagnitudesAndOrientations(FImage image, FImage magnitudes, FImage orientations);
081
082                /**
083                 * Get the minimum angular value (in radians) computed by this mode.
084                 *
085                 * @return the minimum angular value.
086                 */
087                public float minAngle() {
088                        return min;
089                }
090
091                /**
092                 * Get the maximum angular value (in radians) computed by this mode.
093                 *
094                 * @return the maximum angular value.
095                 */
096                public float maxAngle() {
097                        return max;
098                }
099        }
100
101        private final static float PI_FLOAT = (float) Math.PI;
102        private final static float PI_OVER_TWO_FLOAT = (float) Math.PI / 2f;
103        private final static float TWO_PI_FLOAT = (float) (Math.PI * 2);
104
105        /**
106         * The gradient magnitudes
107         */
108        public FImage magnitudes;
109
110        /**
111         * The gradient orientations
112         */
113        public FImage orientations;
114
115        /**
116         * The orientation mode
117         */
118        public Mode mode;
119
120        /**
121         * Default constructor using {@link Mode#Signed} mode.
122         */
123        public FImageGradients() {
124                this.mode = Mode.Signed;
125        }
126
127        /**
128         * Construct using the given {@link Mode}.
129         *
130         * @param mode
131         *            the mode
132         */
133        public FImageGradients(Mode mode) {
134                this.mode = mode;
135        }
136
137        /*
138         * (non-Javadoc)
139         *
140         * @see
141         * org.openimaj.image.analyser.ImageAnalyser#analyseImage(org.openimaj.image
142         * .Image)
143         */
144        @Override
145        public void analyseImage(FImage image) {
146                if (magnitudes == null ||
147                                magnitudes.height != image.height ||
148                                magnitudes.width != image.width)
149                {
150                        magnitudes = new FImage(image.width, image.height);
151                        orientations = new FImage(image.width, image.height);
152                }
153
154                mode.gradientMagnitudesAndOrientations(image, magnitudes, orientations);
155        }
156
157        /**
158         * Static helper to create a new {@link FImageGradients} and call
159         * {@link FImageGradients#analyseImage(FImage)} with the image.
160         *
161         * @param image
162         *            the image
163         * @return a FImageGradients for the image
164         */
165        public static FImageGradients getGradientMagnitudesAndOrientations(FImage image) {
166                final FImageGradients go = new FImageGradients();
167                go.analyseImage(image);
168
169                return go;
170        }
171
172        /**
173         * Static helper to create a new {@link FImageGradients} and call
174         * {@link FImageGradients#analyseImage(FImage)} with the image.
175         *
176         * @param image
177         *            the image
178         * @param mode
179         *            the orientation mode
180         * @return a FImageGradients for the image
181         */
182        public static FImageGradients getGradientMagnitudesAndOrientations(FImage image, Mode mode) {
183                final FImageGradients go = new FImageGradients(mode);
184                go.analyseImage(image);
185
186                return go;
187        }
188
189        /**
190         * Estimate gradients magnitudes and orientations by calculating pixel
191         * differences. Edges get special treatment. The resultant gradients and
192         * orientations are returned though the gradients and orientations parameters
193         * respectively. The images represented by the gradients and orientations
194         * parameters are assumed to be initialized to the same size as the input image.
195         * Gradients are computed using the <code>atan2</code> function and will be in
196         * the range +/-PI.
197         *
198         * @param image
199         *            the input image
200         * @param magnitudes
201         *            the output gradient magnitudes
202         * @param orientations
203         *            the output gradient orientations
204         */
205        public static void gradientMagnitudesAndOrientations(FImage image, FImage magnitudes, FImage orientations) {
206                // Note: unrolling this loop to remove the if's doesn't
207                // actually seem to make it faster!
208                for (int r = 0; r < image.height; r++) {
209                        for (int c = 0; c < image.width; c++) {
210                                float xgrad, ygrad;
211
212                                if (c == 0)
213                                        xgrad = 2.0f * (image.pixels[r][c + 1] - image.pixels[r][c]);
214                                else if (c == image.width - 1)
215                                        xgrad = 2.0f * (image.pixels[r][c] - image.pixels[r][c - 1]);
216                                else
217                                        xgrad = image.pixels[r][c + 1] - image.pixels[r][c - 1];
218                                if (r == 0)
219                                        ygrad = 2.0f * (image.pixels[r][c] - image.pixels[r + 1][c]);
220                                else if (r == image.height - 1)
221                                        ygrad = 2.0f * (image.pixels[r - 1][c] - image.pixels[r][c]);
222                                else
223                                        ygrad = image.pixels[r - 1][c] - image.pixels[r + 1][c];
224
225                                // magnitudes.pixels[r][c] = (float) Math.sqrt( xgrad * xgrad +
226                                // ygrad * ygrad );
227                                // orientations.pixels[r][c] = (float) Math.atan2( ygrad, xgrad
228                                // );
229
230                                // JH - my benchmarking shows that (at least on OSX) Math.atan2
231                                // is really
232                                // slow... FastMath provides an alternative that is much faster
233                                magnitudes.pixels[r][c] = (float) Math.sqrt(xgrad * xgrad + ygrad * ygrad);
234                                orientations.pixels[r][c] = (float) FastMath.atan2(ygrad, xgrad);
235                        }
236                }
237        }
238
239        /**
240         * Estimate gradients magnitudes and orientations by calculating pixel
241         * differences. Edges get special treatment. The resultant gradients and
242         * orientations are returned though the gradients and orientations parameters
243         * respectively. The images represented by the gradients and orientations
244         * parameters are assumed to be initialized to the same size as the input image.
245         * Gradients are computed using the <code>atan</code> function and will be in
246         * the range +/- PI/2.
247         *
248         * @param image
249         *            the input image
250         * @param magnitudes
251         *            the output gradient magnitudes
252         * @param orientations
253         *            the output gradient orientations
254         */
255        public static void gradientMagnitudesAndUnsignedOrientations(FImage image, FImage magnitudes, FImage orientations) {
256                // Note: unrolling this loop to remove the if's doesn't
257                // actually seem to make it faster!
258                for (int r = 0; r < image.height; r++) {
259                        for (int c = 0; c < image.width; c++) {
260                                float xgrad, ygrad;
261
262                                if (c == 0)
263                                        xgrad = 2.0f * (image.pixels[r][c + 1] - image.pixels[r][c]);
264                                else if (c == image.width - 1)
265                                        xgrad = 2.0f * (image.pixels[r][c] - image.pixels[r][c - 1]);
266                                else
267                                        xgrad = image.pixels[r][c + 1] - image.pixels[r][c - 1];
268                                if (r == 0)
269                                        ygrad = 2.0f * (image.pixels[r][c] - image.pixels[r + 1][c]);
270                                else if (r == image.height - 1)
271                                        ygrad = 2.0f * (image.pixels[r - 1][c] - image.pixels[r][c]);
272                                else
273                                        ygrad = image.pixels[r - 1][c] - image.pixels[r + 1][c];
274
275                                magnitudes.pixels[r][c] = (float) Math.sqrt(xgrad * xgrad + ygrad * ygrad);
276                                if (magnitudes.pixels[r][c] == 0)
277                                        orientations.pixels[r][c] = 0;
278                                else
279                                        orientations.pixels[r][c] = (float) FastMath.atan(ygrad / xgrad);
280                        }
281                }
282        }
283
284        /**
285         * Estimate gradients magnitudes and orientations by calculating pixel
286         * differences. Edges get special treatment.
287         * <p>
288         * The orientations are quantised into <code>magnitudes.length</code> bins and
289         * the magnitudes are spread to the adjacent bin through linear interpolation.
290         * The magnitudes parameter must be fully allocated as an array of num
291         * orientation bin images, each of the same size as the input image.
292         *
293         * @param image
294         * @param magnitudes
295         */
296        public static void gradientMagnitudesAndQuantisedOrientations(FImage image, FImage[] magnitudes) {
297                final int numOriBins = magnitudes.length;
298
299                // Note: unrolling this loop to remove the if's doesn't
300                // actually seem to make it faster!
301                for (int r = 0; r < image.height; r++) {
302                        for (int c = 0; c < image.width; c++) {
303                                float xgrad, ygrad;
304
305                                if (c == 0)
306                                        xgrad = 2.0f * (image.pixels[r][c + 1] - image.pixels[r][c]);
307                                else if (c == image.width - 1)
308                                        xgrad = 2.0f * (image.pixels[r][c] - image.pixels[r][c - 1]);
309                                else
310                                        xgrad = image.pixels[r][c + 1] - image.pixels[r][c - 1];
311                                if (r == 0)
312                                        ygrad = 2.0f * (image.pixels[r][c] - image.pixels[r + 1][c]);
313                                else if (r == image.height - 1)
314                                        ygrad = 2.0f * (image.pixels[r - 1][c] - image.pixels[r][c]);
315                                else
316                                        ygrad = image.pixels[r - 1][c] - image.pixels[r + 1][c];
317
318                                // JH - my benchmarking shows that (at least on OSX) Math.atan2
319                                // is really
320                                // slow... FastMath provides an alternative that is much faster
321                                final float mag = (float) Math.sqrt(xgrad * xgrad + ygrad * ygrad);
322                                float ori = (float) FastMath.atan2(ygrad, xgrad);
323
324                                // adjust range
325                                ori = ((ori %= TWO_PI_FLOAT) >= 0 ? ori : (ori + TWO_PI_FLOAT));
326
327                                final float po = numOriBins * ori / TWO_PI_FLOAT; // po is now
328                                                                                                                                        // 0<=po<oriSize
329
330                                final int oi = (int) Math.floor(po);
331                                final float of = po - oi;
332
333                                // reset
334                                for (int i = 0; i < magnitudes.length; i++)
335                                        magnitudes[i].pixels[r][c] = 0;
336
337                                // set
338                                magnitudes[oi % numOriBins].pixels[r][c] = (1f - of) * mag;
339                                magnitudes[(oi + 1) % numOriBins].pixels[r][c] = of * mag;
340                        }
341                }
342        }
343
344        /**
345         * Estimate gradients magnitudes and orientations by calculating pixel
346         * differences. Edges get special treatment.
347         * <p>
348         * The orientations are quantised into <code>magnitudes.length</code> bins.
349         * Magnitudes are optionally spread to the adjacent bin through linear
350         * interpolation. The magnitudes parameter must be fully allocated as an array
351         * of num orientation bin images, each of the same size as the input image.
352         *
353         * @param image
354         * @param magnitudes
355         * @param interp
356         * @param mode
357         */
358        public static void gradientMagnitudesAndQuantisedOrientations(FImage image, FImage[] magnitudes, boolean interp,
359                        Mode mode)
360        {
361                final int numOriBins = magnitudes.length;
362
363                // Note: unrolling this loop to remove the if's doesn't
364                // actually seem to make it faster!
365                for (int r = 0; r < image.height; r++) {
366                        for (int c = 0; c < image.width; c++) {
367                                float xgrad, ygrad;
368
369                                if (c == 0)
370                                        xgrad = 2.0f * (image.pixels[r][c + 1] - image.pixels[r][c]);
371                                else if (c == image.width - 1)
372                                        xgrad = 2.0f * (image.pixels[r][c] - image.pixels[r][c - 1]);
373                                else
374                                        xgrad = image.pixels[r][c + 1] - image.pixels[r][c - 1];
375                                if (r == 0)
376                                        ygrad = 2.0f * (image.pixels[r][c] - image.pixels[r + 1][c]);
377                                else if (r == image.height - 1)
378                                        ygrad = 2.0f * (image.pixels[r - 1][c] - image.pixels[r][c]);
379                                else
380                                        ygrad = image.pixels[r - 1][c] - image.pixels[r + 1][c];
381
382                                // JH - my benchmarking shows that (at least on OSX) Math.atan2
383                                // is really
384                                // slow... FastMath provides an alternative that is much faster
385                                final float mag = (float) Math.sqrt(xgrad * xgrad + ygrad * ygrad);
386
387                                float po;
388                                if (mode == Mode.Unsigned) {
389                                        final float ori = mag == 0 ? PI_OVER_TWO_FLOAT
390                                                        : (float) FastMath.atan(ygrad / xgrad)
391                                                                        + PI_OVER_TWO_FLOAT;
392
393                                        po = numOriBins * ori / PI_FLOAT; // po is now 0<=po<oriSize
394                                } else {
395                                        float ori = (float) FastMath.atan2(ygrad, xgrad);
396
397                                        // adjust range
398                                        ori = ((ori %= TWO_PI_FLOAT) >= 0 ? ori : (ori + TWO_PI_FLOAT));
399
400                                        po = numOriBins * ori / TWO_PI_FLOAT; // po is now
401                                                                                                                        // 0<=po<oriSize
402                                }
403
404                                // reset
405                                for (int i = 0; i < magnitudes.length; i++)
406                                        magnitudes[i].pixels[r][c] = 0;
407
408                                int oi = (int) Math.floor(po);
409                                final float of = po - oi;
410
411                                // set
412                                if (interp) {
413                                        magnitudes[oi % numOriBins].pixels[r][c] = (1f - of) * mag;
414                                        magnitudes[(oi + 1) % numOriBins].pixels[r][c] = of * mag;
415                                } else {
416                                        // if (of > 0.5)
417                                        // magnitudes[(oi + 1) % numOriBins].pixels[r][c] = mag;
418                                        // else
419                                        if (oi > numOriBins - 1)
420                                                oi = numOriBins - 1;
421                                        magnitudes[oi].pixels[r][c] = mag;
422                                }
423                        }
424                }
425        }
426}