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.algorithm;
031
032import org.openimaj.citation.annotation.Reference;
033import org.openimaj.citation.annotation.ReferenceType;
034import org.openimaj.image.FImage;
035import org.openimaj.image.pixel.ConnectedComponent;
036import org.openimaj.image.pixel.ConnectedComponent.ConnectMode;
037import org.openimaj.image.processor.SinglebandImageProcessor;
038
039import net.jafama.FastMath;
040
041/**
042 * Implementation of Perona & Malik's image filtering by anisotropic
043 * diffusion. Enables edge-preserving image smoothing.
044 *
045 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
046 */
047@Reference(
048                type = ReferenceType.Article,
049                author = { "Perona, P.", "Malik, J." },
050                title = "Scale-Space and Edge Detection Using Anisotropic Diffusion",
051                year = "1990",
052                journal = "IEEE Trans. Pattern Anal. Mach. Intell.",
053                pages = { "629", "", "639" },
054                url = "http://dx.doi.org/10.1109/34.56205",
055                month = "July",
056                number = "7",
057                publisher = "IEEE Computer Society",
058                volume = "12",
059                customData = {
060                                "issue_date", "July 1990",
061                                "issn", "0162-8828",
062                                "numpages", "11",
063                                "doi", "10.1109/34.56205",
064                                "acmid", "78304",
065                                "address", "Washington, DC, USA"
066                })
067public class AnisotropicDiffusion implements SinglebandImageProcessor<Float, FImage> {
068        /**
069         * Interface describing a function that computes the conduction coefficient as a
070         * function of space and gradient magnitude.
071         *
072         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
073         *
074         */
075        public static interface ConductionCoefficientFunction {
076                /**
077                 * Compute the conduction coefficient given gradient and position.
078                 *
079                 * @param grad
080                 *            the gradient
081                 * @param x
082                 *            the x-position in the image
083                 * @param y
084                 *            the y-position in the image
085                 * @return the conduction coefficient
086                 */
087                float apply(float grad, int x, int y);
088        }
089
090        /**
091         * The first conduction function proposed by Perona &amp; Malik, that privileges
092         * high contrast edges over low constrast ones.
093         *
094         * <pre>
095         * g(∇I) = exp( - (||∇I/κ||²) )
096         * </pre>
097         *
098         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
099         */
100        public static class HighConstrastEdgeFunction implements ConductionCoefficientFunction {
101                private float kappa;
102
103                /**
104                 * Construct with the given kappa value
105                 *
106                 * @param kappa
107                 *            kappa
108                 */
109                public HighConstrastEdgeFunction(float kappa) {
110                        this.kappa = kappa;
111                }
112
113                @Override
114                public float apply(float val, int x, int y) {
115                        final float t = val / kappa;
116                        return (float) FastMath.exp(-t * t);
117                }
118        };
119
120        /**
121         * The second conduction function proposed by Perona &amp; Malik, that
122         * privileges wider regions over smaller ones.
123         *
124         * <pre>
125         * g(∇I) = 1 / ( 1 + (||∇I/κ||²) )
126         * </pre>
127         *
128         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
129         */
130        public static class WideRegionFunction implements ConductionCoefficientFunction {
131                private float kappa;
132
133                /**
134                 * Construct with the given kappa value
135                 *
136                 * @param kappa
137                 *            kappa
138                 */
139                public WideRegionFunction(float kappa) {
140                        this.kappa = kappa;
141                }
142
143                @Override
144                public float apply(float val, int x, int y) {
145                        final float t = val / kappa;
146                        return 1 / (1 + t * t);
147                }
148        };
149
150        protected int numIterations;
151        protected float lambda = 1f / 7f;
152        protected ConductionCoefficientFunction function;
153        protected ConnectedComponent.ConnectMode neighbourMode;
154
155        /**
156         * Construct using the given parameters.
157         *
158         * @param numIterations
159         *            number of iterations of filtering to apply
160         * @param lambda
161         *            the integration constant (less than or equal to 1/4 for
162         *            4-neighbour mode and 1/8 for 8-neighbour)
163         * @param function
164         *            the conduction coefficient function
165         * @param neighbourMode
166         *            the neighbourhood mode
167         */
168        public AnisotropicDiffusion(int numIterations, float lambda, ConductionCoefficientFunction function,
169                        ConnectMode neighbourMode)
170        {
171                this.numIterations = numIterations;
172                this.lambda = lambda;
173                this.function = function;
174                this.neighbourMode = neighbourMode;
175        }
176
177        /**
178         * Construct using the given parameters. 4-neighbour mode is used as per the
179         * original paper.
180         *
181         * @param numIterations
182         *            number of iterations of filtering to apply
183         * @param lambda
184         *            the integration constant (bigger than 0 and less than or equal to
185         *            1/4)
186         * @param function
187         *            the conduction coefficient function
188         */
189        public AnisotropicDiffusion(int numIterations, float lambda, ConductionCoefficientFunction function) {
190                this.numIterations = numIterations;
191                this.lambda = lambda;
192                this.function = function;
193                this.neighbourMode = ConnectMode.CONNECT_4;
194        }
195
196        @Override
197        public void processImage(FImage image) {
198                for (int i = 0; i < numIterations; i++) {
199                        processImageOneIteration(image);
200                }
201        }
202
203        /**
204         * Perform a single iteration of anisotropic diffusion
205         *
206         * @param image
207         *            the image
208         */
209        public void processImageOneIteration(FImage image) {
210                switch (neighbourMode) {
211                case CONNECT_4:
212                        processImageOneIteration4(image);
213                case CONNECT_8:
214                        processImageOneIteration8(image);
215                }
216        }
217
218        private void processImageOneIteration4(FImage image) {
219                final float[][] tmp = image.clone().pixels;
220
221                for (int y = 0; y < image.height; y++) {
222                        final int ym = Math.max(y - 1, 0);
223                        final int yp = Math.min(y + 1, image.height - 1);
224
225                        for (int x = 0; x < image.width; x++) {
226                                final int xm = Math.max(x - 1, 0);
227                                final int xp = Math.min(x + 1, image.width - 1);
228
229                                final float dN = tmp[ym][x] - tmp[y][x];
230                                final float dS = tmp[yp][x] - tmp[y][x];
231                                final float dE = tmp[y][xm] - tmp[y][x];
232                                final float dW = tmp[y][xp] - tmp[y][x];
233
234                                final float cN = function.apply(Math.abs(dN), x, y);
235                                final float cS = function.apply(Math.abs(dS), x, y);
236                                final float cE = function.apply(Math.abs(dE), x, y);
237                                final float cW = function.apply(Math.abs(dW), x, y);
238
239                                image.pixels[y][x] += lambda * (cN * dN + cS * dS + cE * dE + cW * dW);
240                        }
241                }
242        }
243
244        private void processImageOneIteration8(FImage image) {
245                final float[][] tmp = image.clone().pixels;
246                final float wt = 0.5f;
247
248                for (int y = 0; y < image.height; y++) {
249                        final int ym = Math.max(y - 1, 0);
250                        final int yp = Math.min(y + 1, image.height - 1);
251
252                        for (int x = 0; x < image.width; x++) {
253                                final int xm = Math.max(x - 1, 0);
254                                final int xp = Math.min(x + 1, image.width - 1);
255
256                                final float dN = tmp[ym][x] - tmp[y][x];
257                                final float dS = tmp[yp][x] - tmp[y][x];
258                                final float dE = tmp[y][xm] - tmp[y][x];
259                                final float dW = tmp[y][xp] - tmp[y][x];
260                                final float dNE = tmp[ym][xm] - tmp[y][x];
261                                final float dSE = tmp[yp][xm] - tmp[y][x];
262                                final float dSW = tmp[ym][xp] - tmp[y][x];
263                                final float dNW = tmp[ym][xp] - tmp[y][x];
264
265                                final float cN = function.apply(Math.abs(dN), x, y);
266                                final float cS = function.apply(Math.abs(dS), x, y);
267                                final float cE = function.apply(Math.abs(dE), x, y);
268                                final float cW = function.apply(Math.abs(dW), x, y);
269                                final float cNE = function.apply(Math.abs(dNE), x, y);
270                                final float cSE = function.apply(Math.abs(dSE), x, y);
271                                final float cSW = function.apply(Math.abs(dSW), x, y);
272                                final float cNW = function.apply(Math.abs(dNW), x, y);
273
274                                image.pixels[y][x] += lambda
275                                                * (cN * dN + cS * dS + cE * dE + cW * dW + wt * (cNE * dNE + cSE * dSE + cSW * dSW + cNW * dNW));
276                        }
277                }
278        }
279
280        // public static void main(String[] args) throws Exception {
281        // final MBFImage img = ImageUtilities.readMBF(new
282        // File("/Users/jon/veg.jpg"));
283        // final AnisotropicDiffusion ad = new AnisotropicDiffusion(20, 0.1f, new
284        // HighConstrastEdgeFunction(30f / 255));
285        //
286        // DisplayUtilities.display(img);
287        // DisplayUtilities.display(img.process(ad));
288        // }
289}