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.analysis.algorithm;
031
032import java.io.File;
033import java.io.IOException;
034
035import org.openimaj.image.DisplayUtilities;
036import org.openimaj.image.FImage;
037import org.openimaj.image.ImageUtilities;
038import org.openimaj.image.analyser.ImageAnalyser;
039
040/**
041 * See http://people.cs.uchicago.edu/~pff/papers/dt.pdf
042 * 
043 * An efficient euclidean distance transform applicable to all greyscale images. The distance of each pixel to the closest 
044 * valid pixel is given. In this case a pixel is considered valid when it is less than Float.MAX_VALUE.
045 * 
046 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
047 */
048public class EuclideanDistanceTransform implements ImageAnalyser<FImage> {
049        FImage distances;
050        int [][] indices;
051
052        /* (non-Javadoc)
053         * @see org.openimaj.image.analyser.ImageAnalyser#analyseImage(org.openimaj.image.Image)
054         */
055        @Override
056        public void analyseImage(FImage image) {
057                if (distances == null || distances.height != image.height || distances.width != distances.height) {
058                        distances = new FImage(image.width, image.height);
059                        indices = new int[image.height][image.width];
060                }
061
062                squaredEuclideanDistance(image, distances, indices);
063        }
064
065        /**
066         * @return the distance transformed image (unnormalised)
067         */
068        public FImage getDistances() {
069                return distances;
070        }
071
072        /**
073         * @return the indecies of the closest pixel to any given pixel
074         */
075        public int[][] getIndices() {
076                return indices;
077        }
078
079        protected static void DT1D(float [] f, float [] d, int [] v, int [] l, float [] z) {
080                int k = 0;
081
082                v[0] = 0;
083                z[0] = -Float.MAX_VALUE;
084                z[1] = Float.MAX_VALUE;
085
086                for (int q = 1; q < f.length; q++)
087                {
088                        float s  = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * q - 2 * v[k]);
089
090                        while (s <= z[k]) {
091                                k--;
092                                s  = ((f[q] + q * q) - (f[v[k]] + v[k] * v[k])) / (2 * q - 2 * v[k]);
093                        }
094                        k++;
095                        v[k] = q;
096                        z[k] = s;
097                        z[k + 1] = Float.MAX_VALUE;
098                }
099
100                k = 0;
101                for (int q = 0; q < f.length; q++)
102                {
103                        while (z[k + 1] < q)
104                                k++;
105
106                        d[q] = (q - v[k]) * (q - v[k]) + f[v[k]];
107                        l[q] = v[k];
108                }
109        }
110
111        /**
112         * Calculate the squared euclidean distance transform of a binary image with
113         * foreground pixels set to 1 and background set to 0.
114         * @param image the image to be transformed.
115         * @param distances the distance of each pixel to the closest 1-pixel
116         * @param indices the index of the closes valid pixel
117         */
118        public static void squaredEuclideanDistanceBinary(FImage image, FImage distances, int[][] indices) {
119                float [] f = new float[Math.max(image.height, image.width)];
120                float [] d = new float[f.length];
121                int [] v = new int[f.length];
122                int [] l = new int[f.length];
123                float [] z = new float[f.length + 1];
124
125                for (int x=0; x<image.width; x++) {
126                        for (int y=0; y<image.height; y++) {
127                                f[y] = image.pixels[y][x] == 0 ? Float.MAX_VALUE : 0;
128                        }
129
130                        DT1D(f, d, v, l, z);
131                        for (int y = 0; y < image.height; y++) {
132                                distances.pixels[y][x] = d[y];
133                                indices[y][x] = (l[y] * image.width) + x; //this is now row-major
134                        }
135                }
136                
137                for (int y = 0; y < image.height; y++) {
138                        DT1D(distances.pixels[y], d, v, l, z);
139
140                        for (int x = 0; x < image.width; x++)
141                                l[x] = indices[y][l[x]];
142
143                        for (int x = 0; x < image.width; x++)
144                        {
145                                distances.pixels[y][x] = d[x];
146                                indices[y][x] = l[x];
147                        }
148                }
149        }
150        
151        /**
152         * The static function which underlies EuclideanDistanceTransform. Provide an image, fill distances and indices with 
153         * the distance image and the closest pixel indices. Typically, for the binary case, valid pixels are set to 0 and invalid
154         * pixels are set to Float.MAX_VALUE or Float.POSITIVE_INFINITY.
155         * 
156         * @param image the image to be transformed. Each pixel is considered valid except those of value Float.MAX_VALUE
157         * @param distances the distance of each pixel to the closest non-Float.MAX_VALUE pixel
158         * @param indices the index of the closes valid pixel
159         */
160        public static void squaredEuclideanDistance(FImage image, FImage distances, int[][] indices) {
161                float [] f = new float[Math.max(image.height, image.width)];
162                float [] d = new float[f.length];
163                int [] v = new int[f.length];
164                int [] l = new int[f.length];
165                float [] z = new float[f.length + 1];
166
167                for (int x=0; x<image.width; x++) {
168                        for (int y=0; y<image.height; y++) {
169                                f[y] = Float.isInfinite(image.pixels[y][x]) ? (image.pixels[y][x] > 0 ? Float.MAX_VALUE : -Float.MAX_VALUE) : image.pixels[y][x];
170                        }
171
172                        DT1D(f, d, v, l, z);
173                        for (int y = 0; y < image.height; y++) {
174                                distances.pixels[y][x] = d[y];
175                                indices[y][x] = (l[y] * image.width) + x; //this is now row-major
176                        }
177                }
178                
179                for (int y = 0; y < image.height; y++) {
180                        DT1D(distances.pixels[y], d, v, l, z);
181
182                        for (int x = 0; x < image.width; x++)
183                                l[x] = indices[y][l[x]];
184
185                        for (int x = 0; x < image.width; x++)
186                        {
187                                distances.pixels[y][x] = d[x];
188                                indices[y][x] = l[x];
189                        }
190                }
191        }
192        
193        /**
194         * Test the distance transform
195         * @param args
196         * @throws IOException
197         */
198        public static void main(String args[]) throws IOException{
199                FImage i = ImageUtilities.readF(new File("/Users/ss/Desktop/tache.jpg"));
200                EuclideanDistanceTransform etrans = new EuclideanDistanceTransform();
201//              i.processInplace(new CannyEdgeDetector());
202                i.inverse();
203                for(int x = 0;x < i.width; x++)
204                        for(int y = 0; y < i.height; y++) 
205                                if(i.pixels[y][x] == 1.0f) 
206                                        i.setPixel(x, y, Float.MAX_VALUE);
207                DisplayUtilities.display(i);
208                i.analyseWith(etrans);
209                i = etrans.getDistances();
210                i.normalise();
211                DisplayUtilities.display(i);
212        }
213}