hough.c Source File

00001 #include "sadie.h"
00002 #include "proto.h"
00003 
00004 
00005 
00006 #define FALSE 0
00007 #define TRUE  1
00008 
00009 /* variables necessary for drawing lines with bresh algorithm */
00010 extern int **x_array_buffer;     /* The data where info is stored      */
00011 extern int point;                /* Counter for points used            */
00012 extern int count_point;          /* Flag decide to count or fill array */
00013 
00014 
00015 
00016 
00017 /*----------------------------------------------------------------------------*/
00018 /*-General Information--------------------------------------------------------*/
00019 /*                                                                            */
00020 /* This function implements the HOUGH transform.                              */
00021 /*                                                                            */
00022 /*----------------------------------------------------------------------------*/
00023 /*-Interface Information------------------------------------------------------*/
00024 void HOUGH (
00025 IMAGE  *in,      /*  I   Pointer to the input image.                          */
00026 int    threshold,/*  I   Threshold to find lines of minimum "length".         */
00027 int    rhobins,  /*  I   Number of rho bins in (theta,rho) space.             */
00028 int    thetabins,/*  I   Number of theta bins in (theta,rho) space.           */
00029 IMAGE  **param,  /*  O   Address of a pointer to the parameter space image    */
00030 IMAGE  **out     /*  O   Address of a pointer to the output image             */
00031 /*----------------------------------------------------------------------------*/
00032 ) { register short i, j, k, n;
00033     char msg[SLEN];
00034     PIXEL rho, theta;
00035     PIXEL rhorange;
00036     int rhobucket, thetabucket;
00037     int startx, starty, endx, endy;
00038         
00039 
00040     /* if (TIMES) TIMING(T_HOUGH); */
00041     if (NAMES) {
00042         MESSAGE('I',"");
00043         MESSAGE('I',"HOUGH");
00044         MESSAGE('I',"");
00045         sprintf(msg," Input image:                       %s",in->text);
00046         MESSAGE('I',msg);
00047         sprintf(msg," Number of theta bins:    %d",thetabins);
00048         MESSAGE('I',msg);
00049         sprintf(msg," Number of rho bins:      %d",rhobins);
00050         MESSAGE('I',msg);
00051         sprintf(msg," Length threshold:        %d",threshold);
00052         MESSAGE('I',msg);
00053         MESSAGE('I'," ...............");
00054     }
00055 
00056     /* check input */
00057     if (!CHECKIMG(in)) {
00058         MESSAGE('E'," Can't identify image.");
00059         goto the_end;
00060     } else if (in->nbnd > 1) {
00061         MESSAGE('E'," Image must be single-band.");
00062         goto the_end;
00063     } else if (thetabins < 1) {
00064         MESSAGE('E'," Number of theta bins must be >= 1.");
00065         goto the_end;
00066     } else if (rhobins < 1) {
00067         MESSAGE('E'," Number of rho bins must be >= 1.");
00068         goto the_end;
00069     }
00070 
00071 
00072     /* create images of appropriate size */
00073     if (!CHECKIMG(*out)) GETMEM(in->nbnd,in->nlin,in->npix,out);
00074     if (!*out) goto the_end;
00075     if (!CHECKIMG(*param)) GETMEM((short)1,(short)rhobins,(short)thetabins,param);
00076     if (!*param) goto the_end;
00077 
00078 
00079     /* Initialize the parameter space */
00080     for (j=0; j<rhobins; j++) {
00081         for (k=0; k<thetabins; k++) {
00082             (*param)->data[0][j][k] = (PIXEL)0.0;
00083         }    
00084     }
00085     
00086     /* Initialize the Hough Image */
00087     for (j=0; j<in->nlin; j++) {
00088         for (k=0; k<in->npix; k++) {
00089             (*out)->data[0][j][k] = (PIXEL)0.0;
00090         }    
00091     }
00092 
00093     /* Determine the range of rho */
00094     rhorange = (PIXEL) sqrt((double)((in->nlin)*(in->nlin) + (in->npix)*(in->npix)));
00095     
00096     /* Fill the parameter space */
00097     RANGE(in);
00098     for (j=0; j<in->nlin;j++) {
00099         for (k=0; k<in->npix; k++) {
00100             if(in->data[0][j][k] > in->gmin) {
00101                 for (thetabucket=0; thetabucket<thetabins; thetabucket++) {
00102 
00103                     /* Find theta */
00104                     if (thetabucket <= (thetabins/2)) {
00105                         theta = (PIXEL) (((PIXEL)thetabucket/(thetabins/2.0) - 1.0) * PI);
00106                     } else {
00107                         theta = (PIXEL)((((PIXEL)thetabucket - (PIXEL)(thetabins%2) + 1.0)/(thetabins/2.0) - 1.0) * PI);
00108                     }
00109 
00110                     /* Find rho */
00111                     rho = (PIXEL) (k*cos((double)theta) + j*sin((double)theta));
00112                     
00113                     if (rho>0.0) {
00114                         rhobucket = (rhobins/2) + (int)((rho/rhorange)*(PIXEL)(rhobins/2)) + (rhobins%2) - 1;
00115                     } else {
00116                         rhobucket = (rhobins/2) + (int)((rho/rhorange)*(PIXEL)(rhobins/2));
00117                     }
00118                     
00119                     
00120                     /* Increment the appropriate bucket in theta-rho space */
00121                     (*param)->data[0][rhobucket][thetabucket] += 1.0;
00122                 }
00123             }
00124         }    
00125     }
00126 
00127 
00128     for (j=0; j<(*param)->nlin;j++) {
00129         for (k=0; k<(*param)->npix; k++) {
00130             if((*param)->data[0][j][k] >= threshold) {
00131                 
00132                 /* Find the theta and rho for this point */
00133                 if (k <= (thetabins/2)) {
00134                     theta = (PIXEL) (((PIXEL)k/(thetabins/2.0) - 1.0) * PI);
00135                 } else {
00136                     theta = (PIXEL)((((PIXEL)k - (PIXEL)(thetabins%2) + 1.0)/(thetabins/2.0) - 1.0) * PI);
00137                 }
00138                 
00139                 if (j <= (rhobins/2)) {
00140                     rho = (PIXEL)(((PIXEL)j/(rhobins/2.0) - 1.0) * rhorange);
00141                 } else {
00142                     rho = (PIXEL)((((PIXEL)j - (PIXEL)(rhobins%2) + 1.0)/(rhobins/2.0) - 1.0) * rhorange);
00143                 }
00144                 
00145                 
00146                 if (theta == (PIXEL)0.0) {
00147                     startx = endx = (int) rnd((double)rho);
00148                     starty = 0;
00149                     endy = in->nlin-1;
00150                 } else {
00151                     /* Find the start and end points of the line */
00152                     startx = 0;
00153                     starty = (int) rnd((double)((double)rho - (double)startx*cos((double)theta)) / sin((double)theta));
00154 
00155                     /* Case I: Line touches left edge of image */                
00156                     if (starty >= 0 && starty < in->nlin) {
00157                         /* ...and also touches right edge... */
00158                         endx = in->npix - 1;
00159                         endy = (int) rnd((double)((double)rho - (double)endx*cos((double)theta)) / sin((double)theta));
00160 
00161                         if (endy < 0 || endy >= in->nlin) {
00162                             /* ...and also touches top edge... */
00163                             endy = 0;
00164                             endx = (int) rnd((double)((double)rho - (double)endy*sin((double)theta)) / cos((double)theta));
00165 
00166                             if (endx < 0 || endx >= in->npix) {
00167                                 /* ...and also touches bottom edge... */
00168                                 endy = in->nlin - 1;
00169                                 endx = (int) rnd((double)((double)rho - (double)endy*sin((double)theta)) / cos((double)theta));
00170                             }
00171                         }
00172                     } else {
00173                         starty = 0;
00174                         startx = (int) rnd((double)((double)rho - (double)starty*sin((double)theta)) / cos((double)theta));
00175 
00176                         /* Case II: Line touches top edge of image */                
00177                         if (startx >= 0 && startx < in->npix) {
00178                             /* ...and also touches right edge... */
00179                             endx = in->npix - 1;
00180                             endy = (int) rnd((double)((double)rho - (double)endx*cos((double)theta)) / sin((double)theta));
00181 
00182                             if (endy < 0 || endy >= in->nlin) {
00183                                 /* ...and also touches bottom edge... */
00184                                 endy = in->nlin - 1;
00185                                 endx = (int) rnd((double)((double)rho - (double)endy*sin((double)theta)) / cos((double)theta));
00186                             }
00187                         } else {
00188                             /* Case III: Line touches bottom edge of image */                
00189                             starty = in->nlin - 1;
00190                             startx = (int) rnd((double)((double)rho - (double)starty*sin((double)theta)) / cos((double)theta));
00191 
00192                             if (startx >= 0 && startx < in->npix) {
00193                                 /* ...and also touches right edge... */
00194                                 endx = in->npix - 1;
00195                                 endy = (int) rnd((double)((double)rho - (double)endx*cos((double)theta)) / sin((double)theta));
00196 
00197                                 if (endy < 0 || endy >= in->nlin) {
00198                                     MESSAGE('E'," Could not find line start and end points!");
00199                                     sprintf(msg,"     rhobucket = %d, rho = %f", j, rho);
00200                                     MESSAGE('E',msg);
00201                                     sprintf(msg,"     thetabucket = %d, theta = %f", k, theta);
00202                                     MESSAGE('E',msg);
00203                                     goto the_end;
00204                                 }
00205                             } else {
00206                                     MESSAGE('E'," Could not find line start and end points!");
00207                                     sprintf(msg,"     rhobucket = %d, rho = %f", j, rho);
00208                                     MESSAGE('E',msg);
00209                                     sprintf(msg,"     thetabucket = %d, theta = %f", k, theta);
00210                                     MESSAGE('E',msg);
00211                                     goto the_end;
00212                             }
00213                         }
00214 
00215                     }
00216                 }
00217 
00218                 printf("rhobucket = %d, rho = %f\n", j, rho);
00219                 printf("thetabucket = %d, theta = %f\n", k, theta);
00220                 printf("Line = (%d,%d) --> (%d,%d)\n\n", startx, starty, endx, endy);
00221 
00222 
00223                 /* Count number of points in line */
00224                 count_point = TRUE;
00225                 point = 0;
00226                 brshnm(startx,starty,endx, endy);
00227 
00228                 /* Allocate memory */
00229                 x_array_buffer = (int **) malloc(point*sizeof(int  * ));
00230                 if (x_array_buffer == (int **) NULL){
00231                     printf("No memory for x_array_buffer\n");
00232                 } else{
00233 
00234                     for  (i = 0; i < point; i++){
00235                         x_array_buffer[i] = (int  *) calloc(2, sizeof(int));
00236                         if (x_array_buffer[i] == (int) NULL){
00237                             printf ("No memory for element %d\n",x_array_buffer[i]);
00238                         }
00239                     }
00240 
00241                     /* Fill array */ 
00242                     count_point = FALSE;
00243                     point = 0;
00244                     brshnm(startx,starty,endx, endy);
00245 
00246                         for (i=0; i < point; i++) {
00247                              (*out)->data[0][x_array_buffer[i][1]][x_array_buffer[i][0]] = (PIXEL)1.0;
00248                         }
00249 
00250                     if(x_array_buffer) {
00251                         /* Free Bresh Array */
00252                         for (i=0; i < point; i++){
00253                             if(x_array_buffer[i]) {
00254                                 free(x_array_buffer[i]);
00255                                 x_array_buffer[i] = NULL;
00256                             }
00257                         }
00258 
00259                         free(x_array_buffer);
00260                         x_array_buffer = NULL;
00261                     }
00262 
00263                 }
00264             }
00265         }
00266     }
00267     
00268         
00269     the_end:
00270     /* if (TIMES) TIMING(T_EXIT); */
00271 
00272     if(x_array_buffer) {
00273         /* Free Bresh Array */
00274         for (i=0; i < point; i++){
00275             if(x_array_buffer[i]) {
00276                 free(x_array_buffer[i]);
00277                 x_array_buffer[i] = NULL;
00278             }
00279         }
00280 
00281         free(x_array_buffer);
00282         x_array_buffer = NULL;
00283     }
00284 }
00285 
00286 
00287 
00288 
00289 
00290