C @(#)genxx1.for 13.1.1.4 (ESO-DMD) 02/10/99 15:43:45 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 Massachusetts 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 GENXX1 C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program GENXX1 version 1.00 890704 C K. Banse ESO - Garching C 1.10 891124 C C.KEYWORDS C image processing functions C C.PURPOSE C merger of several functions, see description there C C.ALGORITHM C see local descriptions C C.INPUT/OUTPUT C the following keys are used: C C ACTION/C/1/2 action flag C = ED, for edge detection C = EX, for general extract routines C = EZ, for reference frame extraction C = IN, for inserting an image C = MP, for mapping images C = MA, for matrix operations C C.VERSION C 1.00 merge former individual main programmes C for later versions see SCCS C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,STAT INTEGER UNI(1),NULO C CHARACTER ACTION*2 C C set up MIDAS environment + enable automatic error abort CALL STSPRO('GENXX1') C C get option CALL STKRDC('ACTION',1,1,2,IAV,ACTION,UNI,NULO,STAT) C IF (ACTION.EQ.'ED') THEN CALL SUBEDG ELSE IF (ACTION.EQ.'EX') THEN CALL SUBEXT ELSE IF (ACTION.EQ.'EZ') THEN CALL SUBEZT ELSE IF (ACTION.EQ.'IN') THEN CALL SUBINS ELSE IF (ACTION.EQ.'MP') THEN CALL SUBMAP ELSE IF (ACTION.EQ.'MA') THEN CALL SUBMAT ENDIF C C get away CALL STSEPI END SUBROUTINE SUBEDG C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBEDG version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C edge detection C C.PURPOSE C detect edges in a given data frame C C.ALGORITHM C get the names of input + output frames from IN_A + OUT_A C get the threshold for algorithm from INPUTR C and calculate edge pixels via subroutine FNDEDGE (see description C of algorithm there) C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/80 input data frame C P4/C/1/1 edge detection method C = B, use band - method C = S, use Sobel operator C OUT_A/C/1/80 output data frame C DEFAULT/C/1/1 default flag, if = Y, use half of dynamic range of input frame C as threshold, else read it from INPUTR C INPUTR/R/1/1 threshold for edge detection C C.VERSION C 1.00 from version 2.00 as of 841120 C move to FORTRAN 77 + new ST interfaces C C-------------------------------------------------- C IMPLICIT NONE C INTEGER NAXISA,IAV,STAT INTEGER*8 PNTRA,PNTRC INTEGER IMNOA,IMNOC INTEGER NPIXA(3) INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*80 FRAMEA,FRAMEC CHARACTER CUNITA*64,IDENTA*72 CHARACTER DEFAUL*1,ACTION*1 C REAL THRESH,LHCUTS(2) REAL OLDCUT(2) C DOUBLE PRECISION STEPA(3),STARTA(3) C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA LHCUTS /0.,1./ DATA IDENTA /' '/, CUNITA /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get input frame + map it CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 3,NAXISA,NPIXA,STARTA,STEPA, + IDENTA,CUNITA,PNTRA,IMNOA,STAT) C C get output frame + map it (=edge file) + copy standard descriptors CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXA,STARTA,STEPA, + IDENTA,CUNITA,PNTRC,IMNOC,STAT) C C get threshold + method CALL STKRDC('DEFAULT',1,1,1,IAV,DEFAUL,UNI,NULO,STAT) IF ( (DEFAUL.EQ.'Y') .OR. (DEFAUL.EQ.'y') ) THEN CALL STDRDR(IMNOA,'LHCUTS',3,2,IAV,OLDCUT,UNI,NULO,STAT) THRESH = (OLDCUT(2) - OLDCUT(1)) * .5 ELSE CALL STKRDR('INPUTR',1,1,IAV,THRESH,UNI,NULO,STAT) ENDIF CALL STKRDC('P4',1,1,1,IAV,ACTION,UNI,NULO,STAT) CALL UPCAS(ACTION,ACTION) C C now do the edge detection CALL FNDEDG(ACTION,MADRID(PNTRA),NPIXA(1),NPIXA(2), + THRESH,MADRID(PNTRC)) CALL STDWRR(IMNOC,'LHCUTS',LHCUTS,3,2,UNI,STAT) C RETURN END SUBROUTINE SUBEXT C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBEXT version 1.00 881028 C K. Banse ESO - Garching C 981016 C C.KEYWORDS C bulk data frames C C.PURPOSE C extract a subframe from a given frame C C.ALGORITHM C get output frame from OUT_A C get relevant coordinates as character data from INPUTC or via cursor C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/80 input frame C OUT_A/C/1/80 result frame C ACTION/C/3/2 action flag C action flag = IM: extract subframe from coord interval C = IC: extract subframe around center coord C = LI: extract a line in any direction from a higher C dimensioned frame C if action flag = IM C INPUTC/C/1/80 character string of coordinates C in the form: [x1,y1,z1:x2,y2,z2] C x,y,z according to the MIDAS standard for C specifying coordinates C or if centered extract C INPUTC/C/1/80 [xc,yc,zc] C INQXI/I/I/1/3 xl,yl,zl (0,...,NPIX) C INQXI/I/I/4/3 xr,yr,zr (0,...,NPIX) C if action flag = LI C INPUTC/C/1/80 character string of coordinates C in the form: [x1,y1,z1:x2,y2,z2] C INPUTD/D/1/1 stepsize for line C C.VERSIONS C see SCCS C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,IAX,N,NINDX,NN,REFSET INTEGER NAXISA,NAXISC,STAT INTEGER*8 PNTRA,PNTRC,PNTRX,PNTRY INTEGER NPIXA(3),NPIXC(3) INTEGER RETBUF(3),SUBLO(3),SUBHI(3) INTEGER MADRID(1),UNI(1),NULO INTEGER IMNOA,IMNOC INTEGER NPXY,K2,RMAIND,CHUCHO,CHUCHI INTEGER NOLIN,TRUSIZ,CHUNKI,CHUNKO,K1,KOFF,MOFF,LOOPLM INTEGER PXOFF(6),PCENT(3) C CHARACTER*80 FRAMEA,FRAMEC,EQUAL CHARACTER CUNITA*64,ACTION*2 CHARACTER INPUTC*80 CHARACTER ERROR1*40,ERROR2*40,ERROR3*40 CHARACTER*72 IDENTA C DOUBLE PRECISION STARTA(3),STARTC(3),STEPA(3),STEPC(3),STEPLN DOUBLE PRECISION DIN(4),DOUT(4) C REAL CUTS(4),RM,DIST,WW(2) REAL XA,XB,YA,YB C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C INCLUDE 'MID_INCLUDE:IDIDEV.INC' INCLUDE 'MID_INCLUDE:IDIMEM.INC' C COMMON /VMR/ MADRID C DATA ERROR1 /'invalid coordinate string ... '/ DATA ERROR2 /'problems with dynamic range... '/ DATA ERROR3 /'max. 3-dim images supported... '/ DATA STEPA /3*1.D0/, STARTA /3*0.D0/ DATA STEPC /3*1.D0/, STARTC /3*0.D0/ DATA NPIXA /3*1/, NPIXC /3*1/ DATA SUBLO /1,1,1/, SUBHI /1,1,1/ DATA CUTS /0.,0.,0.,0./ DATA IDENTA /' '/ DATA CUNITA /' '/ DATA PXOFF /6*0/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get name of input frame, result frame, descr_copy_flag + action flag CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) CALL CLNFRA(FRAMEA,FRAMEA,0) CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) CALL STKRDC('ACTION',1,3,2,IAV,ACTION,UNI,NULO,STAT) C C get string with coordinates REFSET = 0 !set to 1, if we update REFPIX descr. CALL STKRDC('INPUTC',1,1,80,IAV,INPUTC,UNI,NULO,STAT) C C check if we shrink same frame CALL GENEQF(FRAMEA,FRAMEC,STAT) IF (STAT.EQ.1) THEN !if same name, use temporary one EQUAL(1:) = FRAMEC(1:) FRAMEC(1:) = 'middummeq.bdf ' ELSE EQUAL(1:) = ' ' ENDIF C C and open input frame or open + already map it IF (ACTION(1:1).EQ.'I') THEN CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT) CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXISA,UNI,NULO,STAT) IF (NAXISA.GT.3) CALL STETER(44,ERROR3) CALL STDRDI(IMNOA,'NPIX',1,NAXISA,IAV,NPIXA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'START',1,NAXISA,IAV,STARTA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'STEP',1,NAXISA,IAV,STEPA,UNI,NULO,STAT) CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENTA,UNI,NULO,STAT) N = (NAXISA+1) * 16 CALL STDRDC(IMNOA,'CUNIT',1,1,N,IAV,CUNITA,UNI,NULO,STAT) C C check if centered extract IF (ACTION(2:2).EQ.'C') THEN !Yes, centered extraction CALL STKRDI('INQXI',1,6,IAV,PXOFF,UNI,NULO,STAT) CALL EXTCO1(IMNOA,INPUTC,NAXISA,NAXISC,PCENT,STAT) IF (STAT.NE.0) CALL STETER(1,ERROR1) SUBLO(1) = PCENT(1) - PXOFF(1) SUBHI(1) = PCENT(1) + PXOFF(4) IF (NAXISC.GT.1) THEN SUBLO(2) = PCENT(2) - PXOFF(2) SUBHI(2) = PCENT(2) + PXOFF(5) IF (NAXISC.GT.2) THEN SUBLO(3) = PCENT(3) - PXOFF(3) SUBHI(3) = PCENT(3) + PXOFF(6) ENDIF ENDIF DO 440, N=1,NAXISC IF ( (SUBLO(N).LT.1) .OR. + (SUBLO(N).GT.SUBHI(N)) .OR. + (SUBHI(N).GT.NPIXA(N)) ) + CALL STETER(4,'invalid pixel offsets ... ') 440 CONTINUE GOTO 480 ENDIF ELSE CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,3, + NAXISA,NPIXA,STARTA,STEPA, + IDENTA,CUNITA,PNTRA,IMNOA,STAT) IF (NAXISA.GT.3) CALL STETER(44,ERROR3) ENDIF CALL FPXWCO(0,IMNOA,DIN,DOUT,STAT) !init wcs conversion C C copy subframe or extract a line CALL EXTCOO(IMNOA,INPUTC,NAXISA,NAXISC,SUBLO,SUBHI,STAT) IF (STAT.NE.0) CALL STETER(1,ERROR1) C C get start + end pixels of subframe within father frame 480 DO 500 N=1,NAXISC DIN(N) = SUBLO(N) !in fp's NPIXC(N) = SUBHI(N) - SUBLO(N) + 1 500 CONTINUE CALL FPXWCO(1,0,DIN,DOUT,STAT) !fp -> wc DO 510 N=1,NAXISC STARTC(N) = STARTA(N) + (SUBLO(N)-1)*STEPA(N) 510 CONTINUE CALL STDFND(IMNOA,'REFPIX',INPUTC,N,NN,STAT) IF (INPUTC(1:1).EQ.'D') THEN CALL STDRDD(IMNOA,'REFPIX',1,NAXISA,IAV,DOUT,UNI,NULO,STAT) DO 8080, N=1,NAXISA DOUT(N) = DOUT(N) + 1.0D0 - DIN(N) 8080 CONTINUE REFSET = 1 ENDIF CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT) C C branch according to line or subframe extraction... IF (ACTION(1:1).EQ.'L') GOTO 5000 C C ------------------------------------------------------- C C here to extract a subframe C C ------------------------------------------------------- C TRUSIZ = 1 !get size of input frame DO 600 N=1,NAXISA TRUSIZ = TRUSIZ * NPIXA(N) 600 CONTINUE C C get buffersize of at least one line but no more than total size CALL HACKUP(NPIXA,D_R4_FORMAT,RETBUF) NOLIN = RETBUF(1) !buffer = NOLIN * NPIXA(1) IF (NOLIN.GT.NPIXC(2)) NOLIN = NPIXC(2) CHUNKI = NOLIN * NPIXA(1) CHUNKO = NOLIN * NPIXC(1) LOOPLM = NPIXC(2) / NOLIN !get no. of passes through loop RMAIND = NPIXC(2) - (LOOPLM*NOLIN) !and remainder C CALL STFXMP(CHUNKI,D_R4_FORMAT,PNTRX,STAT) CALL STFXMP(CHUNKO,D_R4_FORMAT,PNTRY,STAT) C TRUSIZ = 1 DO 800 N=1,NAXISC TRUSIZ = TRUSIZ * NPIXC(N) 800 CONTINUE CALL STFCRE(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + TRUSIZ,IMNOC,STAT) C C now we loop through input frame NPXY = NPIXA(1) * NPIXA(2) K2 = (SUBLO(3)-1) * NPXY !skip unused planes K1 = (SUBLO(2)-1) * NPIXA(1) + 1 !point to 1. used y-lines MOFF = 1 !pointer for output frame C CUTS(3) = 999999. CUTS(4) = -999999. C DO 4000 IAX=1,NPIXC(3) !loop over all planes KOFF = K2 + K1 IF (LOOPLM.LE.0) GOTO 3000 C DO 2000 N=1,LOOPLM CALL STFGET(IMNOA,KOFF,CHUNKI,TRUSIZ,MADRID(PNTRX),STAT) !data in CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),SUBLO(1),1, + NPIXA(1),NOLIN,NPIXA(1),NPIXC(1),CUTS(3)) !rearrange data CALL STFPUT(IMNOC,MOFF,CHUNKO,MADRID(PNTRY),STAT) !data out KOFF = KOFF + CHUNKI !follow with indices... MOFF = MOFF + CHUNKO 2000 CONTINUE C 3000 IF (RMAIND.GT.0) THEN CHUCHI = RMAIND * NPIXA(1) CHUCHO = RMAIND * NPIXC(1) CALL STFGET(IMNOA,KOFF,CHUCHI,TRUSIZ,MADRID(PNTRX),STAT) CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),SUBLO(1),1, + NPIXA(1),RMAIND,NPIXA(1),NPIXC(1),CUTS(3)) CALL STFPUT(IMNOC,MOFF,CHUCHO,MADRID(PNTRY),STAT) MOFF = MOFF + CHUCHO ENDIF C K2 = K2 + NPXY !move to next plane 4000 CONTINUE C C compress NAXISC + write descriptors 4200 IF (NAXISC.GT.1) THEN IF (NPIXC(1).EQ.1) THEN NPIXC(1) = NPIXC(2) NPIXC(2) = NPIXC(3) STARTC(1) = STARTC(2) STARTC(2) = STARTC(3) STEPA(1) = STEPA(2) STEPA(2) = STEPA(3) CUNITA(17:) = CUNITA(33:)//' ' NAXISC = NAXISC - 1 GOTO 4200 ENDIF IF (NPIXC(2).EQ.1) THEN NPIXC(2) = NPIXC(3) STARTC(2) = STARTC(3) STEPA(2) = STEPA(3) CUNITA(33:) = CUNITA(49:)//' ' NAXISC = NAXISC - 1 GOTO 4200 ENDIF IF (NPIXC(NAXISC).EQ.1) NAXISC = NAXISC - 1 ENDIF C CALL STDWRI(IMNOC,'NAXIS',NAXISC,1,1,UNI,STAT) CALL STDWRI(IMNOC,'NPIX',NPIXC,1,NAXISC,UNI,STAT) !still use NAXISC CALL STDWRD(IMNOC,'START',STARTC,1,NAXISC,UNI,STAT) CALL STDWRD(IMNOC,'STEP',STEPA,1,NAXISC,UNI,STAT) CALL STDWRC(IMNOC,'IDENT',1,IDENTA,1,72,UNI,STAT) N = (NAXISC+1) * 16 CALL STDWRC(IMNOC,'CUNIT',1,CUNITA,1,N,UNI,STAT) C GOTO 8000 C C ------------------------------------------------------- C C here to extract a line C C ------------------------------------------------------- C 5000 IF ( (NAXISA.GT.2) .AND. !test that input frame is a 2-d frame... + (NPIXA(3).GT.1) ) THEN NN = INDEX(FRAMEA,' ') - 1 IF (NN.LE.0) NN = LEN(FRAMEA) INPUTC(1:) = 'Input frame '//FRAMEA(1:NN)// + ' has to be 2-dim frame... ' CALL STETER(3,INPUTC) ENDIF C C get stepsize for line CALL STKRDD('INPUTD',1,1,IAV,STEPLN,UNI,NULO,STAT) XA = SUBLO(1) YA = SUBLO(2) XB = SUBHI(1) YB = SUBHI(2) C C determine main axis (depends upon orientation of line in real space) IF (ABS(XA-XB).LT.10.E-5) THEN IAX = 2 !parallel to y-axis... ELSE RM = (YB-YA) / (XB-XA) IF (ABS(RM).GT.1.0001) THEN IAX = 2 !steep ramp, so we use y.. ELSE IAX = 1 !low ramp, so we use x... ENDIF ENDIF IF (ABS(STEPLN).LE.10.E-20) THEN !use orig_step if = 0.0 STEPLN = STEPA(IAX) RM = 1. !RM holds scaled stepsieze ELSE RM = STEPLN/STEPA(IAX) ENDIF C C get no. of points and map index arrays as dummy frame WW(1) = XB - XA WW(2) = YB - YA DIST = SQRT( (WW(1)*WW(1)) + (WW(2)*WW(2)) ) WW(1) = DIST/RM IF (WW(1).LT.0.) WW(1) = -WW(1) NN = IFIX( WW(1) ) WW(1) = 0.01 * WW(1) IF (WW(1).GT.3.0) THEN NN = NN + IFIX( WW(1) ) !add 1 % for safety check ELSE NN = NN + 3 !add 3 for safety check ENDIF C CALL STFXMP(NN,D_R4_FORMAT,PNTRX,STAT) CALL STFXMP(NN,D_R4_FORMAT,PNTRY,STAT) C C extract the indices of the complete line from input frame CALL PIXLIN(XA,YA,XB,YB,RM,MADRID(PNTRX),MADRID(PNTRY),NN,NINDX) IF (NINDX.GE.NN) + CALL STETER(3,'too many elements in index array...') C C create 1-dim output file + get z-values for pixels recorded in XINDX + YINDX CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + 1,NINDX,STARTC,STEPLN, + IDENTA,CUNITA,PNTRC,IMNOC,STAT) CALL ZIMA(MADRID(PNTRA),NPIXA,MADRID(PNTRX),MADRID(PNTRY),NINDX, + MADRID(PNTRC),CUTS(3),CUTS(4)) C C also save the Y-start coordinate STARTC(2) = STARTA(2) + (SUBLO(2)-1)*STEPA(2) CALL STDWRD(IMNOC,'START',STARTC,1,2,UNI,STAT) !STIPUT only fills 1 value C C common section C C copy non-standard descriptors, update LHCUTS + append HISTORY 8000 CALL DSCUPT(IMNOA,IMNOC,' ',STAT) CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT) C IF (REFSET.EQ.1) THEN !don't forget to update REFPIX CALL STDWRD(IMNOC,'REFPIX',DOUT,1,NAXISA,UNI,STAT) ENDIF C IF (EQUAL(1:).NE.' ') THEN CALL STFCLO(IMNOA,STAT) CALL STFCLO(IMNOC,STAT) NN = INDEX(EQUAL,' ') - 1 IF (NN.LT.1) NN = LEN(EQUAL) IDENTA(1:) = ' ' !IDENTA is at least 4 chars. longer IDENTA(1:) = EQUAL(1:) DO 8400 N=NN,1,-1 IF (EQUAL(N:N).EQ.'.') GOTO 8800 IF ((EQUAL(N:N).EQ.'/') .OR. + (EQUAL(N:N).EQ.']')) GOTO 8440 8400 CONTINUE 8440 IDENTA(NN+1:) = '.bdf' 8800 CALL STFRNM(FRAMEC,IDENTA,STAT) ENDIF C RETURN END SUBROUTINE SUBEZT C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBEZT version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frames C C.PURPOSE C extract only "referenced" pixels from a frame C C.ALGORITHM C Pull out only those pixels from inframe which are indicated in C a reference frame. Result frame will be a 1-dim frame. C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/80 input frame C IN_B/C/1/80 reference frame C OUT_A/C/1/80 result frame C THR/R/1/1 threshold value in reference frame C C.VERSIONS C 1.0 from version 1.00 as of 860110 C move to FORTRAN 77 + new ST-Interfaces C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,N,NAXISA,NAXISB INTEGER*8 PNTRA,PNTRB,PNTRC INTEGER IMNOA,IMNOB,STAT INTEGER SIZE,TRUESZ INTEGER NPIXA(3),NPIXB(3) INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*80 FRAMEA,FRAMEB,FRAMEC CHARACTER IDENTA*72,IDENTB*72,CUNITA*64,CUNITB*64 CHARACTER CBUF*100 C DOUBLE PRECISION STARTA(3),STARTB(3),STEPA(3),STEPB(3) C REAL THR C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA NPIXA /1,1,1/, NPIXB /1,1,1/ DATA STARTA /3*0.D0/, STEPA /3*1.D0/ DATA STARTB /3*0.D0/, STEPB /3*1.D0/ DATA IDENTA /' '/, CUNITA /' '/ DATA IDENTB /' '/, CUNITB /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) CALL STKRDC('IN_B',1,1,80,IAV,FRAMEB,UNI,NULO,STAT) CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) CALL STKRDR('THR',1,1,IAV,THR,UNI,NULO,STAT) C C map input + reference CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 3,NAXISA,NPIXA,STARTA,STEPA, + IDENTA,CUNITA,PNTRA,IMNOA,STAT) IF (FRAMEA.NE.FRAMEB) THEN CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 3,NAXISB,NPIXB,STARTB,STEPB, + IDENTB,CUNITB,PNTRB,IMNOB,STAT) ELSE IMNOB = IMNOA PNTRB = PNTRA NAXISB = NAXISA DO 20 N=1,NAXISA NPIXB(N) = NPIXA(N) 20 CONTINUE ENDIF C C get total size + map intermediate frame of that size SIZE = 1 DO 200 N=1,NAXISA SIZE = SIZE * NPIXA(N) 200 CONTINUE TRUESZ = 1 DO 300 N=1,NAXISB TRUESZ = TRUESZ * NPIXB(N) 300 CONTINUE IF (TRUESZ.LT.SIZE) CALL STETER + (8,'refframe smaller than input frame - we abort ...') C CALL STFXMP(SIZE,D_R4_FORMAT,PNTRC,STAT) C C now do the actual operation CALL GXDOIT(MADRID(PNTRA),MADRID(PNTRB),MADRID(PNTRC), + THR,SIZE,TRUESZ) IF (TRUESZ.LE.0) + CALL STETER(1,'result frame empty...') C IF (IMNOA.NE.IMNOB) + CALL STFCLO(IMNOB,STAT) !we don't need the refframe anymore CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + 1,TRUESZ,STARTA,STEPA, + IDENTA,CUNITA,PNTRB,IMNOB,STAT) C CALL COPYF(MADRID(PNTRC),MADRID(PNTRB),TRUESZ) C C update descriptors CALL DSCUPT(IMNOA,IMNOB,' ',STAT) N = INDEX(FRAMEC//' ',' ') - 1 WRITE(CBUF,10000) FRAMEC(1:N),TRUESZ CALL STTPUT(CBUF,STAT) C C that's it folks RETURN C 10000 FORMAT('frame ',A,' created with',I5,' pixels ...') END SUBROUTINE SUBINS C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBINS version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frames C C.PURPOSE C insert a frame into another one C C.ALGORITHM C straight forward C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/80 input frame C OUT_A/C/1/80 result frame C P3/C/1/80 string with start pixels C C.VERSIONS C 1.00 from version 1.20 as of 860109 C move to FORTRAN 77 + new ST Interfaces C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,LSUB,N,NAXISA,NAXISC INTEGER*8 PNTRX,PNTRY INTEGER NPIXA(3),NPIXC(3),IMNOA,IMNOC,STAT INTEGER SUBLO(3),IAX,NN INTEGER UNI(1),NULO,MADRID(1) INTEGER NPXY,K2,RMAIND,CHUCHO,CHUCHI,IDIM(3) INTEGER NOLIN,TRUSIZ,CHUNKI,CHUNKO,K1 INTEGER LOFF,KOFF,MOFF,ZOFF,LOOPLM,NPA(3),PXOFA(3) C CHARACTER*80 FRAMEA,FRAMEC CHARACTER PARM*80,OUTPUT*80 CHARACTER SUBSTR*30,ERROR3*40 C DOUBLE PRECISION STARTA(3),STARTC(3),STEPA(3),STEPC(3) DOUBLE PRECISION EPS,STA(3) C REAL CUTS(4),CUTVLS(2),RBUF C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA ERROR3 /'max. 3-dim images supported... '/ DATA SUBLO /1,1,1/, IDIM /1,1,1/ DATA PARM /' '/ DATA NPIXA /1,1,1/, NPIXC /1,1,1/ DATA NPA /1,1,1/, PXOFA /0,0,0/ DATA CUTVLS /+99999.,-99999./ DATA STARTA /3*0.0D0/, STEPA /3*1.0D0/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C open son + father frame (or daughter + mother if you prefer...) CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) C CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT) !son CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXISA,UNI,NULO,STAT) IF (NAXISA.GT.3) CALL STETER(44,ERROR3) CALL STDRDI(IMNOA,'NPIX',1,NAXISA,IAV,NPIXA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'START',1,NAXISA,IAV,STARTA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'STEP',1,NAXISA,IAV,STEPA,UNI,NULO,STAT) C CALL STFOPN(FRAMEC,D_R4_FORMAT,0,F_IMA_TYPE,IMNOC,STAT) !father CALL STDRDI(IMNOC,'NAXIS',1,1,IAV,NAXISC,UNI,NULO,STAT) IF (NAXISC.GT.3) CALL STETER(44,ERROR3) CALL STDRDI(IMNOC,'NPIX',1,NAXISC,IAV,NPIXC,UNI,NULO,STAT) CALL STDRDD(IMNOC,'START',1,NAXISC,IAV,STARTC,UNI,NULO,STAT) CALL STDRDD(IMNOC,'STEP',1,NAXISC,IAV,STEPC,UNI,NULO,STAT) C C if son has less dimensions, copy stuff from father IF (NAXISA.LT.NAXISC) THEN DO 300, N=NAXISA+1,NAXISC NPIXA(N) = 1 STARTA(N) = STARTC(N) STEPA(N) = STEPC(N) 300 CONTINUE NAXISA = NAXISC ENDIF C C stepsizes of son should not differ from father DO 400, N=1,NAXISA EPS = ABS(0.001 * STEPA(N)) IF (ABS(STEPA(N)-STEPC(N)).GT.EPS) THEN CALL STTPUT + ('Warning: stepsizes do not match...',STAT) GOTO 1000 ENDIF 400 CONTINUE C C get start of insertion in FRAMEC 1000 DO 1020 N=1,NAXISC STA(N) = STARTC(N) 1020 CONTINUE C CALL STKRDC('P3',1,1,80,IAV,PARM,UNI,NULO,STAT) IF (PARM(1:1).NE.'?') THEN IAX = 1 DO 1100 N=1,NAXISC CALL EXTRSS(PARM,',',IAX,SUBSTR,LSUB) IF (LSUB.LE.0) GOTO 1200 IF (SUBSTR(1:1).EQ.'@') THEN CALL GENCNV(SUBSTR(2:),1,1,SUBLO(N),RBUF,STA(N),IAV) IF (IAV.LE.0) + CALL STETER(46,'Invalid syntax in start coords. ') STA(N) = STARTC(N) + ((SUBLO(N)-1)*STEPC(N)) ELSE IF (SUBSTR(1:1).EQ.'<') THEN STA(N) = STARTC(N) ELSE IF (SUBSTR(1:1).EQ.'>') THEN STA(N) = STARTC(N) + ((NPIXC(N)-1)*STEPC(N)) ELSE CALL GENCNV(SUBSTR(1:),4,1,IAV,RBUF,STA(N),IAV) IF (IAV.LE.0) + CALL STETER(46,'Invalid syntax in start coords. ') ENDIF 1100 CONTINUE C C if no start values given, use descriptor START of FRAMEA ELSE DO 1150 N=1,NAXISA !note, NAXISA = NAXISC STA(N) = STARTA(N) 1150 CONTINUE ENDIF C C find offsets inside FRAMEA 1200 DO 1300 N=1,NAXISA IF (ABS(STEPC(N)).LT.1.E-20) THEN WRITE(OUTPUT,10123) N,STEPC(N) CALL STTPUT(OUTPUT,STAT) SUBLO(N) = 0 ELSE SUBLO(N) = NINT( (STA(N) - STARTC(N))/STEPC(N) ) ENDIF IF (SUBLO(N).GE.0) THEN SUBLO(N) = SUBLO(N) + 1 ELSE PXOFA(N) = - SUBLO(N) SUBLO(N) = 1 ENDIF NPA(N) = NPIXA(N) - PXOFA(N) 1300 CONTINUE C C check, if subframe overlaps with result frame DO 2500 N=1,NAXISA IF (SUBLO(N).GT.NPIXC(N)) + CALL STETER(3, + 'invalid coords.: sub- and modframe do not overlap...') NN = SUBLO(N) + NPA(N) - 1 !check for overflow IF ( NN.GT.NPIXC(N) ) THEN IDIM(N) = NPIXC(N) - SUBLO(N) + 1 CALL STTPUT + ('Warning: son does not fit completely in father...',STAT) ELSE IDIM(N) = NPA(N) ENDIF 2500 CONTINUE C C get buffersize of at least one line but no more than total size in y TRUSIZ = NPIXA(1)*IDIM(2) IAX = 20000 !first guess ... NN = IAX IF (NN.LT.NPIXA(1)) NN = NPIXA(1) IF (NN.GT.TRUSIZ) NN = TRUSIZ NOLIN = NN / NPIXA(1) !buffer = NOLIN * NPIXA(1) IF (NOLIN.GT.IDIM(2)) NOLIN = IDIM(2) !avoid overflow in y ... C 3500 CHUNKI = NOLIN * NPIXA(1) CHUNKO = NOLIN * NPIXC(1) IF ((NOLIN.GT.1).AND.(CHUNKO.GT.IAX*2)) THEN NOLIN = NOLIN - 1 GOTO 3500 ENDIF LOOPLM = IDIM(2) / NOLIN !get no. of passes through loop RMAIND = IDIM(2) - (LOOPLM * NOLIN) C CALL STFXMP(CHUNKI,D_R4_FORMAT,PNTRX,STAT) CALL STFXMP(CHUNKO,D_R4_FORMAT,PNTRY,STAT) C C now we loop through input (son) frame NPXY = NPIXC(1) * NPIXC(2) K2 = (SUBLO(3)-1) * NPXY !skip unused planes K1 = (SUBLO(2)-1) * NPIXC(1) + 1 !point to 1. used y-lines ZOFF = PXOFA(3)*NPIXA(2)*NPIXA(1) LOFF = PXOFA(1) + 1 !first pixel in x-line TRUSIZ = NPIXA(1) - LOFF + 1 !no. of pixels per line to move C DO 5000 IAX=1,IDIM(3) !loop over all planes MOFF = ZOFF + (PXOFA(2)*NPIXA(1)) + 1 !pointer for input frame KOFF = K2 + K1 IF (LOOPLM.LE.0) GOTO 4400 C DO 4000 N=1,LOOPLM CALL STFGET(IMNOA,MOFF,CHUNKI,NN,MADRID(PNTRX),STAT) !son in CALL STFGET(IMNOC,KOFF,CHUNKO,IAV,MADRID(PNTRY),STAT) !father in CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),LOFF,SUBLO(1), + TRUSIZ,NOLIN,NPIXA(1),NPIXC(1),CUTVLS) !rearrange data CALL STFPUT(IMNOC,KOFF,IAV,MADRID(PNTRY),STAT) !father out MOFF = MOFF + CHUNKI KOFF = KOFF + CHUNKO !follow with indices... 4000 CONTINUE C 4400 IF (RMAIND.GT.0) THEN CHUCHI = RMAIND * NPIXA(1) CHUCHO = RMAIND * NPIXC(1) CALL STFGET(IMNOA,MOFF,CHUCHI,NN,MADRID(PNTRX),STAT) CALL STFGET(IMNOC,KOFF,CHUCHO,IAV,MADRID(PNTRY),STAT) CALL GXMOVIT(MADRID(PNTRX),MADRID(PNTRY),LOFF,SUBLO(1), + TRUSIZ,RMAIND,NPIXA(1),NPIXC(1),CUTVLS) CALL STFPUT(IMNOC,KOFF,IAV,MADRID(PNTRY),STAT) MOFF = MOFF + CHUCHI ENDIF C K2 = K2 + NPXY !move to next plane ZOFF = ZOFF + (NPIXA(1)*NPIXA(2)) 5000 CONTINUE C C do not copy all descriptors, just update descr HISTORY CALL DSCUPT(IMNOC,IMNOC,' ',STAT) C C and LHCUTS also CALL STDRDR(IMNOC,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT) IAX = 0 IF (CUTVLS(1).LT.CUTS(3)) THEN CUTS(3) = CUTVLS(1) IAX = 1 ENDIF IF (CUTVLS(2).GT.CUTS(4)) THEN CUTS(4) = CUTVLS(2) IAX = 1 ENDIF IF (IAX.EQ.1) CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT) C C that's it folks RETURN C 10123 FORMAT('Warning: step(',I1,') of father frame = ',G10.5, + ' result may be corrupted ...') C END SUBROUTINE SUBMAP C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBMAP version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frame, mapping C C.PURPOSE C map an image via another one (which acts like a LUT) C C.ALGORITHM C straight forward C C.INPUT/OUTPUT C the following keywords are used: C C IN_A/C/1/80 name of input frame C IN_B/C/1/80 name of map frame C OUT_A/C/1/80 name of result frame C INPUTI/I/1/2 map flag: C (1) = 0, indices outside map space go to C start + end C (1) = 1, indices outside avoid mapping, C i.e output pixel = input pixel C (2) = 0, map frame is 1-dim and has C equidistant coord space C (2) = 1, map frame is 2-dim and has C pixel coords stored in 1. line, C pixel intensities in 2. line C C.VERSIONS C C 1.00 from version 1.00 s of 871112 C move to FORTRAN 77 + new ST Interfaces C C-------------------------------------------------- C IMPLICIT NONE C INTEGER N,IAV,NAXISA,NAXISB INTEGER*8 PNTRA,PNTRB,PNTRC INTEGER IMNOA,IMNOB,IMNOC,STAT INTEGER NPIXA(3),NPIXB(3),KMAP(2) INTEGER UNI(1),NULO,MADRID(1) C DOUBLE PRECISION STARTA(3),STEPA(3),STARTB(3),STEPB(3) C CHARACTER CUNITA*64,IDENTA*72,CUNITB*10,IDENTB*10 CHARACTER FRAMEA*80,FRAMEB*80,FRAMEC*80 C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA IDENTA /' '/, CUNITA /' '/ C C get name of input image, map image + output image CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) CALL STKRDC('IN_B',1,1,80,IAV,FRAMEB,UNI,NULO,STAT) CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) CALL STKRDI('INPUTI',1,2,IAV,KMAP,UNI,NULO,STAT) C DO 400 N=1,3 NPIXA(N) = 1 STARTA(N) = 0.D0 STEPA(N) = 1.D0 NPIXB(N) = 1 STARTB(N) = 0.D0 STEPB(N) = 1.D0 400 CONTINUE C C map input frame + map frame for reading CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXISA,NPIXA,STARTA,STEPA,IDENTA, + CUNITA,PNTRA,IMNOA,STAT) CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXISB,NPIXB,STARTB,STEPB,IDENTB, + CUNITB,PNTRB,IMNOB,STAT) IF (KMAP(2).GT.0) THEN IF (NPIXB(2).LE.1) + CALL STETER(1,'we need 2-dim map for this option...') ENDIF C CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISA,NPIXA,STARTA,STEPA,IDENTA, + CUNITA,PNTRC,IMNOC,STAT) C C execute algorithm in a subroutine CALL GXIZAMAP(NPIXA,NPIXB,STARTB,STEPB, + MADRID(PNTRA),MADRID(PNTRB),MADRID(PNTRC),KMAP) C C that's it folks... RETURN END SUBROUTINE SUBMAT C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBMAT version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C bulk data frame, matrix ops C C.PURPOSE C do the usual matrix operations on an image C C.ALGORITHM C straight forward C C.INPUT/OUTPUT C the following keywords are used: C C IN_A/C/1/80 name of input frame C ACTION/C/3/2 action C = 'TR ', transpose a matrix C around major diagonal C = 'LC ', transpose a matrix C around minor diagonal C i.e. change lines into columns C OUT_A/C/1/80 name of result frame C INPUTI/I/1/2 xsize,ysize for 'LC' C C.VERSIONS C 1.00 from version 1.20 as of 870105 C C-------------------------------------------------- C IMPLICIT NONE C INTEGER*8 PNTRA,PNTRC INTEGER IAV,NAXIS,STAT INTEGER IMNOA,IMNOC INTEGER NPIX(2),NPIXC(2),ISIZE(2) INTEGER UNI(1),NULO,MADRID(1) C DOUBLE PRECISION START(2),STEP(2),DTMP C CHARACTER CUNIT*48,IDENT*72 CHARACTER FRAMEA*80,FRAMEC*80,ACTION*2 C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA IDENT /' '/, CUNIT /' '/ C C get name of input, ouput image + action flags CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) CALL STKRDC('ACTION',1,3,2,IAV,ACTION,UNI,NULO,STAT) CALL UPCAS(ACTION,ACTION) C C map input frame CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTRA,IMNOA,STAT) C C check dimensions... IF (NAXIS.NE.2) + CALL STETER(1,'input must be a 2-dim frame ...') C C exchange NPIX, adjust START, STEP + map output frame NPIXC(1) = NPIX(2) NPIXC(2) = NPIX(1) IF (ACTION.EQ.'TR') THEN DTMP = START(1) + (NPIX(1)-1) * STEP(1) !end in X START(1) = START(2) + (NPIX(2)-1) * STEP(2) !end in Y START(2) = DTMP DTMP = - STEP(1) STEP(1) = - STEP(2) STEP(2) = DTMP ELSE DTMP = START(1) START(1) = START(2) START(2) = DTMP DTMP = STEP(1) STEP(1) = STEP(2) STEP(2) = DTMP ENDIF CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIXC,START,STEP,IDENT, + CUNIT,PNTRC,IMNOC,STAT) C C execute algorithm in a subroutine IF (ACTION.EQ.'TR') THEN CALL GXMATRX(MADRID(PNTRA),NPIX,MADRID(PNTRC)) ELSE CALL STKRDI('INPUTI',1,2,IAV,ISIZE,UNI,NULO,STAT) CALL LINCOL(MADRID(PNTRA),NPIX,ISIZE,MADRID(PNTRC)) ENDIF CALL DSCUPT(IMNOA,IMNOC,' ',STAT) C C that's it folks... RETURN END