register.c Source File

00001 #include "sadie.h"
00002 #include "proto.h"
00003 #include <math.h>
00004 #include <sys/time.h>
00005 
00006 
00007 
00008 /*----------------------------------------------------------------------------*/
00009 /*-General Information--------------------------------------------------------*/
00010 /*                                                                            */
00011 /* This function compares two edge maps.                                      */
00012 /*                                                                            */
00013 /*----------------------------------------------------------------------------*/
00014 /*-Interface Information------------------------------------------------------*/
00015 void REGISTER (
00016 IMAGE  *img1,     /*  I   Pointer to the first input image.                    */
00017 IMAGE  *img2,     /*  I   Pointer to the second input image.                   */
00018 int    est_x_off, /*  I   Estimated x-offset of second image                   */
00019 int    est_y_off, /*  I   Estimated y-offset of second image                   */
00020 IMAGE  **target,  /*  O   Pointer to the target image.                         */
00021 IMAGE  **search,  /*  O   Pointer to the search image.                         */
00022 IMAGE  **corr,    /*  O   Pointer to the correllation image.                   */
00023 IMAGE  **out      /*  O   Pointer to the output image.                         */
00024 /*----------------------------------------------------------------------------*/
00025 ) { register short i, j, k, m, n, max_i, max_j;
00026     char msg[SLEN];
00027     int yoffset, xoffset, x_index, y_index;
00028     int x_offset[4], y_offset[4];
00029     int roi1_x1, roi1_x2, roi1_y1, roi1_y2;
00030     int roi2_x1, roi2_x2, roi2_y1, roi2_y2;
00031     int roi1_center_x, roi1_center_y, roi2_center_x, roi2_center_y;
00032     PIXEL mean_s, mean_t;
00033     PIXEL denom_s, denom_t, num;
00034     PIXEL temp_s, temp_t, value, maxval;
00035     //int search_size = 111;
00036     //int target_size = 29;
00037     int search_size = 69;
00038     int target_size = 49;
00039     IMAGE *inimg[2];
00040     int jbgn[2], kbgn[2];
00041     int nlin, npix;
00042     PIXEL bias[2], gain[2];
00043     int xcount, ycount, counter;
00044     struct  timeval start, findreg, doreg, end;
00045         
00046     /* start timing */
00047     gettimeofday(&start, NULL);
00048 
00049     /* check input */
00050     if (!CHECKIMG(img1)) {
00051         MESSAGE('E'," Can't identify first input image.");
00052         goto the_end;
00053         }
00054     if (!CHECKIMG(img2)) {
00055         MESSAGE('E'," Can't identify second input image.");
00056         goto the_end;
00057         }
00058     
00059     inimg[0] = img1;
00060     inimg[1] = img2;
00061         
00062 
00063     if (NAMES) {
00064         MESSAGE('I',"");
00065         MESSAGE('I',"REGISTER");
00066         MESSAGE('I',"");
00067         sprintf(msg," First input image:                   %s",img1->text);
00068         MESSAGE('I',msg);
00069         sprintf(msg," Second input image:                  %s",img2->text);
00070         MESSAGE('I',msg);
00071         MESSAGE('I',"");
00072         MESSAGE('I',"");
00073         sprintf(msg," Estimated x-offset:  %d",est_x_off);
00074                 MESSAGE('I',msg);
00075         sprintf(msg," Estimated y-offset:  %d",est_y_off);
00076                 MESSAGE('I',msg);
00077     }
00078 
00079     
00080     if (est_x_off > img1->npix) {
00081             MESSAGE('E'," Images cannot be registered to one another because they do not overlap!");
00082             goto the_end;
00083     }
00084     if (est_y_off > img1->nlin) {
00085             MESSAGE('E'," Images cannot be registered to one another because they do not overlap!");
00086             goto the_end;
00087     }
00088 
00089 
00090     
00091     /* create sub-images of appropriate size */
00092     if (!CHECKIMG(*target)) GETMEM(1,target_size,target_size,target);
00093     if (!*target) goto the_end;
00094     if (!CHECKIMG(*search)) GETMEM(1,search_size,search_size,search);
00095     if (!*search) goto the_end;
00096     if (!CHECKIMG(*corr)) GETMEM(1,(search_size - target_size)+1,(search_size - target_size)+1,corr);
00097     if (!*corr) goto the_end;
00098     
00099     
00100 
00101 
00102     /* Determine bounds on estimated overlapping regions */
00103     if (est_x_off >= 0) {
00104             roi1_x1 = min(est_x_off,img1->npix-1);
00105             roi1_x2 = min(img1->npix-1,(img2->npix-1)+roi1_x1);
00106 
00107             roi2_x1 = 0;
00108             roi2_x2 = min(img2->npix-1,(img1->npix-1)-est_x_off);
00109     } else {
00110             roi1_x1 = 0;
00111             roi1_x2 = min(img1->npix-1,(img2->npix-1) + est_x_off);
00112 
00113             roi2_x1 = min(img2->npix-1, 0-est_x_off);
00114             roi2_x2 = min(img2->npix-1, (img1->npix-1) - est_x_off);
00115     }
00116 
00117     if (est_y_off >= 0) {
00118             roi1_y1 = min(est_y_off,img1->nlin-1);
00119             roi1_y2 = min(img1->nlin-1,(img2->nlin-1)+roi1_y1);
00120 
00121             roi2_y1 = 0;
00122             roi2_y2 = min(img2->nlin-1, (img1->nlin-1)-est_y_off);
00123     } else {
00124             roi1_y1 = 0;
00125             roi1_y2 = min(img1->nlin-1,(img2->nlin-1) + est_y_off);
00126 
00127             roi2_y1 = min(img2->nlin-1, 0-est_y_off);
00128             roi2_y2 = min(img2->nlin-1, (img1->nlin-1) - est_y_off);
00129     }   
00130 
00131 
00132 printf("Overlap size: %d x %d \n",(roi1_x2 - roi1_x1),(roi1_y2 - roi1_y1));
00133 
00134 
00135     /* Make sure there is enough overlapping data to work with */
00136     if (((roi1_x2 - roi1_x1) < (2*search_size))||((roi1_y2 - roi1_y1) < (2*search_size))) {
00137             MESSAGE('E'," Images cannot be registered to one another because the overlap area is too small!");
00138             goto the_end;
00139     }
00140     if (((roi1_x2 - roi1_x1) != (roi2_x2 - roi2_x1))||((roi1_y2 - roi1_y1) != (roi2_y2 - roi2_y1))) {
00141             MESSAGE('E'," Overlapping areas are not the same size!");
00142             goto the_end;
00143     }
00144 
00145 
00146 
00147 
00148     /* note time at beginning of correlation */
00149     gettimeofday(&findreg, NULL);
00150 
00151     mean_t = (PIXEL)0.0;
00152     mean_s = (PIXEL)0.0;
00153     for (xcount = counter = 0; xcount <2; xcount++) {
00154         for (ycount = 0; ycount < 2; ycount++,counter++) {
00155         
00156             roi1_center_x = roi1_x1 + ((roi1_x2 - roi1_x1)/4) + (xcount*((roi1_x2 - roi1_x1)/2));
00157             roi1_center_y = roi1_y1 + ((roi1_y2 - roi1_y1)/4) + (ycount*((roi1_y2 - roi1_y1)/2));
00158             roi2_center_x = roi2_x1 + ((roi2_x2 - roi2_x1)/4) + (xcount*((roi2_x2 - roi2_x1)/2));
00159             roi2_center_y = roi2_y1 + ((roi2_y2 - roi2_y1)/4) + (ycount*((roi2_y2 - roi2_y1)/2));
00160             
00161             
00162             /* Compute mean target value */
00163             for (m=0; m<target_size; m++) {
00164                     for (n=0; n<target_size; n++) {
00165                                 mean_t += img2->data[0][(roi2_center_y-((target_size-1)/2))+m][(roi2_center_x-((target_size-1)/2))+n];
00166                             (*target)->data[0][m][n] = img2->data[0][(roi2_center_y-((target_size-1)/2))+m][(roi2_center_x-((target_size-1)/2))+n];
00167                 }
00168             }
00169 
00170             /* Compute mean search value */
00171             for (m=0; m<search_size; m++) {
00172                     for (n=0; n<search_size; n++) {
00173                                 mean_s += img1->data[0][(roi1_center_y-((search_size-1)/2))+m][(roi1_center_x-((search_size-1)/2))+n];
00174                             (*search)->data[0][m][n] = img1->data[0][(roi1_center_y-((search_size-1)/2))+m][(roi1_center_x-((search_size-1)/2))+n];
00175                     }
00176             }
00177 
00178         }
00179     }
00180     mean_t /= (PIXEL)(4 * target_size * target_size);
00181     mean_s /= (PIXEL)(4 * search_size * search_size);
00182     
00183     denom_t = (PIXEL)0.0;
00184     for (xcount = counter = 0; xcount <2; xcount++) {
00185         for (ycount = 0; ycount < 2; ycount++,counter++) {
00186         
00187             roi1_center_x = roi1_x1 + ((roi1_x2 - roi1_x1)/4) + (xcount*((roi1_x2 - roi1_x1)/2));
00188             roi1_center_y = roi1_y1 + ((roi1_y2 - roi1_y1)/4) + (ycount*((roi1_y2 - roi1_y1)/2));
00189             roi2_center_x = roi2_x1 + ((roi2_x2 - roi2_x1)/4) + (xcount*((roi2_x2 - roi2_x1)/2));
00190             roi2_center_y = roi2_y1 + ((roi2_y2 - roi2_y1)/4) + (ycount*((roi2_y2 - roi2_y1)/2));
00191             
00192             
00193             /* Compute the target denominator value */
00194             for (m=0; m<target_size; m++) {
00195                     for (n=0; n<target_size; n++) {
00196                                 denom_t += (img2->data[0][(roi2_center_y-((target_size-1)/2))+m][(roi2_center_x-((target_size-1)/2))+n] - mean_t)*(img2->data[0][(roi2_center_y-((target_size-1)/2))+m][(roi2_center_x-((target_size-1)/2))+n] - mean_t);
00197                     }
00198             }
00199 
00200         }
00201     }
00202     denom_t = (PIXEL) sqrt((double)denom_t);
00203 
00204 
00205     printf("mean_t = %f\n",mean_t);
00206     printf("mean_s = %f\n",mean_s);
00207     printf("denom_t = %f\n",denom_t);
00208 
00209     maxval = (PIXEL)0.0;
00210     max_i = max_j = 0;
00211     for (j=0; j<(search_size - target_size)+1; j++) {
00212         for (i=0; i<(search_size - target_size)+1; i++) {
00213                 denom_s = (PIXEL)0.0;
00214             num = (PIXEL)0.0;
00215                 
00216             for (xcount = counter = 0; xcount <2; xcount++) {
00217                 for (ycount = 0; ycount < 2; ycount++,counter++) {
00218                     roi1_center_x = roi1_x1 + ((roi1_x2 - roi1_x1)/4) + (xcount*((roi1_x2 - roi1_x1)/2));
00219                     roi1_center_y = roi1_y1 + ((roi1_y2 - roi1_y1)/4) + (ycount*((roi1_y2 - roi1_y1)/2));
00220                     roi2_center_x = roi2_x1 + ((roi2_x2 - roi2_x1)/4) + (xcount*((roi2_x2 - roi2_x1)/2));
00221                     roi2_center_y = roi2_y1 + ((roi2_y2 - roi2_y1)/4) + (ycount*((roi2_y2 - roi2_y1)/2));
00222 
00223                     for (m=0; m<target_size; m++) {
00224                                 for (n=0; n<target_size; n++) {
00225 
00226                                     temp_s = img1->data[0][(roi1_center_y-((search_size-1)/2))+m+j][(roi1_center_x-((search_size-1)/2))+n+i] - mean_s;
00227                                     temp_t = img2->data[0][(roi2_center_y-((target_size-1)/2))+m][(roi2_center_x-((target_size-1)/2))+n] - mean_t;
00228 
00229                                     denom_s += temp_s*temp_s;
00230                                     num += temp_t*temp_s;
00231                         }
00232                     }
00233                 }
00234                 }
00235             
00236             denom_s = (PIXEL) sqrt((double)denom_s);
00237 
00238             value = num / (denom_t * denom_s);
00239 
00240             (*corr)->data[0][j][i] = value;
00241 
00242             if (value > maxval) {
00243                 maxval = value;
00244                 max_i = i;
00245                 max_j = j;
00246             }
00247 
00248         }
00249     }
00250 
00251     xoffset = est_x_off + (max_i - ((search_size-target_size)/2));
00252     yoffset = est_y_off + (max_j - ((search_size-target_size)/2));
00253 
00254     printf("offset = (%d,%d)\n", xoffset, yoffset);
00255     printf("max = (%d,%d)\n", max_i, max_j);
00256     printf("deltamax = (%d,%d)\n", (max_i - ((search_size-target_size)/2)), (max_j - ((search_size-target_size)/2)));
00257     printf("max corr = %f\n\n", maxval);
00258 
00259 
00260 
00261     /* not time at beginning of actual registration process */
00262     gettimeofday(&doreg, NULL);
00263 
00264 
00265 
00266     FINDGAINADJ(img1, img2, xoffset, yoffset, &gain[1], &bias[1]);
00267     
00268     gain[0] = (PIXEL)1.0;
00269     bias[0] = (PIXEL)0.0;
00270     
00271     if (xoffset > 0) {
00272         kbgn[0] = 0;
00273         kbgn[1] = xoffset;
00274         npix = img1->npix + img2->npix - (img1->npix - xoffset);
00275     } else {
00276         kbgn[0] = 0-xoffset;
00277         kbgn[1] = 0;
00278         npix = img1->npix + img2->npix - (img2->npix + xoffset);
00279     }
00280 
00281     if (yoffset > 0) {
00282         jbgn[0] = 0;
00283         jbgn[1] = yoffset;
00284         nlin = img1->nlin + img2->nlin - (img1->nlin - yoffset);
00285     } else {
00286         jbgn[0] = 0-yoffset;
00287         jbgn[1] = 0;
00288         nlin = img1->nlin + img2->nlin - (img2->nlin + yoffset);
00289     }
00290     
00291     
00292     GAINADJMOS((IMAGE **)inimg, (int)2, (int)nlin, (int)npix, (int *)&jbgn[0], (int *)&kbgn[0], (PIXEL *)&bias[0], (PIXEL *)&gain[0], (PIXEL)0.0, (IMAGE **)out);
00293         
00294 
00295 
00296 
00297     /* note finishing time */
00298     gettimeofday(&end, NULL);
00299 
00300 
00301 
00302     the_end:
00303 
00304     //if (CHECKIMG(target))  RELIMG(&target);
00305     //if (CHECKIMG(search))  RELIMG(&search);
00306 
00307     if (NAMES) {
00308             MESSAGE('I',"");
00309         sprintf(msg," Computed x-offset: %d",xoffset);
00310             MESSAGE('I',msg);
00311         sprintf(msg," Computed y-offset: %d",yoffset);
00312             MESSAGE('I',msg);
00313             MESSAGE('I',"");
00314         sprintf(msg,"Time to compute maximum correlation = %ld ms.", delay(findreg, doreg));
00315             MESSAGE('I',msg);
00316         sprintf(msg,"Time to perform registration = %ld ms.", delay(doreg, end));
00317             MESSAGE('I',msg);
00318         sprintf(msg,"Total time required = %ld ms.", delay(start, end));
00319             MESSAGE('I',msg);
00320         MESSAGE('I'," ...............");
00321     }
00322 
00323     /* if (TIMES) TIMING(T_EXIT); */
00324 
00325 }
00326 
00327 
00328 
00329 
00330 /*
00331  * Compute the delay between t1 and t2 in milliseconds
00332  */
00333 long delay(struct timeval t1, struct timeval t2)
00334 {
00335         long d;
00336         
00337         d = (t2.tv_sec - t1.tv_sec) * 1000;
00338         d += ((t2.tv_usec - t1.tv_usec + 500) / 1000);
00339         return(d);
00340 }
00341 
00342 
00343 
00344