C @(#)flip.for 17.1.1.1 (ESO-DMD) 01/25/02 17:40:57 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM FLIP C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program FLIP version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frame, flipping C C.PURPOSE C flip an image in x and/or y dimension C C.ALGORITHM C straight forward, data is changed in place C C.INPUT/OUTPUT C the following keywords are used: C C IN_A/C/1/60 name of input frame C ACTION/C/1/2 action C = 'X ', flip in x C 'Y ', flip in y C 'XY', flip in x and y C 'SH', shift the image C 'SC', scale the image C 'T ', apply thinning algorithm C 'C ', for correlation C C if ACTION = SH, then C OUT_A/C/1/60 name of result fram C INPUTI/I/1/2 xshift,yshift in [1,NPIX-1] C C if ACTION = SC, then C OUT_A/C/1/60 name of result fram C INPUTI/I/1/1 scaling value, C > 1 - enlarge image by replication C < -1 - decrease image size by omission C ACTION(3:3) A or N for averaging or omission C if ACTION = T, then C OUT_A/C/1/60 name of result fram C C if ACTION = C, then C IN_A/C/1/60 template spectrum C IN_B/C/1/1/60 spectrum C OUT_A/C/1/60 correlation table C INPUTI/I/1/1 max shift in pixels C ACTION(2:2) Y or N for normalizing the pixels C of the output frame C C.VERSIONS C C 020108 last modif C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,NAXIS,STAT,KFLAG,FORMA,IOMOD INTEGER*8 PNTR,PNTRB,PNTRC,WPNTR INTEGER NPIX(2),IMNO,NPIXC(2),IMNOC,SCALE(3),NOUT INTEGER NAXISB,NPIXB(2),IMNOB,ISHI INTEGER SIZE,XSHIFT,YSHIFT INTEGER UNI(1),NULO,MADRID(1) C DOUBLE PRECISION START(2),STEP(2),STARTB(2),STEPB(2) DOUBLE PRECISION END(2),DIFSTP,EPS C REAL CUTS(2),MINVAL,MAXVAL,LHCUTS(4) C CHARACTER CUNIT*48,IDENT*72 CHARACTER FRAME*60,FRAMEB*60,FRAMEC*60,ACTION*3 C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA CUTS /0.,0./ DATA IDENT /' '/, CUNIT /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get into MIDAS CALL STSPRO('FLIP') C C get name of input image + action flags CALL STKRDC('IN_A',1,1,60,IAV,FRAME,UNI,NULO,STAT) CALL STKRDC('ACTION',1,1,3,IAV,ACTION,UNI,NULO,STAT) CALL UPCAS(ACTION,ACTION) C C setup action flag IF ((ACTION(1:1).EQ.'X') .OR. (ACTION(1:1).EQ.'Y')) THEN KFLAG = 1 FORMA = D_R4_FORMAT IOMOD = F_IO_MODE ELSE IF (ACTION(1:2).EQ.'SH') THEN KFLAG = 2 FORMA = D_R4_FORMAT IOMOD = F_I_MODE ELSE IF (ACTION(1:2).EQ.'SC') THEN KFLAG = 3 FORMA = D_R4_FORMAT IOMOD = F_I_MODE ELSE IF (ACTION(1:1).EQ.'T') THEN KFLAG = 4 FORMA = D_I4_FORMAT IOMOD = F_I_MODE ELSE IF (ACTION(1:1).EQ.'C') THEN KFLAG = 5 FORMA = D_R4_FORMAT IOMOD = F_I_MODE ELSE CALL STETER(77,'Invalid action ...') ENDIF C C map input frame CALL STIGET(FRAME,FORMA,IOMOD,F_IMA_TYPE, + 2,NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTR,IMNO,STAT) CALL STDRDR(IMNO,'LHCUTS',1,4,IAV,LHCUTS,UNI,NULO,STAT) C C -------------------------- C branch according to action C -------------------------- C IF (KFLAG.EQ.1) THEN !flipping ... CALL STFXMP(NPIX(1),D_R4_FORMAT,WPNTR,STAT) CALL FLIPPI(MADRID(PNTR),MADRID(WPNTR),NPIX,ACTION) C C update START + STEP descriptors accordingly IF (ACTION(1:1).EQ.'Y') THEN END(2) = START(2) + (NPIX(2)-1)*STEP(2) START(2) = END(2) STEP(2) = -STEP(2) ELSE IF (ACTION(2:2).EQ.'Y') THEN END(1) = START(1) + (NPIX(1)-1)*STEP(1) END(2) = START(2) + (NPIX(2)-1)*STEP(2) START(1) = END(1) START(2) = END(2) STEP(1) = -STEP(1) STEP(2) = -STEP(2) ELSE END(1) = START(1) + (NPIX(1)-1)*STEP(1) START(1) = END(1) STEP(1) = -STEP(1) ENDIF CALL STDWRD(IMNO,'START',START,1,2,UNI,STAT) CALL STDWRD(IMNO,'STEP',STEP,1,2,UNI,STAT) IMNOC = IMNO !needed for DSCUPT ... GOTO 9000 ENDIF C C for shifting, thinning, scaling and correlation get output frame C CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,UNI,NULO,STAT) C C IF (KFLAG.EQ.2) THEN !shifting ... CALL STKRDI('INPUTI',1,2,IAV,NPIXC,UNI,NULO,STAT) XSHIFT = NPIXC(1) YSHIFT = NPIXC(2) 500 IF (XSHIFT.GT.NPIX(1)) THEN XSHIFT = XSHIFT - NPIX(1) GOTO 500 ELSE IF (XSHIFT.LT.0) THEN XSHIFT = NPIX(1) + XSHIFT !remember, XSHIFT is < 0 GOTO 500 ENDIF 600 IF (YSHIFT.GT.NPIX(2)) THEN YSHIFT = YSHIFT - NPIX(2) GOTO 600 ELSE IF (YSHIFT.LT.0) THEN YSHIFT = NPIX(2) + YSHIFT !remember, XSHIFT is < 0 GOTO 600 ENDIF IF ((XSHIFT.EQ.0).AND.(YSHIFT.EQ.0)) GOTO 9000 CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTRC,IMNOC,STAT) C C C if shift in both directions, get virtual memory for working buffer IF ((XSHIFT.NE.0).AND.(YSHIFT.NE.0)) THEN SIZE = NPIX(1) * NPIX(2) CALL STFXMP(SIZE,D_R4_FORMAT,WPNTR,STAT) ELSE WPNTR = PNTR ENDIF C CALL SHIFTI(MADRID(PNTR),MADRID(WPNTR),MADRID(PNTRC), + NPIX,XSHIFT,YSHIFT) CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,1,4,UNI,STAT) C ELSE IF (KFLAG.EQ.3) THEN !scaling ... CALL STKRDI('INPUTI',1,2,IAV,SCALE,UNI,NULO,STAT) IF ((SCALE(1).GE.-1).AND.(SCALE(1).LE.1)) SCALE(1) = 1 IF ((SCALE(2).GE.-1).AND.(SCALE(2).LE.1)) SCALE(2) = 1 IF (ACTION(3:3).EQ.'A') THEN SCALE(3) = 1 ELSE SCALE(3) = 0 ENDIF C C calculate new nopix IF (SCALE(1).GT.1) THEN NPIXC(1) = NPIX(1) * SCALE(1) ELSE IF (SCALE(1).LT.1) THEN NPIXC(1) = ( (NPIX(1)-1) / (-SCALE(1)) ) + 1 ELSE NPIXC(1) = NPIX(1) ENDIF IF (SCALE(2).GT.1) THEN NPIXC(2) = NPIX(2) * SCALE(2) ELSE IF (SCALE(2).LT.1) THEN NPIXC(2) = ( (NPIX(2)-1) / (-SCALE(2)) ) + 1 ELSE NPIXC(2) = NPIX(2) ENDIF CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIXC,START,STEP,IDENT, + CUNIT,PNTRC,IMNOC,STAT) CALL SCALI(SCALE,MADRID(PNTR),NPIX,MADRID(PNTRC),NOUT) CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,1,4,UNI,STAT) C ELSE IF (KFLAG.EQ.4) THEN !thinning ... CALL STIPUT(FRAMEC,D_I4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTRC,IMNOC,STAT) C C first copy data SIZE = NPIX(1) * NPIX(2) C CALL COPYI(MADRID(PNTR),MADRID(PNTRC),SIZE) C C get virtual memory for working buffer CALL STFXMP(SIZE,D_I4_FORMAT,WPNTR,STAT) C C and do it CALL THINNI(MADRID(PNTRC),NPIX,MADRID(WPNTR)) C ELSE !correlation ... CALL STKRDC('IN_B',1,1,60,IAV,FRAMEB,UNI,NULO,STAT) !get spectrum CALL STKRDI('INPUTI',1,1,IAV,ISHI,UNI,NULO,STAT) !get max. shift C CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, !map spectrum + 2,NAXISB,NPIXB,STARTB,STEPB,IDENT, + CUNIT,PNTRB,IMNOB,STAT) EPS = 0.001 * ABS(STEPB(1)) ! 0.1 % of x-step DIFSTP = ABS(STEP(1)) - ABS(STEPB(1)) IF (DIFSTP.GT.EPS) + CALL STETER(1, + 'x-step of template and spec frame not equal...') C C how many points will we have in overlap region? SIZE = NPIX(1) - ISHI - ISHI IF (SIZE.LT.3) THEN CALL STTPUT(' Too little overlap ',STAT) CALL STETER + (2,'Shift too large or template frame too small...') ENDIF C NAXIS = 1 NPIXC(1) = 2*ISHI + 1 START(1) = -(ISHI * STEPB(1)) CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIXC,START,STEPB,IDENT, + CUNIT,PNTRC,IMNOC,STAT) C C and do it IF ((ACTION(2:2).EQ.'Y') .OR. (ACTION(2:2).EQ.'y')) THEN CALL CTCORR(NPIX(1),MADRID(PNTR),NPIXB(1),MADRID(PNTRB), + MADRID(PNTRC),ISHI,MINVAL,MAXVAL) ELSE CALL XSCORR(NPIX(1),MADRID(PNTR),NPIXB(1),MADRID(PNTRB), + MADRID(PNTRC),ISHI,MINVAL,MAXVAL) ENDIF LHCUTS(1) = MINVAL LHCUTS(2) = MAXVAL LHCUTS(3) = MINVAL LHCUTS(4) = MAXVAL CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,1,4,UNI,STAT) ENDIF C C that's it folks... 9000 CALL DSCUPT(IMNO,IMNOC,' ',STAT) !copy all other descriptors CALL STSEPI END SUBROUTINE FLIPPI(A,R,NPIX,ACTION) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subrotine FLIPPI version 1.0 870605 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frame, flipping C C.PURPOSE C flip an image or mask in x and/or y dimension C C.ALGORITHM C straight forward, data is changed in place C C.INPUT/OUTPUT C call as FLIPPI(A,R,NPIX,ACTION) C C A: R*4 array input image buffer C R: R*4 array buffer for single line C NPIX: I*4 array dimensions of A C ACTION: char.exp flip action C = 'X ', flip in x C 'Y ', flip in y C 'XY', flip in x and y C C.VERSIONS C 1.00 from old FLIPPY C C-------------------------------------------------- C C IMPLICIT NONE C INTEGER NPIX(2) INTEGER FROM,GAP,GAPP,INDX,NTMP INTEGER NX,NY,TO,XDIM,XMID,YDIM,YMID C REAL A(*),R(*) REAL SAVE C CHARACTER*(*) ACTION C C init XDIM = NPIX(1) YDIM = NPIX(2) XMID = XDIM/2 YMID = YDIM/2 C C branch according to action IF (ACTION(1:1).EQ.'Y') GOTO 1000 IF (ACTION(1:2).EQ.'XY') GOTO 2000 C C ACTION = 'X ', flip only in x C FROM = 0 !offset to first element in a row DO 500 NY=1,YDIM TO = FROM + XDIM + 1 !point to last + 1 C DO 300 NX=1,XMID !loop till middle of row SAVE = A(NX+FROM) !flip... A(NX+FROM) = A(TO-NX) A(TO-NX) = SAVE 300 CONTINUE FROM = FROM + XDIM 500 CONTINUE RETURN C C ACTION = 'Y ', flip only in y C 1000 FROM = 0 !offset to first element in a row DO 1500 NY=1,YMID !loop till middle of columns TO = YDIM + 1 - NY GAP = (TO-NY) * XDIM NTMP = FROM + GAP C DO 1300 NX=1,XDIM !loop through row INDX = NX + FROM R(NX) = A(INDX) !save row elements in R 1300 CONTINUE DO 1330 NX=1,XDIM INDX = NX + FROM A(INDX) = A(INDX+GAP) !flip complete row... 1330 CONTINUE DO 1360 NX=1,XDIM INDX = NX + NTMP A(INDX) = R(NX) 1360 CONTINUE FROM = FROM + XDIM 1500 CONTINUE RETURN C C ACTION = 'XY', flip in x and y C 2000 FROM = 0 DO 2500 NY=1,YMID !loop till middle of columns TO = YDIM + 1 - NY GAP = (TO-NY) * XDIM GAPP = GAP + XDIM + 1 NTMP = FROM + GAPP C DO 2300 NX=1,XDIM !loop through row INDX = NX + FROM R(NX) = A(INDX) !save row elements in R 2300 CONTINUE DO 2330 NX=1,XDIM INDX = NX + FROM A(INDX) = A(NTMP-NX) !flip complete row... 2330 CONTINUE DO 2360 NX=1,XDIM INDX = NTMP - NX A(INDX) = R(NX) 2360 CONTINUE FROM = FROM + XDIM 2500 CONTINUE C IF (2*YMID.NE.YDIM) THEN !take care of middle row FROM = YMID * XDIM TO = FROM + XDIM + 1 DO 2800 NX=1,XMID INDX = NX + FROM SAVE = A(INDX) A(INDX) = A(TO-NX) A(TO-NX) = SAVE 2800 CONTINUE ENDIF C RETURN END SUBROUTINE SHIFTI(A,B,C,NPIX,XSHIFT,YSHIFT) C IMPLICIT NONE C INTEGER NPIX(2),XSHIFT,YSHIFT INTEGER NY,NX,NOFF,OFF1,OFF2,K1 C REAL A(*),B(*),C(*) C IF (XSHIFT.EQ.0) GOTO 3000 C K1 = NPIX(1) - XSHIFT OFF1 = 0 C C check, if we have to shift also y later on IF (YSHIFT.EQ.0) THEN DO 1800 NY=1,NPIX(2) !shift in x for each line DO 1000 NX=1,K1 NOFF = NX + OFF1 C(NOFF+XSHIFT) = A(NOFF) 1000 CONTINUE DO 1500 NX=1,XSHIFT NOFF = NX + OFF1 C(NOFF) = A(NOFF+K1) 1500 CONTINUE OFF1 = OFF1 + NPIX(1) 1800 CONTINUE RETURN C ELSE DO 2800 NY=1,NPIX(2) !shift in x for each line DO 2000 NX=1,K1 NOFF = NX + OFF1 B(NOFF+XSHIFT) = A(NOFF) 2000 CONTINUE DO 2500 NX=1,XSHIFT NOFF = NX + OFF1 B(NOFF) = A(NOFF+K1) 2500 CONTINUE OFF1 = OFF1 + NPIX(1) 2800 CONTINUE ENDIF C 3000 K1 = NPIX(2) - YSHIFT C C and now shift complete lines (initially B = A ) OFF1 = 0 OFF2 = YSHIFT * NPIX(1) DO 4400 NY=1,K1 DO 4000 NX=1,NPIX(1) C(OFF2+NX) = B(OFF1+NX) 4000 CONTINUE OFF1 = OFF1 + NPIX(1) !move to next line OFF2 = OFF2 + NPIX(1) 4400 CONTINUE C OFF1 = 0 OFF2 = K1 * NPIX(1) DO 5500 NY=1,YSHIFT DO 5000 NX=1,NPIX(1) C(OFF1+NX) = B(OFF2+NX) 5000 CONTINUE OFF1 = OFF1 + NPIX(1) !move to next line OFF2 = OFF2 + NPIX(1) 5500 CONTINUE C RETURN END SUBROUTINE SCALI(SCALE,IN,NPIX,OUT,NOUT) C IMPLICIT NONE C INTEGER SCALE(3),NPIX(2),NOUT INTEGER OUTNDX,INOFF,KSCAL,INDX INTEGER NLINE,NY,NX,K,KK,NPX,M,JOUT C REAL IN(*),OUT(*) REAL VAL,AVFACT C OUTNDX = 1 INOFF = 0 C C split according to y-SCALE IF (SCALE(2).LT.1) GOTO 10000 C C we keep (or increase) the no. of lines C IF (SCALE(1).GE.1) THEN !increase in x C KSCAL = SCALE(1) - 1 DO 2000, NLINE=1,NPIX(2) DO 1800, NY=1,SCALE(2) !replicate lines DO 1600, NX=1,NPIX(1) INDX = INOFF + NX !point to input pixel DO 1000, K=OUTNDX,OUTNDX+KSCAL !replicate pixels OUT(K) = IN(INDX) 1000 CONTINUE OUTNDX = OUTNDX + SCALE(1) !increment out pointer 1600 CONTINUE 1800 CONTINUE INOFF = INOFF + NPIX(1) !point to next line 2000 CONTINUE C ELSE !decrease in x SCALE(1) = -SCALE(1) IF (SCALE(3).EQ.0) THEN !we omit pixels DO 2300, NLINE=1,NPIX(2) DO 2200, NY=1,SCALE(2) !replicate lines DO 2100, NX=1,NPIX(1),SCALE(1) OUT(OUTNDX) = IN(INOFF+NX) OUTNDX = OUTNDX + 1 2100 CONTINUE 2200 CONTINUE INOFF = INOFF + NPIX(1) !point to next line 2300 CONTINUE ELSE AVFACT = 1.0/SCALE(1) DO 3400,NLINE=1,NPIX(2) DO 3300,NY=1,SCALE(2) DO 3200,NX=1,NPIX(1),SCALE(1) VAL = 0. M = INOFF+NX-1 DO 3100,K=1,SCALE(1) VAL = VAL + IN(M+K) 3100 CONTINUE OUT(OUTNDX) = VAL*AVFACT OUTNDX = OUTNDX + 1 3200 CONTINUE 3300 CONTINUE C INOFF = INOFF + NPIX(1) !point to next line 3400 CONTINUE ENDIF ENDIF GOTO 20000 C C we decrease the no. of lines C 10000 SCALE(2) = -SCALE(2) KSCAL = SCALE(1) - 1 C IF (SCALE(1).GE.1) THEN !increase in x KK = NPIX(1)*SCALE(1) - 1 C IF (SCALE(3).EQ.0) THEN DO 12600, NLINE=1,NPIX(2),SCALE(2) DO 12500, NX=1,NPIX(1) INDX = INOFF + NX !point to input pixel DO 12400, K=OUTNDX,OUTNDX+KSCAL !replicate pixels OUT(K) = IN(INDX) 12400 CONTINUE OUTNDX = OUTNDX + SCALE(1) !increment out pointer 12500 CONTINUE INOFF = INOFF + (NPIX(1)*SCALE(2)) !point to next line 12600 CONTINUE C ELSE AVFACT = 1.0/SCALE(2) JOUT = OUTNDX DO 13600, NLINE=1,NPIX(2),SCALE(2) DO 13100, K=OUTNDX,OUTNDX+KK !init sums OUT(K) = 0.0 13100 CONTINUE C DO 13500, NY=1,SCALE(2) DO 13400, NX=1,NPIX(1) INDX = INOFF + NX !point to input pixel DO 13300, K=OUTNDX,OUTNDX+KSCAL !sum up OUT(K) = OUT(K) + IN(INDX) 13300 CONTINUE OUTNDX = OUTNDX + SCALE(1) !increment out pointer 13400 CONTINUE OUTNDX = JOUT INOFF = INOFF + NPIX(1) !point to next line 13500 CONTINUE C DO 13550, K=OUTNDX,OUTNDX+KK OUT(K) = OUT(K) * AVFACT 13550 CONTINUE JOUT = JOUT + (NPIX(1)*SCALE(1)) OUTNDX = JOUT 13600 CONTINUE ENDIF C ELSE !decrease in x SCALE(1) = -SCALE(1) IF (SCALE(3).EQ.0) THEN DO 14400, NLINE=1,NPIX(2),SCALE(2) DO 14200, NX=1,NPIX(1),SCALE(1) OUT(OUTNDX) = IN(INOFF+NX) OUTNDX = OUTNDX + 1 14200 CONTINUE INOFF = INOFF + (NPIX(1)*SCALE(2)) 14400 CONTINUE C ELSE JOUT = OUTNDX AVFACT = 1.0/(SCALE(1)*SCALE(2)) DO 16000, NLINE=1,NPIX(2),SCALE(2) DO 15100, NPX=1,NPIX(1),SCALE(1) OUT(OUTNDX) = 0.0 OUTNDX = OUTNDX + 1 15100 CONTINUE C OUTNDX = JOUT DO 15800,NY=1,SCALE(2) DO 15700,NPX=1,NPIX(1),SCALE(1) INDX = INOFF + NPX - 1 DO 15600, NX=1,SCALE(1) OUT(OUTNDX) = OUT(OUTNDX) + IN(INDX+NX) 15600 CONTINUE OUTNDX = OUTNDX + 1 15700 CONTINUE OUTNDX = JOUT INOFF = INOFF + NPIX(1) 15800 CONTINUE C DO 15900, NPX=1,NPIX(1),SCALE(1) OUT(OUTNDX) = OUT(OUTNDX)*AVFACT OUTNDX = OUTNDX + 1 15900 CONTINUE JOUT = OUTNDX !follow with base index 16000 CONTINUE ENDIF ENDIF C 20000 NOUT = OUTNDX - 1 !return size of output frame C RETURN END SUBROUTINE THINNI(M,NPIX,W) C IMPLICIT NONE C INTEGER M(*),NPIX(2),W(*) INTEGER REMAIN,INC,IOFF,IOFFX INTEGER MX,MY,JN,NX,NY INTEGER JOFF(4),N1,N2,N3,N4,N5,N6,N7,N8,KOFF INTEGER IOFA C NX = NPIX(1) NY = NPIX(2) JOFF(1) = 1 JOFF(2) = NX JOFF(3) = -1 JOFF(4) = -NX REMAIN = 1 C C main loop 200 IF (REMAIN.NE.1) RETURN !We're done. C REMAIN = 0 C DO 5000, JN=1,4 !look at all direct neighbours... KOFF = JOFF(JN) INC = 0 IOFFX = NX DO 4000, MY=2,NY-1 !omit first + last line... DO 3000, MX=2,NX-1 !omit first + last pixel... IOFF = IOFFX + MX IF ( (M(IOFF).EQ.1) .AND. + (M(IOFF+KOFF).EQ.0) ) THEN !get all neighbours... IOFA = IOFF + NX N1 = M(IOFF+1) N2 = M(IOFA+1) N3 = M(IOFA) N4 = M(IOFA-1) IOFA = IOFF - NX N5 = M(IOFF-1) N6 = M(IOFA-1) N7 = M(IOFA) N8 = M(IOFA+1) C C now check conditions IF (N1.EQ.0) THEN C IF (N5.EQ.0) THEN ! 0AAA0BBB IF ((N2+N3+N4.GT.0).AND.(N6+N7+N8.GT.0)) + GOTO 1000 ELSE IF ((N7.EQ.0).AND.(N8.EQ.2)) THEN ! 0AAAAA02 IF (N2+N3+N4+N5+N6.GT.0) GOTO 1000 ELSE IF ((N3.EQ.0).AND.(N2.EQ.2)) THEN ! 020AAAAA IF (N4+N5+N6+N7+N8.GT.0) GOTO 1000 ENDIF ENDIF C IF (N3.EQ.0) THEN C IF (N7.EQ.0) THEN ! BB0AAA0B IF ((N1+N2+N8.GT.0).AND.(N4+N5+N6.GT.0)) + GOTO 1000 ELSE IF ((N5.EQ.0).AND.(N4.EQ.2)) THEN ! AA020AAA IF (N1+N2+N6+N7+N8.GT.0) GOTO 1000 ENDIF ENDIF C IF ((N5+N7.EQ.0).AND.(N6.EQ.2)) THEN ! AAAA020A IF (N1+N2+N3+N4+N8.GT.0) GOTO 1000 ENDIF C M(IOFF) = 3 !no pattern match... REMAIN = 1 INC = INC + 1 W(INC) = IOFF !save index... GOTO 3000 C 1000 M(IOFF) = 2 !yes, pattern match... ENDIF C 3000 CONTINUE !end of MX loop... C IOFFX = IOFFX + NX 4000 CONTINUE !end of MY loop... C C change all 3's to 0's IF (INC.GT.0) THEN DO 4400, MX=1,INC M(W(MX)) = 0 4400 CONTINUE ENDIF C 5000 CONTINUE !end of JN loop... C C end of main loop GOTO 200 C END SUBROUTINE XSCORR(NTEMP,TEMP,NSPEC,SPEC,RES,NSH,MINVAL,MAXVAL) C C This program computes the cross-correlation between a template C and a spectrum over a bandwith of 2*NSH pixels. Identity of C start and step are assumed. We will have 2*NSH + 1 data. C C IMPLICIT NONE C INTEGER NTEMP,NSPEC,NSH INTEGER MPS,MPE,INC,I,J C REAL TEMP(*),SPEC(*),RES(*) REAL MINVAL,MAXVAL C DOUBLE PRECISION SUM C MPS = 1 + NSH IF (NSPEC .GT. NTEMP) THEN !watch out for that MPE = NTEMP - NSH ELSE MPE = NSPEC - NSH ENDIF IF (MPE.LT.MPS) + CALL STETER(2,'Shift too large or frames too small...') INC = 1 C DO 20, I=-NSH,NSH SUM = 0.D0 DO 10, J=MPS,MPE SUM = SUM + ( SPEC(J)*TEMP(J+I) ) 10 CONTINUE C RES(INC) = SNGL(SUM) INC = INC + 1 20 CONTINUE C C get min, max MINVAL = RES(1) MAXVAL = MINVAL DO 100, I=2,INC-1 IF (RES(I).GT.MAXVAL) THEN MAXVAL = RES(I) ELSE IF (RES(I).LT.MINVAL) THEN MINVAL = RES(I) ENDIF 100 CONTINUE C RETURN END SUBROUTINE CTCORR(NTEMP,TEMP,NSPEC,SPEC,RES,NSH,MINVAL,MAXVAL) C C This program computes the cross-correlation between a template C and a spectrum over a bandwith of 2*NSH pixels. Identity of C start and step are assumed. We will have 2*NSH + 1 data. C M. Rosa C CT Modified to compute true correlation coefficient with optimal CT overlap range. CT For two arrays x and y the correlation coefficient is CT r = Covar(x,y)/Sqrt[Var(x)*Var(y)], CT which gives CT r = [n*sum(x*y) - sum(x)*sum(y)]/ CT sqrt{[n*sum(x**2)-sum(x)**2]*[n*sum(y**2)-sum(y)**2]} CT Since r is normalized one can use the total overlap of two CT shifted arrays, varying with the shift variable. CT 29.11.2001 H.-C. Thomas CT Max-Planck Institute for Astrophysics C C IMPLICIT NONE C INTEGER NTEMP,NSPEC,NSH INTEGER MPS,MPE,INC,I,J C REAL TEMP(*),SPEC(*),RES(*) REAL MINVAL,MAXVAL C DOUBLE PRECISION SUMX,SUMY,SUMXX,SUMYY,SUMXY,NPOI,T1,T2,T3 C CT compute the smallest overlap, exit if less than 1 CT MPS = 1 IF (NSPEC .GT. NTEMP) THEN !watch out for that MPE = NTEMP - NSH ELSE MPE = NSPEC - NSH ENDIF IF (MPE.LT.MPS) + CALL STETER(2,'Shift too large or frames too small...') INC = 1 C CT compute the necessary sums over a range MPS to MPE varying with shift I CT DO 20, I=-NSH,NSH SUMX = 0.D0 SUMXX = 0.D0 SUMY = 0.D0 SUMYY = 0.D0 SUMXY = 0.D0 MPS = 1 IF (I .GT. 0) MPS = 1+I MPE = NSPEC+I IF (MPE .GT. NTEMP) MPE = NTEMP DO 10, J=MPS,MPE SUMX = SUMX + TEMP(J) SUMXX = SUMXX + TEMP(J)**2 SUMY = SUMY + SPEC(J-I) SUMYY = SUMYY + SPEC(J-I)**2 SUMXY = SUMXY + ( TEMP(J)*SPEC(J-I) ) 10 CONTINUE C CT compute the correlation coefficient CT NPOI = DFLOAT(1+MPE-MPS) T1 = NPOI*SUMXY-SUMX*SUMY T2 = NPOI*SUMXX-SUMX**2 T3 = NPOI*SUMYY-SUMY**2 RES(INC) = SNGL(T1/DSQRT(T2*T3)) INC = INC + 1 20 CONTINUE C C get min, max MINVAL = RES(1) MAXVAL = MINVAL DO 100, I=2,INC-1 IF (RES(I).GT.MAXVAL) THEN MAXVAL = RES(I) ELSE IF (RES(I).LT.MINVAL) THEN MINVAL = RES(I) ENDIF 100 CONTINUE C RETURN END