~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/geometry/geomImg_canny.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  * 
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as 
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc., 
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  * 
 25  * Program :    TINA
 26  * File    :  
 27  * Date    :  $Date: 2012/10/02 15:50 $
 28  * Version :  $Revision: 1.9 $
 29  * CVS Id  :  
 30  *
 31  * Author  : Legacy TINA
 32  *
 33  * Notes : HR/NAT modified
 34  *
 35  */
 36  /**
 37  * @file
 38  * @brief TINA implementation of the Canny edge detector.
 39  *
 40  * These functions implement the Canny edge detector.  The algorithm requires two thresholds:
 41  * the edge strength must be greater than the upper threshold in order to initiate detection,
 42  * but then the edge is followed until its strength falls below the lower threshold.  This 
 43  * hysteresis takes account of the fact that, if a neighbouring pixel has been confirmed as 
 44  * part of an edge, this offers additional evidence that the current pixel is part of an edge.
 45  * Two versions of the algorithm are offered here: the function canny is a standard implementation,
 46  * and the function iso_canny prioritises edges whose grey-level is consistent with a user provided 
 47  * grey level, allowing the user to pick the edge type detected (e.g. edge between grey matter and 
 48  * white matter in brain MRI).
 49  */
 50 
 51 #include "geomImg_canny.h"
 52 
 53 #if HAVE_CONFIG_H
 54   #include <config.h>
 55 #endif
 56 
 57 #include <math.h>
 58 #include <float.h>
 59 #include <tina/sys/sysDef.h>
 60 #include <tina/sys/sysPro.h>
 61 #include <tina/math/mathDef.h>
 62 #include <tina/math/mathPro.h>
 63 #include <tina/image/imgDef.h>
 64 #include <tina/image/imgPro.h>
 65 #include <tina/geometry/geom_EdgeDef.h>
 66 #include <tina/geometry/geomImg_nonmax.h>
 67 #include <tina/geometry/geomImg_thresh.h>
 68 #include <tina/geometry/geom_EdgePro.h>
 69 
 70 
 71 /**
 72  * @brief TINA implementation of the standard Canny edge detector.
 73  * @param im Image in which you wish to detect edges.
 74  * @param sigma Standard deviation of a Gaussian smoothing kernel applied to the image (0 = no amoothing).
 75  * @param precision The precision at which the Gaussian smoothing kernel is truncated.
 76  * @param lowthres Lower Canny threshold (in grey levels squared).
 77  * @param upthres Upper Canny threshold (in grey levels squared).
 78  * @param lengththres Minimum allowed length of edges.
 79  * @return edrect Image of the detected edges.
 80  *
 81  * Implementation of the standard Canny edge detector.  The image provided is smoothed using a Gaussian
 82  * kernel defined by the standard deviation "sigma" (0 = no smoothing) and truncated at the probability 
 83  * "precision".  An image of the sum-squared gradients in the x and y directions is then calculated, and 
 84  * edges stronger than the lower threshold are detected using non-maximal suppression. Finally, edges that
 85  * do not feature a voxel stronger than the upper threshold, or are shorter than the length threshold, are
 86  * removed. The edge image is then returned to the calling function.
 87  */
 88 Imrect *canny(Imrect * im, double sigma, double precision, double lowthres)
 89 {
 90     Imrect *gim;
 91     Imrect *edrect;
 92     Imrect *gradx;
 93     Imrect *grady;
 94     Imrect *gradsq;
 95     unsigned int label;         /* for blocked allocation */
 96 
 97     if (im == NULL)
 98         return (NULL);
 99 
100     label = ralloc_end_blocked();
101 
102     if (sigma == 0.0)
103         gim = im_copy(im);
104     else
105         gim = imf_gauss(im, sigma, precision);
106     gradx = imf_grad_h(gim);
107     grady = imf_grad_v(gim);
108     gradsq = imf_sumsq(gradx, grady);
109     im_free(gim);
110 
111     if (label)                  /* allocation was blocked */
112         (void) ralloc_start_blocked(label);     /* re-start blocking */
113 
114     edrect = nonmaxsup(gradx, grady, gradsq, lowthres);
115     im_free(gradx);
116     im_free(grady);
117     im_free(gradsq);
118     return (edrect);
119 }
120 
121 Imrect *canny_link(Imrect *edrect, double upthres, int lengththres, Bool simple)
122 {
123 
124     if (simple == true) er_find_simple_edge_strings(edrect);
125     else er_find_edge_strings(edrect);
126     er_rm_edges(edrect, EDGE_GET_CONN_MASK, EDGE_NOLINK);
127     er_edge_strings_thres(edrect, lengththres, upthres);
128     er_set_row_index(edrect);
129     return (edrect);
130 }
131 
132 
133 static Imrect *get_ssqim(Imrect *im, double midpoint)
134 {
135   Imrect *im1 = NULL, *im2 = NULL, *im3 = NULL, *im4 = NULL, *im5 = NULL;
136   double sigma_x = 0.0, sigma_y = 0.0, noise = 0.0;
137   Imregion *roi = NULL;
138   double factor = 0.0;
139   
140   im1 = im_add(midpoint, im);
141   im2 = im_sqr(im1);
142   im3 = im_sqrt(im2);
143 
144   roi = im3->region;
145   sigma_x = imf_diffx_noise(im3, NULL);
146   sigma_y = imf_diffy_noise(im3, NULL);
147     
148   noise = (sigma_x + sigma_y)/2;
149   factor = -1/(noise*10);
150   
151   im4 = im_times(factor, im3);
152   im5 = im_exp(im4);
153   imf_times_inplace(30.0, im5);
154   
155   im_free(im1);
156   im_free(im2);
157   im_free(im3);
158   im_free(im4);
159 
160   return(im5);
161 
162 }
163 
164 Imrect *iso_canny(Imrect * im, double sigma, double precision, double lowthres, double midpoint)
165 {
166     Imrect *gim;
167     Imrect *edrect;
168     Imrect *gradx;
169     Imrect *grady;
170     Imrect *ssqim;
171     unsigned int label;         /* for blocked allocation */
172 
173 
174     if (im == NULL)
175         return (NULL);
176 
177     label = ralloc_end_blocked();
178 
179 
180 /* without gaussian smoothing to limit the gradient we can loose high gradient iso contours NAT 12/10/08 
181     tsim = im_tsmooth(im);
182 */
183 
184     if (sigma == 0.0)
185         gim = im_tsmooth(im);
186     else
187         gim = imf_gauss(im, sigma, precision);
188     gradx = imf_grad_h(gim);
189     grady = imf_grad_v(gim);
190     ssqim = get_ssqim(gim, midpoint);
191 
192     if (label)                  /* allocation was blocked */
193         (void) ralloc_start_blocked(label);     /* re-start blocking */
194 
195     edrect = nonmaxsup(gradx, grady, ssqim, lowthres);
196     im_free(gim);
197     im_free(gradx);
198     im_free(grady);
199     im_free(ssqim);
200     return (edrect);
201 }
202 
203 void gradient_im(Imrect *im, Imrect **grad_im_x, Imrect **grad_im_y)
204 {
205    float             *row1, *row2, *row3, *row4, *row5, *swap; 
206    int               lx,ly,ux,uy;
207    int               i,j;
208    double            pixvala, pixvalb, pixvalc, mag;
209    float             theta;
210    
211    if(im==NULL) return;
212 
213    lx=im->region->lx;
214    ly=im->region->ly;
215    ux=im->region->ux;
216    uy=im->region->uy;
217    
218    row1 = fvector_alloc(lx, ux);
219    row2 = fvector_alloc(lx, ux);
220    row3 = fvector_alloc(lx, ux);
221    row4 = fvector_alloc(lx, ux);
222    row5 = fvector_alloc(lx, ux);
223 
224    *grad_im_x=im_alloc(im->height,im->width,im->region,float_v);
225    *grad_im_y=im_alloc(im->height,im->width,im->region,float_v);
226    
227     im_get_rowf(row1, im, ly, lx, ux);
228     im_get_rowf(row2, im, ly+1, lx, ux);
229     im_get_rowf(row3, im, ly+2, lx, ux);
230     im_get_rowf(row4, im, ly+3, lx, ux);
231     for(i=ly+2;i<uy-2;++i)
232     {
233         im_get_rowf(row5, im, i+2, lx, ux);
234         for(j=lx+2;j<ux-2;++j)
235         {
236             pixvala = -0.25*(2.0*row3[j]-row3[j-2]-row3[j+2]);
237             pixvalc = -0.25*(2.0*row3[j]-row1[j]-row5[j]);
238 /* cross derivative is sqrt(2) bigger than expected.. factor of 2 -> sqrt(2) here .. NAT */
239             pixvalb = 0.3535*(row2[j-1]-row2[j+1]-row4[j-1]+row4[j+1]);
240 
241             theta = atan2(pixvalb,pixvala-pixvalc);
242             if (theta < 0.0) theta += 2.0*PI;
243             theta /= 2.0;
244 /* etimate signal to noise on orientation and scale gradient to be consistent with step edge
245    estimates for later use NAT */
246             mag = 3.0*sqrt(pixvalb*pixvalb+(pixvala-pixvalc)*(pixvala-pixvalc));
247             IM_PIX_SET(*grad_im_x,i,j,mag*sin(theta));
248             IM_PIX_SET(*grad_im_y,i,j,-mag*cos(theta));
249         }
250         swap = row1;
251         row1 = row2;
252         row2 = row3;
253         row3 = row4;
254         row4 = row5;
255         row5 = swap;
256     }
257     fvector_free(row1, lx);
258     fvector_free(row2, lx);
259     fvector_free(row3, lx);
260     fvector_free(row4, lx);
261     fvector_free(row5, lx);
262 }
263 
264 Imrect *var_canny(Imrect * im, double sigma, double precision, double lowthres)
265 {
266     Imrect *im1, *im2, *im3, *im4, *im5, *im6, *im7;
267     Imrect *edrect;
268     Imrect *gradx, *grady;
269     unsigned int label;         /* for blocked allocation */
270 
271     if (im == NULL)
272         return (NULL);
273 
274     label = ralloc_end_blocked();
275 
276     im1 = imf_sqr(im);
277     im2 = imf_gauss(im1, sigma, precision);
278     im3 = imf_gauss(im, sigma, precision);
279     im4 = imf_sqr(im3);
280     im5 = im_diff(im2,im4);
281     im6 = imf_thresh(0.0, im5);
282     im_free(im1);
283     im_free(im2);
284     im_free(im4);
285     im_free(im5);
286     im7 = imf_sqrt(im6);
287     gradient_im(im7, &gradx, &grady);
288     im_free(im3);
289 
290     if (label)                  /* allocation was blocked */
291         (void) ralloc_start_blocked(label);     /* re-start blocking */
292 
293     edrect = nonmaxsup(gradx, grady, im6, lowthres);
294 
295     im_free(im7);
296     im_free(gradx);
297     im_free(grady);
298     im_free(im6);
299 
300     return (edrect);
301 }
302 
303 Imrect *rank_canny(Imrect * im, double sigma, double precision, double lowthres)
304 {
305     Imrect *im1,*im2;
306     Imrect *edrect;
307     Imrect *gradx;
308     Imrect *grady;
309     unsigned int label;         /* for blocked allocation */
310 
311     if (im == NULL)
312         return (NULL);
313 
314     label = ralloc_end_blocked();
315 
316     im1 = imf_gauss(im, sigma, precision);
317     im2 = im_rank(im1, 5, lowthres*10);
318     
319     gradient_im(im2, &gradx, &grady);
320 
321     if (label)                  /* allocation was blocked */
322         (void) ralloc_start_blocked(label);     /* re-start blocking */
323 
324     edrect = nonmaxsup(gradx, grady, im2, 2.0*lowthres);
325     im_free(gradx);
326     im_free(grady);
327     im_free(im1);
328     im_free(im2);
329     return (edrect);
330 }
331 
332 Imrect *im_canny(Imrect * im, double sigma, double precision, double lowthres)
333 {
334     Imrect *im1;
335     Imrect *edrect;
336     Imrect *gradx;
337     Imrect *grady;
338     unsigned int label;         /* for blocked allocation */
339 
340     if (im == NULL)
341         return (NULL);
342 
343     label = ralloc_end_blocked();
344 
345     im1 = imf_gauss(im, sigma, precision);
346 
347     gradient_im(im1, &gradx, &grady);
348 
349     if (label)                  /* allocation was blocked */
350         (void) ralloc_start_blocked(label);     /* re-start blocking */
351 
352     edrect = nonmaxsup(gradx, grady, im1, 2.0*lowthres);
353     im_free(gradx);
354     im_free(grady);
355     im_free(im1);
356     return (edrect);
357 }
358 
359 /**************************************************************************************************/
360 
361 /* HR: taken from skeleton.c for the application of pairwise tool on the flywing data */
362 
363 Imrect *ridge_canny(Imrect *im, double sigma, double precision, double lowthres) /* HR: changed */
364 {
365     Imrect *edrect, *gradx, *grady, *sig_im;
366     unsigned int label;
367 
368     if (im==NULL)
369         return(NULL);
370 
371     label = ralloc_end_blocked();
372 
373     if (label)
374         (void) ralloc_start_blocked(label);
375 
376     sig_im = imf_gauss_diff(im, sigma, precision, 2.0); /* HR: used this in place of significance_im */
377 
378     gradient_im(sig_im, &gradx, &grady);
379   
380     edrect = nonmaxsup(gradx, grady, sig_im, lowthres);
381 
382     im_free(gradx);
383     im_free(grady);
384     im_free(sig_im);
385 
386     return(edrect);
387 }
388 
389 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.