hough.c语言

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; j00081 for (k=0; k00082 (*param)->data[0][j][k] = (PIXEL)0.0;
00083 }
00084 }
00085
00086 /* Initialize the Hough Image */
00087 for (j=0; jnlin; j++) {
00088 for (k=0; knpix; 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; jnlin;j++) {
00099 for (k=0; knpix; k++) {
00100 if(in->data[0][j][k] > in->gmin) {
00101 for (thetabucket=0; thetabucket00102
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 {
001

36 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