Source code

001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 *  image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause 
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 
006 * Copyright (c) 2006-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 
007 * Visit http://imagingbook.com for additional details.
008 *******************************************************************************/
009
010package imagingbook.pub.color.edge;
011
012import ij.plugin.filter.Convolver;
013import ij.process.ByteProcessor;
014import ij.process.ColorProcessor;
015import ij.process.FloatProcessor;
016import ij.process.ImageProcessor;
017
018import java.awt.Point;
019import java.util.LinkedList;
020import java.util.List;
021import java.util.Stack;
022
023/**
024 * This class implements a Canny edge detector for grayscale and RGB images.
025 * @author W. Burger
026 * @version 2016/05/04
027 */
028public class CannyEdgeDetector extends ColorEdgeDetector {
029        
030        public static class Parameters {
031                /** Gaussian sigma (scale) */
032                public double gSigma = 2.0f;
033                /** High threshold (20% of max. edge magnitude) */
034                public double hiThr  = 20.0f;
035                /** Low threshold (5% of max. edge magnitude) */
036                public double loThr = 5.0f;
037                /** Set true to normalize gradient magnitude */
038                public boolean normGradMag = true;      
039                
040                /**
041                 * Checks the parameter set.
042                 * @return true if any invalid condition is found
043                 */
044                public boolean isInValid () { 
045                        return gSigma < 0.1f || loThr > hiThr;
046                }
047        }
048        
049        final Parameters params;
050        final ImageProcessor I;                         // original image
051        final int M, N;                                         // width and height of I
052        
053        FloatProcessor Emag = null;                             // gradient magnitude
054        FloatProcessor Enms = null;                             // non-max suppressed gradient magnitude
055        FloatProcessor Ex = null, Ey = null;    // edge normal vectors
056        ByteProcessor Ebin = null;                              // final (binary) edge image
057        List<List<Point>> traceList = null;             // list of edge traces
058        
059        // Constructor with default parameters:
060        public CannyEdgeDetector(ImageProcessor ip) {
061                this(ip, new Parameters());
062        }
063        
064        // Constructor with parameter object:
065        public CannyEdgeDetector(ImageProcessor I, Parameters params) {
066                if (params.isInValid()) throw new IllegalArgumentException();
067                this.params = params;
068                this.I = I;
069                M = I.getWidth();
070                N = I.getHeight();
071                findEdges();
072        }
073        
074        // do the work ...
075        void findEdges() {
076                if (I instanceof ColorProcessor) 
077                        makeGradientsAndMagnitudeColor();
078                else
079                        makeGradientsAndMagnitudeGray();
080                nonMaxSuppression();
081                //detectAndTraceEdges();
082        }
083        
084        //---------------------------------------------------------------------------
085        
086        void makeGradientsAndMagnitudeGray() {
087//              FloatProcessor If = (I instanceof FloatProcessor) ? 
088//                              (FloatProcessor) I.duplicate() :
089//                              (FloatProcessor) I.convertToFloat();
090                
091                FloatProcessor If = I.convertToFloatProcessor();        // always makes a copy
092                                
093                // apply a separable Gaussian filter to I
094                float[] gaussKernel = makeGaussKernel1d(params.gSigma);
095                Convolver conv = new Convolver();
096                conv.setNormalize(true);
097                conv.convolve(If, gaussKernel, gaussKernel.length, 1);
098                conv.convolve(If, gaussKernel, 1, gaussKernel.length);
099                
100                // calculate the gradients in X- and Y-direction
101                Ex = If;
102                Ey = (FloatProcessor) If.duplicate();
103                float[] gradKernel = {-0.5f, 0, 0.5f};
104                conv.setNormalize(false);
105                conv.convolve(Ex, gradKernel, gradKernel.length, 1);
106                conv.convolve(Ey, gradKernel, 1, gradKernel.length);
107                
108                Emag = new FloatProcessor(M, N);
109                float emax = 0;
110                for (int v = 0; v < N; v++) {
111                        for (int u = 0; u < M; u++) {
112                                float dx = Ex.getf(u,v);
113                                float dy = Ey.getf(u,v);
114                                float mag = (float) Math.hypot(dx, dy); // = (float) Math.sqrt(dx*dx + dy*dy);
115                                if (mag > emax) 
116                                        emax = mag;
117                                Emag.setf(u, v, mag);
118                        }
119                }
120                //IJ.log("Gray emax = " + emax);
121                
122                // normalize gradient magnitude 
123                if (params.normGradMag && emax > 0.001) 
124                        Emag.multiply(100.0/emax);
125        }
126        
127        void makeGradientsAndMagnitudeColor() {
128                FloatProcessor[] Irgb = rgbToFloatChannels((ColorProcessor)I);
129                FloatProcessor[] Ixrgb = new FloatProcessor[3];
130                FloatProcessor[] Iyrgb = new FloatProcessor[3];
131                
132                // apply a separable Gaussian filter to each RGB channel
133                float[] gaussKernel = makeGaussKernel1d(params.gSigma);
134                Convolver conv = new Convolver();
135                conv.setNormalize(true);
136                for (int i=0; i<Irgb.length; i++) {
137                        FloatProcessor I = Irgb[i];
138                        conv.convolve(I, gaussKernel, gaussKernel.length, 1);
139                        conv.convolve(I, gaussKernel, 1, gaussKernel.length);
140                        Ixrgb[i] = I;
141                        Iyrgb[i] = (FloatProcessor) I.duplicate();
142                }
143                
144                // calculate the gradients in X- and Y-direction for each RGB channel
145                float[] gradKernel = {-0.5f, 0, 0.5f};
146                conv.setNormalize(false);
147                for (int i = 0; i < Irgb.length; i++) {
148                        FloatProcessor Ix = Ixrgb[i];
149                        FloatProcessor Iy = Iyrgb[i];
150                        conv.convolve(Ix, gradKernel, gradKernel.length, 1);
151                        conv.convolve(Iy, gradKernel, 1, gradKernel.length);
152                }
153
154                // calculate gradient magnitude
155                Ex = new FloatProcessor(M, N);
156                Ey = new FloatProcessor(M, N);
157                Emag = new FloatProcessor(M, N);
158                float emax = 0;
159                for (int v = 0; v < N; v++) {
160                        for (int u = 0; u < M; u++) {
161                                float rx = Ixrgb[0].getf(u,v), ry = Iyrgb[0].getf(u,v);
162                                float gx = Ixrgb[1].getf(u,v), gy = Iyrgb[1].getf(u,v);
163                                float bx = Ixrgb[2].getf(u,v), by = Iyrgb[2].getf(u,v);
164                                float A = rx*rx + gx*gx + bx*bx;
165                                float B = ry*ry + gy*gy + by*by;
166                                float C = rx*ry + gx*gy + bx*by;
167                                float D = (float) Math.sqrt((A - B)*(A - B) + 4*C*C);
168                                
169                                float mag = (float) Math.sqrt(0.5*(A+B+D));
170                                if (mag > emax) emax = mag;
171                                Emag.setf(u, v, mag);
172                                Ex.setf(u, v, A - B + D);
173                                Ey.setf(u, v, 2*C);
174                        }
175                }
176                //IJ.log("RGB emax = " + emax);
177                // normalize gradient magnitude 
178                if (params.normGradMag && emax > 0.001) 
179                        Emag.multiply(100.0/emax);
180        }
181        
182        
183        //---------------------------------------------------------------------------
184        
185        void nonMaxSuppression() {
186                // perform non-maximum suppression along gradient direction
187                Enms = new FloatProcessor(M, N);
188                for (int v = 1; v < N - 1; v++) {
189                        for (int u = 1; u < M - 1; u++) {
190                                int s_theta = getOrientationSector(Ex.getf(u, v), Ey.getf(u, v));
191                                if (isLocalMaximum(Emag, u, v, s_theta, (float) params.loThr)) {
192                                        Enms.setf(u, v, Emag.getf(u, v)); // keep local maximum only
193                                }
194                        }
195                }
196        }
197        
198        void detectAndTraceEdges() {
199                if (Enms == null) {
200                        nonMaxSuppression();
201                }
202                Ebin = new ByteProcessor(M, N);
203                int color = 255;
204                traceList = new LinkedList<List<Point>>();
205                for (int v = 0; v < N; v++) {
206                        for (int u = 0; u < M; u++) {
207                                if (Enms.getf(u, v) >= params.hiThr && Ebin.get(u, v) == 0) { // unmarked edge point
208                                        List<Point> trace = traceAndThreshold(u, v, (float) params.loThr, color);
209                                        traceList.add(trace);
210                                }
211                        }
212                }
213        }
214        
215        // Determines if the gradient magnitude is a local maximum at position (u,v)
216        // in direction s_theta.
217        boolean isLocalMaximum(FloatProcessor gradMagnitude, int u, int v, int s_theta, float mMin) {
218                float mC = gradMagnitude.getf(u, v);
219                if (mC < mMin) {
220                        return false;
221                }
222                else {
223                        float mL = 0, mR = 0;
224                        switch (s_theta) {
225                        case 0 : 
226                                mL = gradMagnitude.getf(u-1, v);
227                                mR = gradMagnitude.getf(u+1, v);
228                                break;
229                        case 1 : 
230                                mL = gradMagnitude.getf(u-1, v-1);
231                                mR = gradMagnitude.getPixelValue(u+1, v+1);
232                                break;
233                        case 2 : 
234                                mL = gradMagnitude.getf(u, v-1);
235                                mR = gradMagnitude.getf(u, v+1);
236                                break;
237                        case 3 : 
238                                mL = gradMagnitude.getf(u-1, v+1);
239                                mR = gradMagnitude.getf(u+1, v-1);
240                                break;
241                        }
242                        return (mL <= mC && mC >= mR);
243                }
244        }
245        
246        ///Recursively collect and mark all pixels of an edge that are 8-connected to 
247        // (u0,v0) and have a gradient magnitude above loThr.
248        /**
249         * Recursively collect and mark all pixels of an edge that are 8-connected to 
250         * (u0,v0) and have a gradient magnitude above loThr.
251         * @param u0 start coordinate (x)
252         * @param v0 start coordinate (y)
253         * @param loThr low threshold (min. edge magnitude to continue tracing)
254         * @param markColor color used for marking pixels on edge trace
255         * @return a list of Point objects.
256         */
257        List<Point> traceAndThreshold(int u0, int v0, float loThr, int markColor) {
258                //IJ.log("  working point " + u + " " + v);
259                Stack<Point> pointStack = new Stack<Point>();
260                List<Point> pointList = new LinkedList<Point>();
261                pointStack.push(new Point(u0, v0));
262                while (!pointStack.isEmpty()) {
263                        Point p = pointStack.pop();
264                        int up = p.x;
265                        int vp = p.y;
266                        Ebin.putPixel(up, vp, markColor);       // mark this edge point
267                        pointList.add(p);
268                                
269                        int uL = Math.max(up - 1, 0);           // (up > 0) ? up-1 : 0;
270                        int uR = Math.min(up + 1, M - 1);       // (up < M-1) ? up+1 : M-1;
271                        int vT = Math.max(vp - 1, 0);           // (vp > 0) ? vp - 1 : 0;
272                        int vB = Math.min(vp + 1, N - 1);       // (vp < N-1) ? vp+1 : N-1;
273                        
274                        for (int u = uL; u <= uR; u++) {
275                                for (int v = vT; v <= vB; v++) {
276                                        if (Ebin.get(u, v) == 0 && Enms.getf(u, v) >= loThr) {
277                                                pointStack.push(new Point(u,v));
278                                        }
279                                }
280                        }
281                }
282                return pointList;
283        }
284        
285        private final float cosPi8 = (float) Math.cos(Math.PI/8);
286        private final float sinPi8 = (float) Math.sin(Math.PI/8);
287        
288        // returns the quantized orientation sector for vector (dx, dy)
289        int getOrientationSector(float dx, float dy) {
290                // rotate the gradient vector by PI/8
291                float dxR = cosPi8 * dx - sinPi8 * dy;
292                float dyR = sinPi8 * dx + cosPi8 * dy;  
293                // mirror vector (dxR,dyR) to [0,PI]
294                if (dyR < 0) {
295                        dxR = -dxR;
296                        dyR = -dyR;
297                }
298                // determine the octant for (dx, dy)
299                int s_theta;
300                if (dxR >= 0) { // dxR >= 0, dyR >= 0
301                        if (dxR >= dyR)
302                                s_theta = 0; // return 0;
303                        else
304                                s_theta = 1; // return 1;
305                }
306                else {  // dxR < 0, dyR >= 0
307                        if (-dxR < dyR)
308                                s_theta = 2; // return 2;
309                        else
310                                s_theta = 3; //return 3;
311                }
312                return s_theta;
313        }
314        
315        public FloatProcessor getEdgeMagnitude() {
316                return Emag;
317        }
318
319        public FloatProcessor getEdgeOrientation() {
320                FloatProcessor E_theta = new FloatProcessor(M, N);
321                for (int u = 0; u < M; u++) {
322                        for (int v = 0; v < N; v++) {
323                                double ex = Ex.getf(u, v);
324                                double ey = Ey.getf(u, v);
325                                float theta = (float) Math.atan2(ey, ex);
326                                E_theta.setf(u, v, theta);
327                        }
328                }
329                return E_theta;
330        }
331        
332        public ByteProcessor getEdgeBinary() {
333                if (Ebin == null) {
334                        detectAndTraceEdges();
335                }
336                return Ebin;
337        }
338        
339        public List<List<Point>> getEdgeTraces() {
340                if (traceList == null) {
341                        detectAndTraceEdges();
342                }
343                return traceList;
344        }
345        
346        //---------------------------------------------------------------------------
347
348        float[] makeGaussKernel1d(double sigma) {
349                // make 1D Gauss filter kernel large enough
350                int rad = (int) (3.5 * sigma);
351                int size = rad + rad + 1;
352                float[] kernel = new float[size]; // center cell = kernel[rad]
353                double sigma2 = sigma * sigma;
354                double scale = 1 / Math.sqrt(2 * Math.PI * sigma2);
355                for (int i = 0; i < size; i++) {
356                        double x = rad - i;
357                        kernel[i] = (float) (
358                        scale * Math.exp(-0.5 * (x * x) / sigma2)
359                        );
360                }
361                return kernel;
362        }
363        
364        // extract RGB channels of 'cp' as 3 float processors
365        FloatProcessor[] rgbToFloatChannels(ColorProcessor cp) {
366                int w = cp.getWidth();
367                int h = cp.getHeight();
368                FloatProcessor rp = new FloatProcessor(w, h);
369                FloatProcessor gp = new FloatProcessor(w, h);
370                FloatProcessor bp = new FloatProcessor(w, h);
371                int[] pixels = (int[]) cp.getPixels(); 
372                float[] rpix = (float[]) rp.getPixels();
373                float[] gpix = (float[]) gp.getPixels();
374                float[] bpix = (float[]) bp.getPixels();
375        for (int i=0; i<pixels.length; i++) {
376            int c = pixels[i];
377            rpix[i] = (c&0xff0000)>>16;
378                gpix[i] = (c&0xff00)>>8;
379                        bpix[i] =  c&0xff;
380        }
381        return new FloatProcessor[] {rp, gp, bp};
382        }
383        
384}