canny.c Source File

00001 #include "sadie.h"
00002 #include "proto.h"
00003 
00004 static PIXEL firstdiff_y[3][3] = { { 0.0,  0.0,  0.0 }, 
00005                                    { 0.0,  1.0,  1.0 }, 
00006                                    { 0.0, -1.0, -1.0 } };
00007 
00008 static PIXEL firstdiff_x[3][3] = { { 0.0,  0.0,  0.0 }, 
00009                                    { 0.0, -1.0,  1.0 }, 
00010                                    { 0.0, -1.0,  1.0 } };
00011 
00012 int calcsector(PIXEL);
00013 
00014 
00015 /*----------------------------------------------------------------------------*/
00016 /*-General Information--------------------------------------------------------*/
00017 /*                                                                            */
00018 /* This function implements the Canny Edge Detector.                          */
00019 /*                                                                            */
00020 /*----------------------------------------------------------------------------*/
00021 /*-Interface Information------------------------------------------------------*/
00022 void CANNY (
00023 IMAGE  *in,      /*  I   Pointer to the input image.                          */
00024 double sigma,    /*  I   Std. deviation of gaussian smoothing filter.         */
00025 int    size,     /*  I   Size of gaussian smoothing filter.                   */
00026 IMAGE  **out     /*  O   Address of a pointer to the output image             */
00027                  /*      containing an edge map.                              */
00028 /*----------------------------------------------------------------------------*/
00029 ) { register short i, j, k, n;
00030     char msg[SLEN];
00031     IMAGE *gaussian, *smoothed;
00032     IMAGE *FoG_mag, *FoG_dir;
00033     PIXEL *psf=0, factor;
00034     PIXEL theta;
00035     int sectorval;
00036         
00037 
00038     /* if (TIMES) TIMING(T_CANNY); */
00039     if (NAMES) {
00040         MESSAGE('I',"");
00041         MESSAGE('I',"CANNY");
00042         MESSAGE('I',"");
00043         sprintf(msg," Input image:                       %s",in->text);
00044         MESSAGE('I',msg);
00045         sprintf(msg," Gaussian smoothing filter size:    %d",size);
00046         MESSAGE('I',msg);
00047         sprintf(msg,"             Standard Deviation:    %f",sigma);
00048         MESSAGE('I',msg);
00049         MESSAGE('I'," ...............");
00050     }
00051 
00052     /* check input */
00053     if (!CHECKIMG(in)) {
00054         MESSAGE('E'," Can't identify image.");
00055         goto the_end;
00056     } else if (in->nbnd > 1) {
00057         MESSAGE('E'," Image must be single-band.");
00058         goto the_end;
00059     }
00060 
00061 
00062     /* create image of appropriate size */
00063     if (!CHECKIMG(*out)) GETMEM(in->nbnd,in->nlin,in->npix,out);
00064     if (!*out) goto the_end;
00065 
00066 
00067     /* create gaussian smoothing filter */
00068     CREATEGAUSS(size, sigma, &gaussian);
00069     if (CHECKIMG(gaussian)) {
00070         sprintf(gaussian->text, "%dx%d Gaussian -- sigma=%f", size,size,sigma);
00071     } else {
00072         MESSAGE('E'," Can't identify gaussian filter.");
00073         goto the_end;
00074     }
00075     
00076     /* perform gaussian smoothing on input image */
00077     if (!(psf = (PIXEL *) malloc(size*size*sizeof(PIXEL)))) {
00078         MESSAGE('E'," Memory allocation failed.");
00079         goto the_end;
00080     }
00081     for (n=0, j=0, factor=0.0; j<size; j++) {
00082         for (k=0; k<size; k++) {
00083             psf[n] = gaussian->data[0][j][k];
00084             factor += psf[n++];
00085         }
00086     }
00087     for (n=0, j=0; j<size; j++) {
00088         for (k=0; k<size; k++) {
00089             psf[n++] /= factor;
00090         }
00091     }
00092     SCONVL (in, psf, size, size, &smoothed);
00093     if (CHECKIMG(smoothed)) {
00094         sprintf(smoothed->text, "Smoothed image");
00095     } else {
00096         MESSAGE('E'," Can't identify smoothed image.");
00097         goto the_end;
00098     }
00099     
00100     /* create the First Order Gaussian filtered image */
00101     GRADIENT(smoothed,(PIXEL *)firstdiff_y,(PIXEL *)firstdiff_x,3,&FoG_mag,&FoG_dir);
00102     
00103 
00104     /* perform non-maxima supression */
00105     for (j=0; j<in->nlin; j++) {
00106         for (k=0; k<in->npix; k++) {
00107             theta = FoG_dir->data[0][j][k];
00108             if ((sectorval = calcsector(theta)) < 0) {
00109                 MESSAGE('E'," Invalid sector number!.");
00110                 goto the_end;
00111             }
00112 
00113             if (sectorval == 0) {
00114             
00115                 if (((k>0)&&(FoG_mag->data[0][j][k-1] > FoG_mag->data[0][j][k]))
00116                     || ((k<in->npix-1)&&(FoG_mag->data[0][j][k+1] > FoG_mag->data[0][j][k]))) {
00117                         (*out)->data[0][j][k] = 0;
00118                 } else {
00119                     (*out)->data[0][j][k] = FoG_mag->data[0][j][k];
00120                 }
00121             
00122             } else if (sectorval == 1) {
00123             
00124                 if (((j>0 && k<in->npix-1)&&(FoG_mag->data[0][j-1][k+1] > FoG_mag->data[0][j][k]))
00125                     || ((j<in->nlin-1 && k>0)&&(FoG_mag->data[0][j+1][k-1] > FoG_mag->data[0][j][k]))) {
00126                         (*out)->data[0][j][k] = 0;
00127                 } else {
00128                     (*out)->data[0][j][k] = FoG_mag->data[0][j][k];
00129                 }
00130             
00131             } else if (sectorval == 2) {
00132             
00133                 if (((j>0)&&(FoG_mag->data[0][j-1][k] > FoG_mag->data[0][j][k]))
00134                     || ((j<in->nlin-1)&&(FoG_mag->data[0][j+1][k] > FoG_mag->data[0][j][k]))) {
00135                         (*out)->data[0][j][k] = 0;
00136                 } else {
00137                     (*out)->data[0][j][k] = FoG_mag->data[0][j][k];
00138                 }
00139             
00140             } else { /*if (sectorval == 3)*/
00141             
00142                 if (((j>0 && k>0)&&(FoG_mag->data[0][j-1][k-1] > FoG_mag->data[0][j][k]))
00143                     || ((j<in->nlin-1 && k<in->npix-1)&&(FoG_mag->data[0][j+1][k+1] > FoG_mag->data[0][j][k]))) {
00144                         (*out)->data[0][j][k] = 0;
00145                 } else {
00146                     (*out)->data[0][j][k] = FoG_mag->data[0][j][k];
00147                 }
00148             
00149             }
00150         }
00151     }
00152         
00153     the_end:
00154     if (psf) free(psf);
00155     if (CHECKIMG(gaussian))  RELIMG(&gaussian);
00156     if (CHECKIMG(smoothed))  RELIMG(&smoothed);
00157     if (CHECKIMG(FoG_mag))  RELIMG(&FoG_mag);
00158     if (CHECKIMG(FoG_dir))  RELIMG(&FoG_dir);
00159     /* if (TIMES) TIMING(T_EXIT); */
00160 }
00161 
00162 
00163 
00164 /*----------------------------------------------------------------------------*/
00165 /*-General Information--------------------------------------------------------*/
00166 /*                                                                            */
00167 /*   This function calculates the sector of an angle as shown:                */
00168 /*                                                                            */
00169 /*          337.5 -  22.5 degrees   --      Sector 0                          */
00170 /*          22.5  -  67.5 degrees   --      Sector 1                          */
00171 /*          67.5  - 112.5 degrees   --      Sector 2                          */
00172 /*          112.5 - 157.5 degrees   --      Sector 3                          */
00173 /*          157.5 - 202.5 degrees   --      Sector 0                          */
00174 /*          202.5 - 247.5 degrees   --      Sector 1                          */
00175 /*          247.5 - 292.5 degrees   --      Sector 2                          */
00176 /*          292.5 - 337.5 degrees   --      Sector 3                          */
00177 /*                                                                            */
00178 /*                                                                            */
00179 /*----------------------------------------------------------------------------*/
00180 /*-Interface Information------------------------------------------------------*/
00181 int calcsector (
00182 PIXEL  theta      /*  I   Input angle in radians                              */
00183 /*----------------------------------------------------------------------------*/
00184 ) {    
00185     int sectornum;
00186 
00187     /* Verify that this is a gradient direction */
00188     if (theta < -2.0*PI  ||  2.0*PI < theta) {
00189         MESSAGE('E'," calcsector was given an invalid direction angle!");
00190         MESSAGE('E'," Direction must be between -360 and +360 degrees.");
00191         return -1;
00192     }
00193 
00194     /* Convert to positive angle */
00195     if (theta < 0.0) {
00196         theta = (2.0*PI) + theta;
00197     }
00198 
00199     if (theta < (PI / 8.0)  &&  0 <= theta) {
00200         sectornum = 0;
00201     } else if (theta <= (2*PI)  &&  ((15.0*PI)/8.0) <= theta) {
00202         sectornum = 0;
00203     } else if (theta < ((3.0*PI)/8.0)  &&  (PI/8.0) <= theta) {
00204         sectornum = 1;
00205     } else if (theta < ((5.0*PI)/8.0)  &&  ((3.0*PI)/8.0) <= theta) {
00206         sectornum = 2;
00207     } else if (theta < ((7.0*PI)/8.0)  &&  ((5.0*PI)/8.0) <= theta) {
00208         sectornum = 3;
00209     } else if (theta < ((9.0*PI)/8.0)  &&  ((7.0*PI)/8.0) <= theta) {
00210         sectornum = 0;
00211     } else if (theta < ((11.0*PI)/8.0)  &&  ((9.0*PI)/8.0) <= theta) {
00212         sectornum = 1;
00213     } else if (theta < ((13.0*PI)/8.0)  &&  ((11.0*PI)/8.0) <= theta) {
00214         sectornum = 2;
00215     } else {  /* if (theta < ((15.0*PI)/8.0)  &&  ((13.0*PI)/8.0) <= theta) */
00216         sectornum = 3;
00217     }
00218 
00219     return sectornum;
00220 }
00221