/* @(#)mo_match.c 17.1.1.1 (ES0-DMD) 01/25/02 17:49:06 */ /*=========================================================================== Copyright (C) 1995 European Southern Observatory (ESO) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, MA 02139, USA. Corresponding concerning ESO-MIDAS should be addressed as follows: Internet e-mail: midas@eso.org Postal address: European Southern Observatory Data Management Division Karl-Schwarzschild-Strasse 2 D 85748 Garching bei Muenchen GERMANY ===========================================================================*/ /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ .COPYRIGHT (c) 1995 European Southern Observatory .IDENTIFIER mo_match2d.c .AUTHOR R.H. Warmels IPG-ESO Garching .KEYWORDS intensity matching software .LANGUAGE C .PURPOSE Routine to compute the intensity matching parameters .ENVIRONment MIDAS #include Symbols used by the ccd package .VERSION 1.0 16-May-1995 creation ------------------------------------------------------------*/ /* * Define _POSIX_SOURCE to indicate * that this is a POSIX program */ #define _POSIX_SOURCE 1 /* * definition of the used functions in this module */ #include #include #include #include #include /* MO_OVERLAP -- Routine to compute the intensity matching paremeters */ MO_OVERLAP(pc1out, pc2out, pl1out, pl2out, c1out, c2out, l1out, l2out, oc1out, oc2out, ol1out, ol2out) int pc1out, pc2out; /* previous subraster column limits */ int pl1out, pl2out; /* previous subraster line limits */ int c1out, c2out; /* current subraster column limits */ int l1out, l2out; /* current subraster line limits */ int *oc1out, *oc2out; /* overlap column limits */ int *ol1out, *ol2out; /* overlap line limits */ { /* Check for the case where no intersection is present. */ if (c1out > pc2out || c2out < pc1out || l1out > pl2out || l2out < pl1out) return 0; /* Compute the column overlap limits. */ if (pc1out <= c1out) *oc1out = c1out; else *oc1out = pc1out; if (pc2out <= c2out) *oc2out = pc2out; else *oc2out = c2out; /* Compute the line overlap limits. */ if (pl1out <= l1out) *ol1out = l1out; else *ol1out = pl1out; if (pl2out <= l2out) *ol2out = pl2out; else *ol2out = l2out; return 1; } /* MO_MED -- Routine to compute the median value of an array */ MO_MED(n, a, amed) int n; float *a; float *amed; { int i, imed; float denom; float *w; w = (float *) osmmget(n * sizeof( float )); for (i = 0; i < n; i++) w[i] = a[i]; (void) sortr(n,w); imed = n/2 + 1; denom = 2.0; if (fmod((float) n, denom) == 0) *amed = .5 * (w[imed] + w[imed-1]); else *amed = w[imed]; osmmfree((char *) w); } /* MO_M2MATCH -- Routine to compute the intensity matching paremeters */ MO_M2MATCH(inim, ranges, nranges, ic1, ic2, il1, il2, deltax, deltay, deltai) int inim; int *ranges; int nranges; int *ic1; int *ic2; int *il1; int *il2; float *deltax; float *deltay; float *deltai; { /* Initialize the intensity subraster. */ MO_VECINIT(deltai, MO_NXSUB * MO_NYSUB, ranges, nranges); if (*ranges == (int) MO_NULL) return; /* Match the intensities in the direction of observation. */ MO_OMATCH(inim, ic1, ic2, il1, il2, deltax, deltay, deltai); /* Match the intensities in the other direction. */ MO_NMATCH(inim, ic1, ic2, il1, il2, deltax, deltay, deltai); } /* MO_VECINIT -- Procedure to initialize the intensity matching algorithm. If the ranges are undefined and no matching is to take place the ishifts are set to INDEFR and the routine returns. Otherwise the shifts are all initialized to zero and shifts for the missing subrasters are set to INDEFR. */ MO_VECINIT(deltai, nsubrasters, ranges, nranges) float *deltai; /* intensity shifts */ int nsubrasters; /* number of subrasters */ int *ranges; /* ranges of missing subrasters */ int nranges; /* total number matching subrasters */ { int i; /* Initialize the shifts to INDEFR. */ for (i = 0; i < nsubrasters; i++) deltai[i] = MO_INDEFR; if (*ranges == (int) MO_NULL) return; for (i = 0; i < nsubrasters; i++) deltai[i] = 0.0; } /* MO_OMATCH -- Procedure to match images in the direction of observation direction. */ MO_OMATCH(im, ic1, ic2, il1, il2, deltax, deltay, deltai) int im; /* pointer to the input image */ int *ic1; /* beginning column limits */ int *ic2; /* ending column limits */ int *il1; /* beginning line limits */ int *il2; /* ending line limits */ float *deltax; /* array of x shifts */ float *deltay; /* array of y shifts */ float *deltai; /* array of intensity shifts */ { int num, nimages, nrasters; int stat; int npixi[3]; int uni, nulo; int pc1, pc2, pl1, pl2; int c1, c2, l1, l2; int pideltax, pideltay, ideltax, ideltay; int oc1, oc2, ol1, ol2; float pmedian, median, dif; float rem, *buf; int iav, ism, size; float image[4]; /* Get the number of pixels */ stat = SCDRDI(im,"NPIX",1,3,&iav,npixi,&uni,&nulo); ism = 0; size = npixi[0] * npixi[1]; buf = (float *) osmmget( size * sizeof( float )); /* Compute the do loop parameters. */ nimages = MO_NXSUB * MO_NYSUB; if (strcmp(MO_ORDER,MO_ROW) == 0) nrasters = MO_NXSUB; else nrasters = MO_NYSUB; /* Loop over the subrasters to be matched. */ for (num = 1; num <= nimages; num++) { rem = fmod ((float) num, (float) nrasters); if (rem > 0 && rem < 2) { /* Get the position and shift for the first subraster in the column */ pideltax = NINT (deltax[num-1]); pideltay = NINT (deltay[num-1]); pc1 = ic1[num-1]; pc2 = ic2[num-1]; pl1 = il1[num-1]; pl2 = il2[num-1]; num = num + 1; dif = 0.0; /* Get the position and shift for the next subraster in the column */ ideltax = NINT (deltax[num-1]); ideltay = NINT (deltay[num-1]); c1 = ic1[num-1]; c2 = ic2[num-1]; l1 = il1[num-1]; l2 = il2[num-1]; } else { /* Reset the coordinates of the previous subraster. */ pc1 = c1; pc2 = c2; pl1 = l1; pl2 = l2; pideltax = ideltax; pideltay = ideltay; /* Get the positions and shifts of the next subraster. */ ideltax = NINT (deltax[num-1]); ideltay = NINT (deltay[num-1]); c1 = ic1[num-1]; c2 = ic2[num-1]; l1 = il1[num-1]; l2 = il2[num-1]; } /* Compute the overlap region. */ stat = MO_OVERLAP(pc1 + pideltax, pc2 + pideltax, pl1 + pideltay, pl2 + pideltay, c1 + ideltax, c2 + ideltax, l1 + ideltay, l2 + ideltay, &oc1, &oc2, &ol1, &ol2); if (stat == 1) { image[0] = (float) MYMAX (pc1, MYMIN (oc1 - pideltax, pc2)); image[1] = (float) MYMIN (pc2, MYMAX (oc2 - pideltax, pc1)); image[2] = (float) MYMAX (pl1, MYMIN (ol1 - pideltay, pl2)); image[3] = (float) MYMIN (pl2, MYMAX (ol2 - pideltay, pl1)); GETDAT(im, size, npixi, image, ism, buf); MO_MED((int) (image[1] - image[0] + 1) * (int) (image[3] - image[2] + 1), buf, &pmedian); image[0] = (float) MYMAX (c1, MYMIN (oc1 - ideltax, c2)); image[1] = (float) MYMIN (c2, MYMAX (oc2 - ideltax, c1)); image[2] = (float) MYMAX (l1, MYMIN (ol1 - ideltay, l2)); image[3] = (float) MYMIN (l2, MYMAX (ol2 - ideltay, l1)); GETDAT (im, size, npixi, image, ism, buf); MO_MED((int) (image[1] - image[0] + 1) * (int) (image[3] - image[2] + 1), buf, &median); dif = dif + median - pmedian; if ( deltai[num-1] != MO_INDEFR) deltai[num-1] = deltai[num-1] - dif; } } osmmfree((char *) buf); } /* MO_NMATCH -- Procedure to match images in the direction of observation direction. */ MO_NMATCH(im, ic1, ic2, il1, il2, deltax, deltay, deltai) int im; /* pointer to the input image */ int *ic1; /* beginning column limits */ int *ic2; /* ending column limits */ int *il1; /* beginning line limits */ int *il2; /* ending line limits */ float *deltax; /* array of x shifts */ float *deltay; /* array of y shifts */ float *deltai; /* array of intensity shifts */ { int num, nimages, nrasters, fac; int stat; int npixi[3]; int uni, nulo; int pc1, pc2, pl1, pl2; int c1, c2, l1, l2; int pideltax, pideltay, ideltax, ideltay; int oc1, oc2, ol1, ol2; float pmedian, median, dif; float pdif, tdif; float *buf; int count; int iav, ism; float image[4]; int size; /* Get the number of pixels */ stat = SCDRDI(im,"NPIX",1,3,&iav,npixi,&uni,&nulo); ism = 0; size = npixi[0] * npixi[1]; buf = (float *) osmmget( size * sizeof( float )); /* Compute the do loop parameters. */ nimages = MO_NXSUB * MO_NYSUB; if (strcmp(MO_ORDER,MO_ROW) == 0) nrasters = MO_NXSUB; else nrasters = MO_NYSUB; fac = 2 * nrasters; /* Loop over the subrasters to be matched. */ num = 1; count = 1; do { /* Get the position and shift for the first subraster in the column */ if (num < nrasters) { pideltax = NINT (deltax[num-1]); pideltay = NINT (deltay[num-1]); pc1 = ic1[num-1]; pc2 = ic2[num-1]; pl1 = il1[num-1]; pl2 = il2[num-1]; if (deltai[num-1] == MO_INDEFR) pdif= 0.0; else pdif = deltai[num-1]; tdif = 0.0; if (MO_RASTER == "YES") { num = fac - num + 1; fac = fac + fac; } else num = num + nrasters; /* Get the position and shift for the next */ ideltax = NINT (deltax[num-1]); ideltay = NINT (deltay[num-1]); c1 = ic1[num-1]; c2 = ic2[num-1]; l1 = il1[num-1]; l2 = il2[num-1]; if (deltai[num-1] == MO_INDEFR) dif = 0.0; else dif = deltai[num-1]; } else { /* Reset the coordinates of the previous subraster. */ pc1 = c1; pc2 = c2; pl1 = l1; pl2 = l2; pideltax = ideltax; pideltay = ideltay; pdif = dif; /* Get the positions and shifts of the next subraster. */ ideltax = NINT (deltax[num-1]); ideltay = NINT (deltay[num-1]); c1 = ic1[num-1]; c2 = ic2[num-1]; l1 = il1[num-1]; l2 = il2[num-1]; if (deltai[num-1] == MO_INDEFR) dif = 0.0; else dif = deltai[num-1]; } /* Compute the overlap region. */ stat = MO_OVERLAP(pc1 + pideltax, pc2 + pideltax, pl1 + pideltay, pl2 + pideltay, c1 + ideltax, c2 + ideltax, l1 + ideltay, l2 + ideltay, &oc1, &oc2, &ol1, &ol2); if (stat == 1) { image[0] = (float) MYMAX (pc1, MYMIN (oc1 - pideltax, pc2)); image[1] = (float) MYMIN (pc2, MYMAX (oc2 - pideltax, pc1)); image[2] = (float) MYMAX (pl1, MYMIN (ol1 - pideltay, pl2)); image[3] = (float) MYMIN (pl2, MYMAX (ol2 - pideltay, pl1)); GETDAT(im, size, npixi, image, ism, buf); GETDAT(im, size, npixi, image, ism, buf); MO_MED( (int) (image[1] - image[0] + 1) * (int) (image[3] - image[2] + 1), buf, &pmedian); image[0] = (float) MYMAX (c1, MYMIN (oc1 - ideltax, c2)); image[1] = (float) MYMIN (c2, MYMAX (oc2 - ideltax, c1)); image[2] = (float) MYMAX (l1, MYMIN (ol1 - ideltay, l2)); image[3] = (float) MYMIN (l2, MYMAX (ol2 - ideltay, l1)); GETDAT (im, size, npixi, image, ism, buf); MO_MED( (int) (image[1] - image[0] + 1) * (int) (image[3] - image[2] + 1), buf, &median); tdif = tdif + median + dif - pmedian - pdif; if ( deltai[num-1] != MO_INDEFR) deltai[num-1] = deltai[num-1] - tdif; } if (MO_RASTER == "YES") { num = fac - num + 1; fac = fac + fac; } else num = num + nrasters; if (num > nimages) { count = count + 1; num = count; fac = 2 * nrasters; } } while (count <= nrasters); osmmfree((char *) buf); }