combine.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This function computes a weighted linear combination of two images.      */
00009 /*                                                                            */
00010 /*----------------------------------------------------------------------------*/
00011 /*-Interface Information------------------------------------------------------*/
00012 void COMBINE (
00013 IMAGE  *in1,    /*  I   Pointer to the first input image.                     */
00014 IMAGE  *in2,    /*  I   Pointer to the second input image.                    */
00015 short  option,  /*  I   Switch, indicating type of combination operation:     */
00016                 /*      ADD      - Sum of the input images is computed.       */
00017                 /*      SUBTRACT - Difference of the input images is computed.*/
00018                 /*      MULTIPLY - Product of the input images is computed.   */
00019                 /*      DIVIDE   - Quotient of the input images is computed.  */
00020                 /*      MINIMUM  - Minimum of the input images is computed.   */
00021                 /*      MAXIMUM  - Maximum of the input images is computed.   */
00022 double wgt1,    /*  I   Weight on the first input image.                      */
00023 double wgt2,    /*  I   Weight on the second input image.                     */
00024 IMAGE  **out    /*  O   Address of a pointer to the output image.             */
00025                 /*      An input image may be overwritten (out = &in1,&in2).  */
00026 /*----------------------------------------------------------------------------*/
00027 ) { register short i, j, k;
00028     double factor;
00029     char msg[SLEN];
00030 
00031     if (TIMES) TIMING(T_COMBINE);
00032     if (NAMES) {
00033         MESSAGE('I',"");
00034         MESSAGE('I',"COMBINE");
00035         MESSAGE('I',"");
00036         sprintf(msg," First input image:        %s",in1->text);
00037         MESSAGE('I',msg);
00038         sprintf(msg," Second input image:       %s",in2->text);
00039         MESSAGE('I',msg);
00040         if (option == ADD) {
00041             MESSAGE('I'," Operation:                add");
00042         } else if (option == SUBTRACT) {
00043             MESSAGE('I'," Operation:                subtract");
00044         } else if (option == MULTIPLY) {
00045             MESSAGE('I'," Operation:                multiply");
00046         } else if (option == DIVIDE) {
00047             MESSAGE('I'," Operation:                divide");
00048         } else if (option == MINIMUM) {
00049             MESSAGE('I'," Operation:                minimum");
00050         } else if (option == MAXIMUM) {
00051             MESSAGE('I'," Operation:                maximum");
00052         }
00053         sprintf(msg," Weight on first image:  %12.4e",wgt1);
00054         MESSAGE('I',msg);
00055         sprintf(msg," Weight on second image: %12.4e",wgt2);
00056         MESSAGE('I',msg);
00057         MESSAGE('I'," ...............");
00058     }
00059 
00060     /* check input */
00061     if (!CHECKIMG(in1)) {
00062         MESSAGE('E'," Can't identify first image.");
00063         goto the_end;
00064     } else if (!CHECKIMG(in2)) {
00065         MESSAGE('E'," Can't identify second image.");
00066         goto the_end;
00067     } else if (in1->nbnd != in2->nbnd) {
00068         MESSAGE('E'," Number of bands must be identical in both images.");
00069         goto the_end;
00070     } else if (in1->nlin != in2->nlin) {
00071         MESSAGE('E'," Number of lines must be identical in both images.");
00072         goto the_end;
00073     } else if (in1->npix != in2->npix) {
00074         MESSAGE('E'," Number of pixels/line must be identical in both images.");
00075         goto the_end;
00076     } else if (option != ADD  &&  option != SUBTRACT  &&  option != MULTIPLY  &&  option != DIVIDE  &&  option != MINIMUM  &&  option != MAXIMUM) {
00077         MESSAGE('E'," Option must be one of ADD, SUBTRACT, MULTIPLY, DIVIDE, MINIMUM, or MAXIMUM.");
00078         goto the_end;
00079     } else if (option == DIVIDE  &&  wgt2 == 0.0) {
00080         MESSAGE('E'," Denominator weight must not be equal to zero.");
00081         goto the_end;
00082     }
00083 
00084     /* create image of appropriate size */
00085     if (!CHECKIMG(*out)) GETMEM(in1->nbnd,in1->nlin,in1->npix,out);
00086     if (!*out) goto the_end;
00087 
00088     switch (option) {
00089     case ADD:
00090         /* add two images */
00091         for (i=0; i<in1->nbnd; i++) {
00092             for (j=0; j<in1->nlin; j++) {
00093                 for (k=0; k<in1->npix; k++) {
00094                     (*out)->data[i][j][k] = (PIXEL)((double)in1->data[i][j][k]*wgt1+(double)in2->data[i][j][k]*wgt2);
00095                 }
00096             }
00097         }
00098         break;
00099     case SUBTRACT:
00100         /* subtract two images */
00101         for (i=0; i<in1->nbnd; i++) {
00102             for (j=0; j<in1->nlin; j++) {
00103                 for (k=0; k<in1->npix; k++) {
00104                     (*out)->data[i][j][k] = (PIXEL)((double)in1->data[i][j][k]*wgt1-(double)in2->data[i][j][k]*wgt2);
00105                 }
00106             }
00107         }
00108         break;
00109     case MULTIPLY:
00110         /* multiply two images */
00111         factor = wgt1*wgt2;
00112         for (i=0; i<in1->nbnd; i++) {
00113             for (j=0; j<in1->nlin; j++) {
00114                 for (k=0; k<in1->npix; k++) {
00115                     (*out)->data[i][j][k] = (PIXEL)((double)in1->data[i][j][k]*(double)in2->data[i][j][k]*factor);
00116                 }
00117             }
00118         }
00119         break;
00120     case DIVIDE:
00121         /* divide two images */
00122         factor = wgt1/wgt2;
00123         for (i=0; i<in1->nbnd; i++) {
00124             for (j=0; j<in1->nlin; j++) {
00125                 for (k=0; k<in1->npix; k++) {
00126                     if (in2->data[i][j][k] == (PIXEL)0) {
00127                         (*out)->data[i][j][k] = k == 0 ? (PIXEL)0 : (*out)->data[i][j][k-1];
00128                     } else {
00129                         (*out)->data[i][j][k] = (PIXEL)((double)in1->data[i][j][k]/(double)in2->data[i][j][k]*factor);
00130                     }  
00131                 }
00132             }
00133         }
00134         break;
00135     case MINIMUM:
00136         /* find minimum of two images */
00137         for (i=0; i<in1->nbnd; i++) {
00138             for (j=0; j<in1->nlin; j++) {
00139                 for (k=0; k<in1->npix; k++) {
00140                     (*out)->data[i][j][k] = (PIXEL)min((double)in1->data[i][j][k]*wgt1,(double)in2->data[i][j][k]*wgt2);
00141                 }
00142             }
00143         }
00144         break; 
00145     case MAXIMUM:
00146         /* find maximum of two images */
00147         for (i=0; i<in1->nbnd; i++) {
00148             for (j=0; j<in1->nlin; j++) {
00149                 for (k=0; k<in1->npix; k++) {
00150                     (*out)->data[i][j][k] = (PIXEL)max((double)in1->data[i][j][k]*wgt1,(double)in2->data[i][j][k]*wgt2);
00151                 }
00152             }
00153         }
00154         break;
00155     }
00156 
00157     the_end:
00158     if (TIMES) TIMING(T_EXIT);
00159 }