vstretch.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This function performs a spatially-variable graylevel transformation     */
00009 /*   on a single-band image.                                                  */
00010 /*                                                                            */
00011 /*----------------------------------------------------------------------------*/
00012 /*-Background Information-----------------------------------------------------*/
00013 /*                                                                            */
00014 /*   Fahnestock, J.D. and Schowengerdt, R.A.:                                 */
00015 /*   "Spatially variant contrast enhancement using local range modification." */
00016 /*   Optical Engineering, Vol. 22, No. 3, May/June 1983                       */
00017 /*                                                                            */
00018 /*----------------------------------------------------------------------------*/
00019 /*-Interface Information------------------------------------------------------*/
00020 void VSTRETCH (
00021 IMAGE *in,      /*  I   Pointer to the input image.                           */
00022 short size,     /*  I   Number of lines and pixels in image blocks.           */
00023 IMAGE **out     /*  O   Address of a pointer to the output image.             */
00024 /*----------------------------------------------------------------------------*/
00025 ) { register short j, k, u, v;
00026     char   msg[SLEN];
00027     long   total, x, y, idx1, idx2, blkx, blky;
00028     double rat1, rat2, rat3, rat4, pinc, psum;
00029     PIXEL  *hi=0, *lo=0, *mp=0, *np=0, m, n;
00030 
00031     if (TIMES) TIMING(T_VSTRETCH);
00032     if (NAMES) {
00033         MESSAGE('I',"");
00034         MESSAGE('I',"VSTRETCH");
00035         MESSAGE('I',"");
00036         sprintf(msg," Input image:                    %s",in->text);
00037         MESSAGE('I',msg);
00038         sprintf(msg," Block size: lines and pixels:   %d",size);
00039         MESSAGE('I',msg);
00040         MESSAGE('I'," ...............");
00041     }
00042 
00043     /* check input */
00044     if (!CHECKIMG(in)) {
00045         MESSAGE('E'," Can't identify image.");
00046         goto the_end;
00047     } else if (in->nbnd > 1) {
00048         MESSAGE('E'," Image must be single-band.");
00049         goto the_end;
00050     } else if (size < 1) {
00051         MESSAGE('E'," Size of subblock must be equal to or greater than one.");
00052         goto the_end;
00053     } else if (size > in->nlin  ||  size > in->npix) {
00054         MESSAGE('E'," Image size must be equal to or greater than subblock size.");
00055         goto the_end;
00056     }
00057 
00058     /* create image of appropriate size */
00059     if (!CHECKIMG(*out)) GETMEM(1,in->nlin,in->npix,out);
00060     if (!*out) goto the_end;
00061     blkx = (in->npix-1)/size + 1L;
00062     blky = (in->nlin-1)/size + 1L;
00063 
00064     /* allocate space for arrays */
00065     total = blkx * blky;
00066     if (!(hi=(PIXEL *)malloc(total*sizeof(PIXEL)))) {
00067         MESSAGE('E'," Memory request failed.");
00068         goto the_end;
00069     }
00070     if (!(lo=(PIXEL *)malloc(total*sizeof(PIXEL)))) {
00071         MESSAGE('E'," Memory request failed.");
00072         goto the_end;
00073     }
00074     total = (blkx+1L) * (blky+1L);
00075     if (!(mp=(PIXEL *)malloc(total*sizeof(PIXEL)))) {
00076         MESSAGE('E'," Memory request failed.");
00077         goto the_end;
00078     }
00079     if (!(np=(PIXEL *)malloc(total*sizeof(PIXEL)))) {
00080         MESSAGE('E'," Memory request failed.");
00081         goto the_end;
00082     }
00083 
00084     /* initialize progress indicator */
00085     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00086     pinc = 0.5/(double)blkx/(double)blky;
00087 
00088     /* first pass: obtain transformation parameters from subblocks */
00089     /* search for maximum and minimum */
00090     for (y=0L; y<blky; y++) {
00091         for (x=0L; x<blkx; x++) {
00092             for (m=n=in->data[0][y*size][x*size],j=y*size; j<(y+1L)*size && j<in->nlin; j++) {
00093                 for (k=x*size; k<(x+1L)*size && k<in->npix; k++) {
00094                     if (in->data[0][j][k] > m) m = in->data[0][j][k];
00095                     if (in->data[0][j][k] < n) n = in->data[0][j][k];
00096                 }
00097             }
00098             hi[y*blkx+x] = m;
00099             lo[y*blkx+x] = n;
00100             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00101         }
00102     }
00103     
00104     /* obtain values for parameter arrays */
00105     /* four corners */
00106     mp[0] = hi[0];
00107     np[0] = lo[0];
00108     mp[blkx] = hi[blkx-1L];
00109     np[blkx] = lo[blkx-1L];
00110     mp[blky*(blkx+1L)] = hi[(blky-1L)*blkx];
00111     np[blky*(blkx+1L)] = lo[(blky-1L)*blkx];
00112     mp[(blkx+1L)*(blky+1L)-1L] = hi[blkx*blky-1L];
00113     np[(blkx+1L)*(blky+1L)-1L] = lo[blkx*blky-1L];
00114     
00115     /* four sides */
00116     for (y=1L; y<blky; y++) {
00117         idx1 = (y-1L)*blkx;
00118         idx2 = y*blkx;
00119         mp[y*(blkx+1L)] = max(hi[idx1],hi[idx2]);
00120         np[y*(blkx+1L)] = min(lo[idx1],lo[idx2]);
00121         idx1 = (y+1L)*blkx - 1L;
00122         mp[(y+1L)*(blkx+1L)-1L] = max(hi[idx2-1L],hi[idx1]);
00123         np[(y+1L)*(blkx+1L)-1L] = min(lo[idx2-1L],lo[idx1]);
00124     }
00125     for (x=1L; x<blkx; x++) {
00126         mp[x] = max(hi[x-1L],hi[x]);
00127         np[x] = min(lo[x-1L],lo[x]);
00128         idx1 = blky*(blkx+1L)+x;
00129         idx2 = (blky-1L)*blkx+x;
00130         mp[idx1] = max(hi[idx2-1L],hi[idx2]);
00131         np[idx1] = min(lo[idx2-1L],lo[idx2]);
00132     }
00133     
00134     /* inside */
00135     for (y=1L; y<blky; y++) {
00136         for (x=1L; x<blkx; x++) {
00137             idx1 = (y-1L)*blkx + x;
00138             idx2 = y*blkx + x;
00139             mp[y*(blkx+1L)+x] = max(max(hi[idx1],hi[idx1-1L]),max(hi[idx2],hi[idx2-1L]));
00140             np[y*(blkx+1L)+x] = min(min(lo[idx1],lo[idx1-1L]),min(lo[idx2],lo[idx2-1L]));
00141         }
00142     }
00143 
00144     /* second pass: perform transformation */
00145     for (x=0L; x<blkx; x++) {
00146         for (y=0L; y<blky; y++) {
00147             for (u=0,j=y*size; j<(y+1L)*size && j<in->nlin; u++,j++) {
00148                 for (v=0,k=x*size; k<(x+1L)*size && k<in->npix; v++,k++) {
00149                     idx1 = (y+1L)*(blkx+1L)+x;
00150                     idx2 = y*(blkx+1L)+x;
00151                     rat1 = (double)u/(double)size;
00152                     rat2 = (double)(size-u)/(double)size;
00153                     rat3 = (double)(size-v)/(double)size;
00154                     rat4 = (double)v/(double)size;
00155                     m = (rat1*(double)mp[idx1]+rat2*(double)mp[idx2])*rat3 + (rat1*(double)mp[idx1+1L]+rat2*(double)mp[idx2+1L])*rat4;
00156                     n = (rat1*(double)np[idx1]+rat2*(double)np[idx2])*rat3 + (rat1*(double)np[idx1+1L]+rat2*(double)np[idx2+1L])*rat4;
00157                     (*out)->data[0][j][k] = m == n ? in->data[0][j][k] : (PIXEL)((double)(in->gmax-in->gmin)/(double)(m-n)*((double)in->data[0][j][k]-n)) + in->gmin;
00158                 }
00159             }
00160             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00161         }
00162     }
00163 
00164     the_end:
00165     if (hi)    free(hi);
00166     if (lo)    free(lo);
00167     if (mp)    free(mp);
00168     if (np)    free(np);
00169     if (TIMES) TIMING(T_EXIT);
00170 }