destripe.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1992 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   There are three methods to use in destriping an image.                   */
00009 /*   All use either a single detector or the average of all detectors as      */
00010 /*   a reference.  The first method adjusts the means of all the detectors    */
00011 /*   to match that of the reference.  The second method adjusts the means and */
00012 /*   standard deviations to match those of the reference.  The third method   */
00013 /*   computes cumulative distribution functions using a specified number of   */
00014 /*   levels for the reference and each of the detectors.  These are then used */
00015 /*   in matching the contrast of the detectors to that of the reference.      */
00016 /*                                                                            */
00017 /*----------------------------------------------------------------------------*/
00018 /*-Interface Information------------------------------------------------------*/
00019 void DESTRIPE (
00020 IMAGE *in,       /*  I  Pointer to the image.                                 */
00021 short option,    /*  I  Switch, indicating type of line matching              */
00022                  /*     DESMNE - Destripe using mean only levels              */
00023                  /*     DESSTD - Destripe using mean and standard dev.        */
00024                  /*     DESCDF - Destripe using cumulative dist. function     */
00025 short roption,   /*  I  Switch, indicating type of reference used             */
00026                  /*     DESSNG - Use a single line (detector) as reference    */
00027                  /*     DESAVG - Average all detectors for reference          */
00028 short lines,     /*  I  Number of lines to destripe - ie, number of           */
00029                  /*     detectors.  Must be a factor of total image line num  */
00030 short refline,   /*  I  The line (detector) to be used for destriping.        */
00031                  /*     This value does not matter if DESAVG is selected.     */
00032 short nlev,      /*  I  The number of levels used in the cumulative           */
00033                  /*     distribution function for the reference and           */
00034                  /*     transformation lines.                                 */
00035 IMAGE **out      /*  O  Address of a pointer to the output image.             */
00036                  /*     The input image may be overwritten (out = &in1).      */
00037 /*----------------------------------------------------------------------------*/
00038 ) { register short i, j, k, l;
00039     short  skip, index, index2;
00040     double mean, stddev, oldmean, oldstddev, factor, *rcdf=0, *tcdf=0;
00041     PIXEL  *table=0;
00042 
00043     if (TIMES) TIMING(T_DESTRIPE);
00044     if (NAMES) {
00045         MESSAGE('I',"");
00046         MESSAGE('I',"DESTRIPE");
00047     }
00048 
00049     /* check input */
00050     if (!CHECKIMG(in)) {
00051         MESSAGE('E'," Can't identify image.");
00052         goto the_end;
00053     } else if (option != DESMNE  &&  option != DESSTD  &&  option != DESCDF) {
00054         MESSAGE('E'," Option must be DESMNE, DESSTD, or DESCDF.");
00055         goto the_end;
00056     } else if (roption != DESSNG  &&  roption != DESAVG) {
00057         MESSAGE('E'," Reference option must be either DESSNG or DESAVG.");
00058         goto the_end;
00059     }
00060 
00061     /* create image of appropriate size */
00062     if (!CHECKIMG(*out)) GETMEM(in->nbnd,in->nlin,in->npix,out);
00063     if (!*out) goto the_end;
00064 
00065     /* Set up image hopping parameters */
00066     if (roption == DESAVG) {
00067         refline = 1;
00068         skip = 1;
00069     } else {
00070         skip = lines;
00071     }
00072 
00073     if (option == DESMNE  ||  option == DESSTD) {
00074 
00075         for (i=0; i<in->nbnd; i++) {
00076 
00077             /* compute reference mean */
00078             mean = 0.0;
00079             for (j=refline-1; j<in->nlin; j+=skip) {
00080                 for (k=0; k<in->npix; k++) {
00081                     mean += (double)in->data[i][j][k];
00082                 }
00083             }
00084             if (roption == DESSNG) {
00085                 mean /= (double)in->npix * (double)(j-refline+1) / (double)lines;
00086             } else {
00087                 mean /= (double)in->npix * (double)in->nlin;
00088             }
00089 
00090             /* compute reference standard deviation */
00091             if (option == DESSTD) {
00092                 stddev = 0.0;
00093                 for (j=refline-1; j<in->nlin; j+=skip) {
00094                     for (k=0; k<in->npix; k++) {
00095                         stddev += ((double)in->data[i][j][k]-mean) * ((double)in->data[i][j][k]-mean);
00096                     }
00097                 }
00098                 if (roption == DESSNG) {
00099                     stddev /= (double)in->npix * (double)(j-refline+1) / (double)lines;
00100                 } else {
00101                     stddev /= (double)in->npix * (double)in->nlin;
00102                 }
00103                 stddev = sqrt(stddev);
00104             }
00105 
00106             /* Calculate mean and/or standard deviation of each line to be transformed. Copy reference line to output. */
00107             for (j=0; j<lines; j++) {
00108 
00109                 if (roption == DESSNG  &&  j == refline-1) {
00110                     /* copy reference lines (detector) directly to output image */
00111                     for (k=j; k<in->nlin; k+=lines) {
00112                         for (l=0; l<in->npix; l++) {
00113                             (*out)->data[i][k][l] = in->data[i][k][l];
00114                         }
00115                     }
00116                 } else {
00117                     /* compute mean of one set of lines (one detector) */
00118                     oldmean = 0.0;
00119                     for (k=j; k<in->nlin; k+=lines) {
00120                         for (l=0; l<in->npix; l++) {
00121                             oldmean += (double)in->data[i][k][l];
00122                         }
00123                     }
00124                     oldmean /= (double)in->npix * (double)(k-j) / (double)lines;
00125 
00126                     /* compute standard deviation of one set of lines (one detector) */
00127                     if (option == DESSTD) {
00128                         oldstddev = 0.0;
00129                         for (k=j; k<in->nlin; k+=lines) {
00130                             for (l=0; l<in->npix; l++) {
00131                                 oldstddev += ((double)in->data[i][k][l]-oldmean) * ((double)in->data[i][k][l]-oldmean);
00132                             }
00133                         }
00134                         oldstddev /= (double)in->npix * (double)(k-j) / (double)lines;
00135                         oldstddev = sqrt(oldstddev);
00136                     }
00137 
00138                     /* copy transformed line to output, adjusting gray levels to achieve mean and/or standard deviation of reference */
00139                     for (k=j; k<in->nlin; k+=lines) {
00140                         for (l=0; l<in->npix; l++) {
00141                             if (option == DESMNE) {
00142                                 (*out)->data[i][k][l] = (PIXEL)((double)in->data[i][k][l] - oldmean + mean);
00143                             } else {
00144                                 (*out)->data[i][k][l] = (PIXEL)(((double)in->data[i][k][l]-oldmean) * stddev / oldstddev + mean);
00145                             }
00146                         }
00147                     }
00148                 }
00149             }
00150         }
00151 
00152     } else {
00153 
00154         if (!(rcdf=(double*)malloc(nlev*sizeof(double)))) {
00155             MESSAGE ('E'," Memory request failed.");
00156             goto the_end;
00157         }
00158         if (!(tcdf=(double*)malloc(nlev*sizeof(double)))) {
00159             MESSAGE ('E'," Memory request falied.");
00160             goto the_end;
00161         }
00162         if (!(table=(PIXEL*)malloc(nlev*sizeof(PIXEL)))) {
00163             MESSAGE ('E'," Memory request failed.");
00164             goto the_end;
00165         }
00166 
00167         RANGE(in);
00168         factor = (double)(nlev-1)/(double)(in->gmax-in->gmin);
00169 
00170         for (i=0; i<in->nbnd; i++) {
00171 
00172             /* compute histogram of reference */
00173             for (j=0; j<nlev; j++) {
00174                 rcdf[j] = 0.0;
00175             }
00176             for (j=refline-1; j<in->nlin; j+=skip) {
00177                 for (k=0; k<in->npix; k++) {
00178                     rcdf[(short)((double)(in->data[i][j][k]-in->gmin)*factor + 0.5)] += 1.0;
00179                 }
00180             }
00181 
00182             /* convert histogram to cumulative distribution function and scale to [0,1] */
00183             for (j=1; j<nlev; j++) {
00184                 rcdf[j] += rcdf[j-1];
00185             }
00186             for (j=0; j<nlev; j++) {
00187                 rcdf[j] /= rcdf[nlev-1];
00188             }
00189 
00190             /* loop through each line, calculate cdf, and use it to transform the line */
00191             for (j=0; j<lines; j++) {
00192 
00193                 /* copy reference lines (detector) directly to output image */
00194                 if (roption == DESSNG  &&  j == refline-1) {
00195                     for (k=j; k<in->nlin; k+=lines) {
00196                         for (l=0; l<in->npix; l++) {
00197                             (*out)->data[i][k][l] = in->data[i][k][l];
00198                         }
00199                     }
00200                 } else {
00201 
00202                     /* compute histogram of one set of lines (one detector) */
00203                     for (k=0; k<nlev; k++) {
00204                         tcdf[k] = 0.0;
00205                     }
00206                     for (k=j; k<in->nlin; k+=lines) {
00207                         for (l=0; l<in->npix; l++) {
00208                             tcdf[(short)((double)(in->data[i][k][l]-in->gmin)*factor + 0.5)] += 1;
00209                         }
00210                     }
00211 
00212                     /* convert histogram to cumulative distribution function amd scale to [0,1] */
00213                     for (k=1; k<nlev; k++) {
00214                         tcdf[k] += tcdf[k-1];
00215                     }
00216                     for (k=0; k<nlev; k++) {
00217                         tcdf[k] /= tcdf[nlev-1];
00218                     }
00219 
00220                     /* loop through histograms making lookup table entries */
00221                     for (index=0; index<nlev; index++) {
00222                         for (index2=0; index2<nlev; index2++) {
00223                             if (rcdf[index2] == tcdf[index]) {
00224                                 table[index] = (PIXEL)index2/(PIXEL)factor + in->gmin;
00225                                 break;
00226                             } else if (rcdf[index2] > tcdf[index]) {
00227                                 if (index2 == 0) {
00228                                     table[index] = in->gmin;
00229                                 } else {
00230                                     if (rcdf[index2]-tcdf[index] <= tcdf[index]-rcdf[index2 - 1]) {
00231                                         table[index] = (PIXEL)index2/(PIXEL)factor + in->gmin;
00232                                     } else {
00233                                         table[index] = (PIXEL)(index2-1)/(PIXEL)factor + in->gmin;
00234                                     }
00235                                 }
00236                                 break;
00237                             }
00238                         }
00239                     }
00240 
00241                     /* create new set of lines (detector) for output image */
00242                     for (k=j; k<in->nlin; k+=lines) {
00243                         for (l=0; l<in->npix; l++) {
00244                             index = (short)((double)(in->data[i][k][l]-in->gmin)*factor + 0.5);
00245                             (*out)->data[i][k][l] = table[index];
00246                         }
00247                     }
00248                 }
00249             }
00250         }
00251     }
00252 
00253     the_end:
00254     if (rcdf)  free(rcdf);
00255     if (tcdf)  free(tcdf);
00256     if (table) free(table);
00257     if (TIMES) TIMING(T_EXIT);
00258 }