---

findgainadj.c

Go to the documentation of this file.

#include        "sadie.h"
#include        "proto.h"



/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function is based on stats.c from the SADIE library.  It takes two  */
/*   overlapping images and calculates necessary gain adjustments between the */
/*   two images.                                                              */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void FINDGAINADJ (
IMAGE  *imga,     /*  Pointer to the first input image.                       */
IMAGE  *imgb,     /*  Pointer to the second input image.                      */
int    xoffset,   /*  x-offset of second image                                */
int    yoffset,   /*  y-offset of second image                                */
PIXEL  *gain,     /*  Gain adjustment of second input image relative to       */
                  /*    the first input image                                 */
PIXEL  *bias      /*  Bias adjustment of second input image relative to       */
                  /*    the first input image                                 */
/*----------------------------------------------------------------------------*/
) { register int i, j, k;
    char   msg[SLEN];
    long   *hist=0, maxval, npts;
    double sum, sumsq;
    int leftlimita, rightlimita, toplimita, bottomlimita;
    int leftlimitb, rightlimitb, toplimitb, bottomlimitb;
    int overlapwidth, overlapheight;
    double gmina, gmaxa, meana, deva;
    double gminb, gmaxb, meanb, devb;


    /* check input */
    if (!CHECKIMG(imga)) {
        MESSAGE('E'," Can't identify image A.");
        goto the_end;
    } else if (!CHECKIMG(imgb)) {
        MESSAGE('E'," Can't identify image B.");
        goto the_end;
    }
    
    
    /* compute the overlap regions of each image */
    if (xoffset >= 0) {
            leftlimita = min(xoffset,imga->npix-1);
            rightlimita = min(imga->npix-1,(imgb->npix-1)+leftlimita);

            leftlimitb = 0;
            rightlimitb = min(imgb->npix-1,(imga->npix-1)-xoffset);
    } else {
            leftlimita = 0;
            rightlimita = min(imga->npix-1,(imgb->npix-1) + xoffset);

            leftlimitb = min(imgb->npix-1, 0-xoffset);
            rightlimitb = min(imgb->npix-1, (imga->npix-1) - xoffset);
    }

    if (yoffset >= 0) {
            toplimita = min(yoffset,imga->nlin-1);
            bottomlimita = min(imga->nlin-1,(imgb->nlin-1)+toplimita);

            toplimitb = 0;
            bottomlimitb = min(imgb->nlin-1, (imga->nlin-1)-yoffset);
    } else {
            toplimita = 0;
            bottomlimita = min(imga->nlin-1,(imgb->nlin-1) + yoffset);

            toplimitb = min(imgb->nlin-1, 0-yoffset);
            bottomlimitb = min(imgb->nlin-1, (imga->nlin-1) - yoffset);
    }   
    
    overlapheight = (bottomlimita - toplimita) + 1;
    overlapwidth = (rightlimita - leftlimita) + 1;
    
    /* compute the mean and deviation statistics for the overlap regions of each image */
    gmina = imga->data[0][toplimita][leftlimita];
    gmaxa = imga->data[0][toplimita][leftlimita];
    for (sum=sumsq=0.0,j=toplimita; j<=bottomlimita; j+=1) {
        for (k=leftlimita; k<=rightlimita; k+=1) {
            sum    += (double)imga->data[0][j][k];
            sumsq  += (double)(imga->data[0][j][k]*imga->data[0][j][k]);
            gmina = min(gmina,imga->data[0][j][k]);
            gmaxa = max(gmaxa,imga->data[0][j][k]);
        }
    }
    meana = sum/((double)(overlapheight+(overlapheight%1 ? 1:0))*(double)(overlapwidth+(overlapwidth%1 ? 1:0)));
    deva  = sqrt(sumsq/((double)(overlapheight+(overlapheight%1 ? 1:0))*(double)(overlapwidth+(overlapwidth%1 ? 1:0)))-meana*meana);

    gminb = imgb->data[0][toplimitb][leftlimitb];
    gmaxb = imgb->data[0][toplimitb][leftlimitb];
    for (sum=sumsq=0.0,j=toplimitb; j<=bottomlimitb; j+=1) {
        for (k=leftlimitb; k<=rightlimitb; k+=1) {
            sum    += (double)imgb->data[0][j][k];
            sumsq  += (double)(imgb->data[0][j][k]*imgb->data[0][j][k]);
            gminb = min(gminb,imgb->data[0][j][k]);
            gmaxb = max(gmaxb,imgb->data[0][j][k]);
        }
    }
    meanb = sum/((double)(overlapheight+(overlapheight%1 ? 1:0))*(double)(overlapwidth+(overlapwidth%1 ? 1:0)));
    devb  = sqrt(sumsq/((double)(overlapheight+(overlapheight%1 ? 1:0))*(double)(overlapwidth+(overlapwidth%1 ? 1:0)))-meanb*meanb);

    *gain = (PIXEL) (deva/devb);
    *bias = (PIXEL) (meana - (*gain * meanb));



    MESSAGE('I',"");
    MESSAGE('I',"----------");
    MESSAGE('I',"");
    MESSAGE('I',"FINDGAINADJ");
    MESSAGE('I',"");
    sprintf(msg," Input image #1:               %s",imga->text);
    MESSAGE('I',msg);
    MESSAGE('I',"   Minimum     Maximum       Mean     Deviation       ");
    sprintf(msg,"%12.4e%12.4e%12.4e%12.4e",(double)gmina,(double)gmaxa,meana,deva);
    MESSAGE('I',msg);
    sprintf(msg,"   Overlap Region: (%d,%d) to (%d,%d)",leftlimita+1,toplimita+1,rightlimita+1,bottomlimita+1);
    MESSAGE('I',msg);
    MESSAGE('I',"");
    sprintf(msg," Input image #2:               %s",imgb->text);
    MESSAGE('I',msg);
    MESSAGE('I',"   Minimum     Maximum       Mean     Deviation       ");
    sprintf(msg,"%12.4e%12.4e%12.4e%12.4e",(double)gminb,(double)gmaxb,meanb,devb);
    MESSAGE('I',msg);
    sprintf(msg,"   Overlap Region: (%d,%d) to (%d,%d)",leftlimitb+1,toplimitb+1,rightlimitb+1,bottomlimitb+1);
    MESSAGE('I',msg);
    MESSAGE('I',"");
    sprintf(msg," Gain Adjustment: %12.4e ,  Bias Adjustment: %12.4e",*gain,*bias);
    MESSAGE('I',msg);
    MESSAGE('I',"");
    MESSAGE('I'," ...............");

    the_end:
    MESSAGE('I',"");
}

---