C @(#)genyy1.for 13.1.1.2 (ESO-DMD) 06/23/98 11:57:42 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 GENYY1 C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.IDENTIFICATION C program GENYY1 951012 C K. Banse ESO - Garching 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 = GR, for GROW/IMAGE C = IM, for converting from image to table C = TB, for converting from table to image C = CC, for aligning center pixels C = SC, for SCALE C = WE, for WEIGHT C = RG, for ROTGEN C = RO, for ROTATE C C.VERSION C see SCCS C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,STAT INTEGER UNI(1),NULO C CHARACTER ACTION*3 C C set up MIDAS environment + enable automatic error abort CALL STSPRO('GENYY1') C C get option CALL STKRDC('ACTION',1,1,3,IAV,ACTION,UNI,NULO,STAT) C IF (ACTION(1:2).EQ.'GR') THEN CALL GROW(ACTION(3:3)) !grow/image ELSE IF (ACTION(1:2).EQ.'IM') THEN CALL IM2TAB(1) !image -> table ELSE IF (ACTION(1:2).EQ.'TB') THEN CALL IM2TAB(2) !table -> image ELSE IF (ACTION(1:2).EQ.'CC') THEN CALL ALICNT ELSE IF (ACTION(1:2).EQ.'SC') THEN CALL SUBSCL ELSE IF (ACTION(1:2).EQ.'WE') THEN CALL SUBWEI ELSE IF (ACTION(1:2).EQ.'RG') THEN CALL SUBRTG ELSE IF (ACTION(1:2).EQ.'RO') THEN CALL SUBROT ENDIF C C finished CALL STSEPI END SUBROUTINE GROW(OPT) C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine GROW version 1.00 871217 C J.D. Ponz ESO - Garching C C.KEYWORDS C image expansion, pixel replication C C.PURPOSE C expand input image by replicating a selected row row- or columnwise C C.ALGORITHM C straight forward C C.INPUT/OUTPUT C call as GROW(OPT) C C input par: C OPT: char. 'L' or 'C' for row replication as line or column C C.VERSION C-------------------------------------------------------------------- C IMPLICIT NONE C INTEGER MADRID INTEGER KUN,KNUL,IMNOA,IMNOB INTEGER IN1,KK,NAXISA,NAXISB,ACTVAL,SUBLO(3) INTEGER NPIXA(2),NPIXB(2),STATUS INTEGER*8 PNTRA,PNTRB,PNTRW C DOUBLE PRECISION STEPA(2),STEPB(2),DPAR(3) DOUBLE PRECISION STARTA(2),STARTB(2) C CHARACTER OPT*1 CHARACTER*60 FRAMEA,FRAMEB CHARACTER CUNIT*48,IDENT*72,CBUF*40,MOPT*2,NEWSTR*80 C COMMON /VMR/MADRID(1) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA IDENT /' '/, CUNIT /' '/ DATA NPIXA /1,1/ DATA STARTA /0.D0,0.D0/ DATA STEPA /1.D0,1.D0/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get input, output frame C CALL STKRDC('IN_A',1,1,60,ACTVAL,FRAMEA,KUN,KNUL,STATUS) CALL STKRDC('OUT_A',1,1,60,ACTVAL,FRAMEB,KUN,KNUL,STATUS) CALL GENEQF(FRAMEA,FRAMEB,ACTVAL) IF (ACTVAL.NE.0) + CALL STETER(9,'inframe and outframe are equal...') C CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXISA,NPIXA,STARTA,STEPA, + IDENT,CUNIT,PNTRA,IMNOA,STATUS) CALL STKRDD('INPUTD',1,3,ACTVAL,DPAR,KUN,KNUL,STATUS) C C get line/column we want to replicate C CALL STKRDC('P5',1,1,40,ACTVAL,CBUF,KUN,KNUL,STATUS) MOPT(2:2) = 'L' !default to line input IF (NAXISA.EQ.1) THEN IN1 = 1 ELSE KK = INDEX(CBUF,',') IF (KK.GT.1) THEN KK = KK + 1 IF ((CBUF(KK:KK).EQ.'c') .OR. (CBUF(KK:KK).EQ.'C')) + MOPT(2:2) = 'C' ENDIF IF (MOPT(2:2).EQ.'C') THEN NEWSTR(1:) = '[ ' NEWSTR(2:) = CBUF KK = INDEX(NEWSTR,' ') NEWSTR(KK:) = ',@1] ' CALL EXTCO1(IMNOA,NEWSTR,2,KK,SUBLO,STATUS) IF (STATUS.NE.0) CALL STSEPI IN1 = SUBLO(1) ELSE NEWSTR(1:) = '[@1, ' NEWSTR(5:) = CBUF KK = INDEX(NEWSTR,' ') NEWSTR(KK:) = '] ' CALL EXTCO1(IMNOA,NEWSTR,2,KK,SUBLO,STATUS) IF (STATUS.NE.0) CALL STSEPI IN1 = SUBLO(2) ENDIF ENDIF C NAXISB = 2 ACTVAL = NINT(DPAR(3)) IF (ACTVAL.LE.0) ACTVAL = NPIXA(1) C C if we extract a column, some more work necessary IF (MOPT(2:2).EQ.'C') THEN CALL STFXMP(NPIXA(2),D_R4_FORMAT,PNTRW,STATUS) CALL GETCOL(MADRID(PNTRA),NPIXA,IN1,MADRID(PNTRW)) IN1 = 1 NPIXA(1) = NPIXA(2) !simulate 1-dim frame input STARTA(1) = STARTA(2) STEPA(1) = STEPA(2) PNTRA = PNTRW !and point to extracted column ENDIF C C map output frame C IF ((OPT.EQ.'c') .OR. (OPT.EQ.'C')) THEN NPIXB(1) = ACTVAL !we replicate the input row as column NPIXB(2) = NPIXA(1) STARTB(1) = DPAR(1) STARTB(2) = STARTA(1) STEPB(1) = DPAR(2) STEPB(2) = STEPA(1) CBUF(1:16) = CUNIT(17:32) CUNIT(17:) = 'none ' CUNIT(33:48) = CBUF(1:16) MOPT(1:1) = 'C' ELSE !we replicate the input row as row NPIXB(1) = NPIXA(1) NPIXB(2) = ACTVAL STARTB(1) = STARTA(1) STARTB(2) = DPAR(1) STEPB(1) = STEPA(1) STEPB(2) = DPAR(2) CUNIT(33:) = 'none ' MOPT(1:1) = 'L' ENDIF C CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISB,NPIXB,STARTB,STEPB, + IDENT,CUNIT,PNTRB,IMNOB,STATUS) C C and do it CALL GROWIT + (MOPT(1:1),NPIXA,MADRID(PNTRA),NPIXB,MADRID(PNTRB),IN1) CALL STDCOP(IMNOA,IMNOB,4,'LHCUTS',STATUS) CALL DSCUPT(IMNOA,IMNOB,' ',STATUS) C C return RETURN END SUBROUTINE GETCOL(A,NPIXA,NO,B) C C IMPLICIT NONE C INTEGER NPIXA(2),NO INTEGER I,KIN C REAL A(*),B(*) C KIN = NO DO 150, I=1,NPIXA(2) B(I) = A(KIN) KIN = KIN + NPIXA(1) 150 CONTINUE C RETURN END SUBROUTINE IM2TAB(FLAG) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine IM2TAB version 1.00 891213 C K. Banse ESO - Garching C C.KEYWORDS C image, table conversion C C.PURPOSE C convert a table holding subareas of an image into a black+white mask C and vice versa C C.ALGORITHM C see below... C C.INPUT/OUTPUT C call as IM2TAB(FLAG) C C Input: C FLAG Integer 1 for imae -> table C 2 for table -> image C C.VERSION C 1.00 creation C C-------------------------------------------------- C IMPLICIT NONE C INTEGER MAXDIM PARAMETER (MAXDIM=4000) C INTEGER FLAG INTEGER MADRID INTEGER KUN,KNUL,TID,IMNO,IMNO2,MPIX(2) INTEGER KCOUNT,NAXIS,IAV,WSIZE INTEGER*8 PNTRW,PNTRV INTEGER NPIX(3),STAT INTEGER STAX(MAXDIM),STAY(MAXDIM),ENDX(MAXDIM),ENDY(MAXDIM) INTEGER TBCLNM(4),NCOLS C REAL THRESH,RBUF(12),FGR,BGR C DOUBLE PRECISION START(3),STEP(3),DBUF(4) C CHARACTER*60 FRAME,TABLE CHARACTER CUNIT*48,IDENT*80,MIDUNI*2,TMPFIL*14 CHARACTER TBUNIT*16,TBLABL(4)*16 C LOGICAL TABNUL(4),SELFLG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C COMMON /VMR/MADRID(1) C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA NPIX /1,1,1/, START /0.,0.,0./, STEP /1.,1.,1./ DATA TBUNIT /'WORLD COORD '/ DATA TBLABL /'XSTART ','YSTART ','XEND ','YEND '/ DATA TBCLNM /0,0,0,0/ C C get image frame, table CALL STKRDC('IN_A',1,1,60,IAV,FRAME,KUN,KNUL,STAT) CALL STKRDC('IN_B',1,1,60,IAV,TABLE,KUN,KNUL,STAT) NCOLS = 4 CALL STKRDC('MID$SESS',1,11,2,IAV,MIDUNI,KUN,KNUL,STAT) TMPFIL = 'imtab'//MIDUNI//'dum.bdf' C C branch according to FLAG IF (FLAG.EQ.2) GOTO 5000 C C here we handle building of table from image C CALL STKRDR('INPUTR',1,1,IAV,THRESH,KUN,KNUL,STAT) C C open data file CALL STFOPN(FRAME,D_R4_FORMAT,0,F_IMA_TYPE,IMNO,STAT) CALL STDRDI(IMNO,'NAXIS',1,1,IAV,NAXIS,KUN,KNUL,STAT) CALL STDRDI(IMNO,'NPIX',1,NAXIS,IAV,NPIX,KUN,KNUL,STAT) CALL STDRDD(IMNO,'START',1,NAXIS,IAV,START,KUN,KNUL,STAT) CALL STDRDD(IMNO,'STEP',1,NAXIS,IAV,STEP,KUN,KNUL,STAT) C C map working space for 50/100 lines and 3 working lines IF (NPIX(1).GT.512) THEN WSIZE = NPIX(1) * 50 ELSE WSIZE = NPIX(1) * 100 ENDIF CALL STFXMP(WSIZE,D_R4_FORMAT,PNTRW,STAT) KCOUNT = 3 * (NPIX(1) + 2) !1 more pixel in front and the end CALL STFXMP(KCOUNT,D_I4_FORMAT,PNTRV,STAT) C MPIX(1) = NPIX(1) + 2 MPIX(2) = NPIX(2) + 2 KCOUNT = MPIX(1) * MPIX(2) CALL STFCRE(TMPFIL,D_I4_FORMAT,F_O_MODE,F_IMA_TYPE, + KCOUNT,IMNO2,STAT) CALL STDWRI(IMNO2,'NAXIS',NAXIS,1,1,KUN,STAT) CALL STDWRI(IMNO2,'NPIX',MPIX,1,2,KUN,STAT) CALL STDWRD(IMNO2,'START',START,1,2,KUN,STAT) CALL STDWRD(IMNO2,'STEP',STEP,1,2,KUN,STAT) C KCOUNT = 0 CALL DOI2T(IMNO,NPIX,MADRID(PNTRW),WSIZE,MADRID(PNTRV), + THRESH,STAX,STAY,ENDX,ENDY,MAXDIM,KCOUNT,IMNO2) C C now create the table and fill it CALL TBTINI(TABLE,0,F_O_MODE,NCOLS,KCOUNT,TID,STAT) DO 4000, IAV=1,NCOLS CALL TBCINI(TID,D_R4_FORMAT,1,'G12.6',TBUNIT, + TBLABL(IAV),TBCLNM(IAV),STAT) 4000 CONTINUE C DO 4400, IAV=1,KCOUNT RBUF(1) = START(1) + ((STAX(IAV) - 1) * STEP(1)) RBUF(2) = START(2) + ((STAY(IAV) - 1) * STEP(2)) RBUF(3) = START(1) + ((ENDX(IAV) - 1) * STEP(1)) RBUF(4) = START(2) + ((ENDY(IAV) - 1) * STEP(2)) CALL TBRWRR(TID,IAV,NCOLS,TBCLNM,RBUF,STAT) 4400 CONTINUE C CALL TBTCLO(TID,STAT) WRITE(IDENT,10000) KCOUNT CALL STTPUT(IDENT,STAT) RETURN C C C here we handle building of image from table C 5000 CALL STKRDR('INPUTR',1,2,IAV,RBUF,KUN,KNUL,STAT) !bgr,fgr CALL STKRDD('INPUTD',1,4,IAV,DBUF,KUN,KNUL,STAT) !start,step CALL STKRDI('INPUTI',1,2,IAV,NPIX,KUN,KNUL,STAT) !NPIX of image BGR = RBUF(1) FGR = RBUF(2) START(1) = DBUF(1) START(2) = DBUF(2) STEP(1) = DBUF(3) STEP(2) = DBUF(4) IDENT(1:) = 'Image created from table ' CUNIT(1:) = ' ' NAXIS = 2 C CALL STIPUT(FRAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START,STEP, + IDENT,CUNIT,PNTRW,IMNO,STAT) KCOUNT = NPIX(1) * NPIX(2) CALL CONFIL(BGR,MADRID(PNTRW),KCOUNT) !fill frame with background C CALL TBTOPN(TABLE,F_I_MODE,TID,STAT) DO 6000, IAV=1,NCOLS CALL TBLSER(TID,TBLABL(IAV),TBCLNM(IAV),STAT) IF (TBCLNM(IAV).LE.0) CALL STETER + (9,'column not found in input table...') 6000 CONTINUE CALL TBIGET(TID,IAV,KCOUNT,IAV,IAV,IAV,STAT) !get no. of rows C DO 6600, IAV=1,KCOUNT CALL TBSGET(TID,IAV,SELFLG,STAT) !only use selected rows IF (.NOT.SELFLG) GOTO 6600 C CALL TBRRDR(TID,IAV,4,TBCLNM,RBUF,TABNUL,STAT) STAX(1) = (RBUF(1) - START(1)) / STEP(1) + 1 STAX(2) = (RBUF(2) - START(2)) / STEP(2) + 1 ENDX(1) = (RBUF(3) - START(1)) / STEP(1) + 1 ENDX(2) = (RBUF(4) - START(2)) / STEP(2) + 1 ENDY(1) = ENDX(1) - STAX(1) IF (ENDY(1).LT.0) ENDY(1) = -ENDY(1) ENDY(1) = ENDY(1) + 1 ENDY(2) = ENDX(2) - STAX(2) IF (ENDY(2).LT.0) ENDY(2) = -ENDY(2) ENDY(2) = ENDY(2) + 1 CALL COPYF2(FGR,MADRID(PNTRW),NPIX,STAX,ENDY) 6600 CONTINUE CALL TBTCLO(TID,STAT) WRITE(IDENT,10001) KCOUNT CALL STTPUT(IDENT,STAT) C RETURN C 10000 FORMAT(I6,' entries written into output table ...') 10001 FORMAT(I6,' entries processed from input table ...') END SUBROUTINE DOI2T(IMNO,NPIX,W,WSIZE,IV,THR, + STAX,STAY,ENDX,ENDY,MAXDIM,KCOUNT,IM2) C IMPLICIT NONE C INTEGER IMNO,NPIX(*),WSIZE,STAX(*),STAY(*),ENDX(*),ENDY(*) INTEGER IV(1500),MAXDIM,KCOUNT,IM2 INTEGER XDIM,YDIM,YDIM2,N,NX,NY,K,FELEM,ACTSZ,STAT INTEGER VFIRST,VLAST,NV,VDIM,WOFF,IAV INTEGER FOUT C REAL W(*),THR C XDIM = NPIX(1) YDIM = NPIX(2) YDIM2 = YDIM - 2 VDIM = XDIM + 2 C FOUT = 1 C FELEM = 1 !get first chunk of data CALL STFGET(IMNO,FELEM,WSIZE,IAV,W,STAT) FELEM = FELEM + IAV ACTSZ = IAV C VFIRST = 2*VDIM + 1 VLAST = 3 * VDIM C C fill first 3 working lines with zero DO 500, K=1,VLAST IV(K) = 0 500 CONTINUE C N = VDIM + 1 DO 700, K=1,XDIM !fill lines 2, 3 IF (W(K).GE.THR) IV(K+N) = 1 IF (W(K+XDIM).GE.THR) IV(K+VFIRST) = 1 700 CONTINUE WOFF = 2 * XDIM C CALL STFPUT(IM2,FOUT,VDIM,IV,STAT) !write 1. line of temp file FOUT = FOUT + VDIM C C move line by line DO 5000, NY=1,YDIM C DO 3000, NV=2,VDIM-1 N = NV + VDIM !move through middle work line C IF (IV(N).EQ.1) THEN C C test if a candidate for lower left start pixel IF (IV(N-1).EQ.0) THEN IAV = N - VDIM IF ((IV(IAV).EQ.0) .OR. + ((IV(IAV).EQ.1).AND.(IV(IAV+1).EQ.0))) THEN KCOUNT = KCOUNT + 1 C IF (KCOUNT.GT.MAXDIM) + CALL STETER(1, + 'Internal data overflow - we abort...') C NX = NV - 1 STAX(KCOUNT) = NX STAY(KCOUNT) = NY ENDX(KCOUNT) = 0 ENDY(KCOUNT) = 0 ENDIF ENDIF C C test if a candidate for lower right end pixel IF (IV(N+1).EQ.0) THEN IAV = N - VDIM IF ((IV(IAV).EQ.0) .OR. + ((IV(IAV).EQ.1).AND.(IV(IAV-1).EQ.0))) THEN ENDX(KCOUNT) = NV - 1 !last entry... ENDIF C C test if a candidate for upper right end pixel IAV = N + VDIM IF ((IV(IAV).EQ.0) .OR. + ((IV(IAV).EQ.1).AND.(IV(IAV-1).EQ.0))) THEN NX = NV - 1 DO 2000, K=KCOUNT,1,-1 !search for match in NX IF (ENDX(K).EQ.NX) THEN ENDY(K) = NY GOTO 3000 ENDIF 2000 CONTINUE ENDIF ENDIF ENDIF 3000 CONTINUE C C update working lines DO 3300, K=1,2*VDIM !copy line 2,3 -> line 1,2 IV(K) = IV(K+VDIM) 3300 CONTINUE DO 3330, K=VFIRST,VLAST !and set line 3 to zero IV(K) = 0 3330 CONTINUE C IF (NY .LE. YDIM2) THEN DO 3400, K=1,XDIM IF (W(K+WOFF).GE.THR) IV(VFIRST+K) = 1 3400 CONTINUE WOFF = WOFF + XDIM ENDIF C C test, if we need more data IF (NY .LT. YDIM2) THEN IF (WOFF.GE.ACTSZ) THEN CALL STFGET(IMNO,FELEM,WSIZE,IAV,W,STAT) FELEM = FELEM + IAV ACTSZ = IAV WOFF = 0 !reset pointer in data buffer ENDIF ENDIF C CALL STFPUT(IM2,FOUT,VDIM,IV,STAT) FOUT = FOUT + VDIM C 5000 CONTINUE C CALL STFPUT(IM2,FOUT,VDIM,IV(VDIM+1),STAT) !last middle line C RETURN END SUBROUTINE ALICNT C C-------------------------------------------------------------------- C IMPLICIT NONE C INTEGER MADRID INTEGER KUN,KNUL,IMNOA,IMNOB INTEGER NAXISA,NAXISB,ACTVAL INTEGER N,SCPIX(3),FCPIX(3),DIFPIX INTEGER NPIXA(3),NPIXB(3),STATUS C DOUBLE PRECISION STEPA(3),STEPB(3) DOUBLE PRECISION STARTA(3),STARTB(3),STARTC(3) C CHARACTER*60 FRAMEA,FRAMEB CHARACTER LINE*80 C COMMON /VMR/MADRID(1) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA STARTC /-1.D0,-1.D0,-1.D0/ INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get son and father frame C CALL STKRDC('IN_A',1,1,60,ACTVAL,FRAMEA,KUN,KNUL,STATUS) CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STATUS) CALL STDRDI(IMNOA,'NAXIS',1,1,ACTVAL,NAXISA,KUN,KNUL,STATUS) CALL STDRDI(IMNOA,'NPIX',1,NAXISA,ACTVAL,NPIXA,KUN,KNUL,STATUS) CALL STDRDD(IMNOA,'START',1,NAXISA,ACTVAL,STARTA,KUN,KNUL,STATUS) CALL STDRDD(IMNOA,'STEP',1,NAXISA,ACTVAL,STEPA,KUN,KNUL,STATUS) C CALL STKRDC('IN_B',1,1,60,ACTVAL,FRAMEB,KUN,KNUL,STATUS) CALL STFOPN(FRAMEB,D_R4_FORMAT,0,F_IMA_TYPE,IMNOB,STATUS) CALL STDRDI(IMNOB,'NAXIS',1,1,ACTVAL,NAXISB,KUN,KNUL,STATUS) IF (NAXISB .LT. NAXISA) + CALL STETER(1,'Father must have at least dimension of son...') C CALL STDRDI(IMNOB,'NPIX',1,NAXISB,ACTVAL,NPIXB,KUN,KNUL,STATUS) CALL STDRDD(IMNOB,'START',1,NAXISB,ACTVAL,STARTB,KUN,KNUL,STATUS) CALL STDRDD(IMNOB,'STEP',1,NAXISB,ACTVAL,STEPB,KUN,KNUL,STATUS) C C get center pixels of son CALL STKRDC('P3',1,1,80,ACTVAL,LINE,KUN,KNUL,STATUS) CALL EXTCO1(IMNOA,LINE,3,N,SCPIX,STATUS) IF (STATUS .NE. 0) + CALL STETER(2,'invalid center coords of son...') C C get center pixels of father CALL STKRDC('P4',1,1,80,ACTVAL,LINE,KUN,KNUL,STATUS) CALL EXTCO1(IMNOB,LINE,3,N,FCPIX,STATUS) IF (STATUS .NE. 0) + CALL STETER(2,'invalid center coords of father...') C C now calculate the relevant offsets DO 3000, N=1,NAXISA DIFPIX = FCPIX(N) - SCPIX(N) STARTC(N) = STARTB(N) + (DIFPIX * STEPB(N)) 3000 CONTINUE C C show results IF (NAXISA.EQ.1) THEN WRITE(LINE,10000) STARTC(1) ELSE IF (NAXISA.EQ.2) THEN WRITE(LINE,10001) STARTC(1),STARTC(2) ELSE IF (NAXISA.GE.3) THEN WRITE(LINE,10002) STARTC(1),STARTC(2),STARTC(3) ENDIF CALL STTPUT(LINE,STATUS) CALL STKWRD('OUTPUTD',STARTC,1,3,KUN,STATUS) C RETURN C 10000 FORMAT('new start in x: ',G20.8) 10001 FORMAT('new start in x,y: ',2G20.8) 10002 FORMAT('new start in x,y,z: ',3G20.8) C END SUBROUTINE SUBSCL C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C subroutine SUBSCL version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS: C scaling C C.PURPOSE: C Scale input image into given interval. C C.ALGORITHM: C Same as in WIMG + IPACK C C.INPUT/OUTPUT: C the following keywords are used: C C IN_A/C/1/60 input image C OUT_A/C/1/60 output image C INPUTR/R/1/2 low,hi (inclusive) end of interval C C.VERSIONS C 1.00 from version 1.20 as of 870429 C use FORTRAN 77 + new ST interfaces C C------------------------------------------------------------------- C IMPLICIT NONE C INTEGER STAT,IAV,NAXIS,NPIX(3) INTEGER*8 PNTRA,PNTRB INTEGER IMNOA,IMNOB,SIZE INTEGER UNI(1),NULO,MADRID(1) C CHARACTER FRAMEA*60,FRAMEB*60 CHARACTER IDENT*72,CUNIT*64 C DOUBLE PRECISION START(3),STEP(3) C REAL RANGE,F,RIN(2),CUTS(4),RMAP C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA IDENT /' '/, CUNIT /' '/ C C get frames + scaling interval CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,UNI,NULO,STAT) CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEB,UNI,NULO,STAT) CALL STKRDR('INPUTR',1,2,IAV,RIN,UNI,NULO,STAT) C C map the frames CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 3,NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTRA,IMNOA,STAT) CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTRB,IMNOB,STAT) C C read cuts of input frame CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT) IF (CUTS(1).GE.CUTS(2)) THEN CUTS(1) = CUTS(3) CUTS(2) = CUTS(4) IF (CUTS(1).GE.CUTS(2)) + CALL STETER(1,'invalid cut values...') ENDIF C C get dynamic range RMAP = RIN(2) - RIN(1) IF (RMAP.LE.10.E-25) CALL STETER(2,'invalid mapping interval...') C RANGE = CUTS(2) - CUTS(1) IF (RANGE.LT.10.E-25) THEN F = 1.0 CUTS(1) = 0. ELSE F = RMAP / RANGE ENDIF C C now do it SIZE = NPIX(1) * NPIX(2) * NPIX(3) CALL SCALA(MADRID(PNTRA),MADRID(PNTRB),SIZE,F,CUTS,RIN) C C modify cut values CUTS(1) = 0. CUTS(2) = 0. CUTS(3) = RIN(1) CUTS(4) = RIN(2) CALL STDWRR(IMNOB,'LHCUTS',CUTS,1,4,UNI,STAT) C C That's it folks... RETURN END SUBROUTINE SCALA(A,B,NDIM,F,CUTS,RIN) C IMPLICIT NONE C INTEGER NDIM INTEGER N C REAL A(*),B(*),F,CUTS(4),RIN(2) C DO 1000, N=1,NDIM IF (A(N).LE.CUTS(1)) THEN B(N) = RIN(1) ELSE IF (A(N).GE.CUTS(2)) THEN B(N) = RIN(2) ELSE B(N) = (F * (A(N) - CUTS(1)) ) + RIN(1) ENDIF 1000 CONTINUE C RETURN END SUBROUTINE SUBWEI C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBWEI version 1.00 890123 C K. Banse ESO - Garching C 1.10 910116 C C.KEYWORDS C bulk data frames, statistical weights C C.PURPOSE C determine statistical weights for frames (e.g. for use in AVERAGE/WEIGHTS) C C.ALGORITHM C extract input frames either from given parameter string or from catalogue, C get mean + std. of "sky-field" of image frames, compute weights for all C frames by dividing each mean by max. mean, i.e. frame with max. sky C gets weight = 1. C C.INPUT/OUTPUT C the following keys are used: C C P1/C/1/80 list of frames to be averaged or catalog C INPUTC/C/1/80 CURSOR or xsta,ysta,xend,yend of window C C.VERSIONS C 1.00 created from version 3.00 as of 841121 C use FORTRAN 77 + new ST interfaces C 1.10 make it work also for 1-dim frames and pixel coords. C C-------------------------------------------------- C IMPLICIT NONE C INTEGER BEGIN INTEGER FRMCNT,IAV,LL,N,M,STAT,ISEQ INTEGER NAXIS,ITER,KODE,LEAP,ISTP,ILAST,MINP INTEGER NPIX(2),IMNOA INTEGER*8 PNTRA INTEGER APIX(2,2),NPNTS(5) INTEGER PIXEL(2,2),SUBLO(2),SUBHI(2),SUBDIM INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*60 FRAME(50) CHARACTER CUNIT*20,IDENT*20,OUTPUT*80 CHARACTER LINE*80,CATFIL*60 CHARACTER ERROR1*60,ERROR2*50,ERROR3*40,ERROR4*40 CHARACTER ERROR5*40 C DOUBLE PRECISION STEP(2),START(2) DOUBLE PRECISION OSTART(2),OEND(2) C REAL DIF,W,BAGR(50),SIGM(50),BMAX REAL BGVL(5),RMSV(5),COORDS(10) C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA ERROR1 /'operands do not match in stepsize... '/ DATA ERROR2 /'operands do not overlap...'/ DATA ERROR3 /'stepsizes of different sign...'/ DATA ERROR4 /'catalog empty...'/ DATA ERROR5 /'more than 50 entries in catalog ...'/ DATA APIX /4*1/, BMAX /-9.E10/ DATA SUBLO /1,1/, SUBHI /1,1/ DATA IDENT /' '/, CUNIT /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up stuff for BGVAL KODE = 1 MINP = 2 ISTP = 0 LEAP = 0 ITER = 0 ISEQ = 0 C C get necessary keys CALL STKRDC('P1',1,1,80,IAV,LINE,UNI,NULO,STAT) !list of frames... C C test, if we have list of frames or catalog LL = INDEX(LINE,'.cat') IF (LL.LE.0) THEN LL = INDEX(LINE,'.CAT') IF (LL.LE.0) GOTO 500 !No. Treat list input... ENDIF C C here we handle input from a catalog C get name of first image file in catalog C CATFIL(1:) = LINE(1:) CALL STCGET(CATFIL,0,FRAME(1),OUTPUT,ISEQ,STAT) IF (FRAME(1)(1:1).EQ.' ') !there must be at least 1 image ... + CALL STETER(4,ERROR4) C C now loop through catalogue DO 200, N=2,50 !max. 50 frames CALL STCGET(CATFIL,0,FRAME(N),OUTPUT,ISEQ,STAT) IF (FRAME(N)(1:1).EQ.' ') THEN !indicates the end ... FRMCNT = N - 1 GOTO 1000 ENDIF 200 CONTINUE C CALL STCGET(CATFIL,0,LINE,OUTPUT,ISEQ,STAT) IF (LINE(1:1).NE.' ') + CALL STTPUT + ('Max. 50 frames are used, following frames ignored...',STAT) FRMCNT = 50 GOTO 1000 C C here we handle input from single input line C pull out operands (max. 50) C 500 BEGIN = 1 DO 800, N=1,50 CALL EXTRSS(LINE,',',BEGIN,FRAME(N),LL) IF (LL.LE.0) THEN FRMCNT = N - 1 GOTO 1000 ELSE CALL CLNFRA(FRAME(N),FRAME(N),0) ENDIF 800 CONTINUE FRMCNT = 50 C C get cursor area 1000 CALL STKRDC('INPUTC',1,1,80,IAV,LINE,UNI,NULO,STAT) IF (LINE(1:1).EQ.'C') THEN CALL STKRDR('OUTPUTR',10,10,IAV,COORDS,UNI,NULO,STAT) OSTART(1) = COORDS(3) OSTART(2) = COORDS(4) OEND(1) = COORDS(8) OEND(2) = COORDS(9) SUBDIM = 2 ENDIF C C now loop DO 3000, N=1,FRMCNT CALL STIGET(FRAME(N),D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIX,START,STEP, + IDENT,CUNIT,PNTRA,IMNOA,STAT) C C if not cursor input, we have to get the window coords now IF ((N.EQ.1) .AND. (LINE(1:1).EQ.'[')) THEN M = INDEX(LINE,',') IF (NAXIS.GT.1) THEN IAV = INDEX(LINE(M+1:),',') M = IAV + M ENDIF IF (M.LE.1) CALL STETER(34,'Invalid coordinate string...') LINE(M:M) = ':' CALL EXTCOO(IMNOA,LINE,NAXIS,SUBDIM,SUBLO,SUBHI,STAT) PIXEL(1,1) = SUBLO(1) PIXEL(2,1) = SUBLO(2) PIXEL(1,2) = SUBHI(1) PIXEL(2,2) = SUBHI(2) OSTART(1) = START(1) + (PIXEL(1,1)-1) * STEP(1) OEND(1) = START(1) + (PIXEL(1,2)-1) * STEP(1) IF (SUBDIM.GT.1) THEN OSTART(2) = START(2) + (PIXEL(2,1)-1) * STEP(2) OEND(2) = START(2) + (PIXEL(2,2)-1) * STEP(2) ENDIF ENDIF C C convert start + end of cursor region into pixel no.'s DIF = (OSTART(1)-START(1))/STEP(1) APIX(1,1) = INT(DIF) + 1 !BGVAL wants xstart, xend DIF = (OEND(1)-START(1))/STEP(1) !+ ystart, yend...! APIX(2,1) = INT(DIF) + 1 IF (SUBDIM.GT.1) THEN DIF = (OSTART(2)-START(2))/STEP(2) APIX(1,2) = INT(DIF) + 1 DIF = (OEND(2)-START(2))/STEP(2) APIX(2,2) = INT(DIF) + 1 ELSE APIX(1,2) = 1 APIX(2,2) = 1 ENDIF C C and get background stuff CALL BGVAL(KODE,MADRID(PNTRA),NPIX(1),NPIX(2),APIX,ITER, + MINP,LEAP,BGVL,RMSV,NPNTS,ISTP,ILAST) COORDS(1) = BGVL(1) !store mean + sigma continuously... COORDS(2) = RMSV(1) BAGR(N) = BGVL(1) IF (BAGR(N).GT.BMAX) BMAX = BAGR(N) !look for maximum... SIGM(N) = RMSV(1) IF (SIGM(N).LE.-1.0) THEN !check returned sigma OUTPUT(1:) = 'routine BGVAL failed on '//FRAME(N) CALL STETER(2,OUTPUT) ELSE LL = INDEX(FRAME(N),' ') IF (LL.LE.0) LL = LEN(FRAME(N)) WRITE(OUTPUT,10001) FRAME(N) CALL STTPUT(OUTPUT,STAT) WRITE(OUTPUT,10002) COORDS(1),COORDS(2) CALL STTPUT(OUTPUT,STAT) CALL STDWRR(IMNOA,'FLAT_BKG',COORDS,1,2,UNI,STAT) ENDIF C CALL STFCLO(IMNOA,STAT) !and release input frame 3000 CONTINUE C C now calculate weights and write them to descriptor WEIGHT DO 3300, N=1,FRMCNT W = BAGR(N)/BMAX CALL STFOPN(FRAME(N),D_OLD_FORMAT,0,F_IMA_TYPE,IMNOA,STAT) CALL STDWRR(IMNOA,'WEIGHT',W,1,1,UNI,STAT) 3300 CONTINUE C C That's it folks... RETURN C 10001 FORMAT('frame: ',A) 10002 FORMAT('mean, sigma of area = ',G12.6,', ',G12.6) END SUBROUTINE SUBRTG C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBRTG version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C artificial image C C.PURPOSE C rotate a 1-dim function around its start point to create a 2-dim image C C.ALGORITHM C rather straight forward C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/60 input frame C INPUTI/I/1/1 EVEN/ODD (1/0) option C OUT_A/C/1/60 output frame C NULL/R/1/1 user defined "null value" C C.VERSIONS C 1.00 from version 1.30 as of 861211 C move to FORTAN 77 + new St interfaces C C-------------------------------------------------- C IMPLICIT NONE C INTEGER NAXIS INTEGER*8 PNTRA,PNTRC INTEGER IMNOA,IMNOC INTEGER STAT,IAV,EVEN INTEGER NPIXA(3),NPIXC(3) INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*60 FRAMEA,FRAMEC CHARACTER CUNIT*48,IDENT*72 C REAL RNULL C DOUBLE PRECISION STEP(3),STARTA(3),STARTC(3) C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA IDENT /' '/, CUNIT /' '/ C C get input frame + map it CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,UNI,NULO,STAT) CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 3,NAXIS,NPIXA,STARTA,STEP, + IDENT,CUNIT,PNTRA,IMNOA,STAT) C C make sure, it's a 1-dim frame... IF (NAXIS.GT.1) + CALL STETER(1,'input should be 1-dim frame...') C C get EVEN or ODD option (variable EVEN = 1 or 0) CALL STKRDI('INPUTI',1,1,IAV,EVEN,UNI,NULO,STAT) C C get output frame + check, that it's different from input frame CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,UNI,NULO,STAT) IF (FRAMEA.EQ.FRAMEC) CALL STETER + (2,'input and result frame have to be different...') C C calculate new dimensions + map output file NPIXC(1) = (2 * NPIXA(1)) - 1 + EVEN NPIXC(2) = NPIXC(1) NPIXC(3) = 1 NAXIS = 2 C C start pixel of 1-dim frame is at world coord. = 0.0 STARTC(1) = STARTA(1) - ( (NPIXA(1) - 1 + EVEN) * STEP(1) ) STARTC(2) = STARTC(1) STARTC(3) = 0. C STEP(2) = STEP(1) STEP(3) = 1. C CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIXC,STARTC,STEP, + IDENT,CUNIT,PNTRC,IMNOC,STAT) C C now do it CALL STKRDR('NULL',1,1,IAV,RNULL,UNI,NULO,STAT) CALL DOROTA(MADRID(PNTRA),MADRID(PNTRC),NPIXA,NPIXC,EVEN,RNULL) C C update descriptors + add history CALL DSCUPT(IMNOA,IMNOC,' ',STAT) C RETURN END SUBROUTINE DOROTA(A,C,NPIXA,NPIXC,EVEN,RNULL) C IMPLICIT NONE C INTEGER NPIXA,NPIXC,EVEN INTEGER YOFF,YOFF2,NX,NY,I1,I2,NPX C REAL A(*),C(*),RNULL REAL DX,DY,DY2,DIST,F REAL SQRT C YOFF = 0 NPX = NPIXA - 1 C C *** C loop over first NPIXA - 1 lines C *** C DO 1000, NY=1,NPX DY = NY - NPIXA DY2 = DY * DY C C only calculate for NPIXA-1 values DO 400, NX=1,NPX DX = NX - NPIXA DIST = SQRT( (DX*DX) + DY2 ) + 1. !center pixel = 1 ... C C interpolate function value at coord. DIST I1 = INT(DIST) IF (I1.LT.NPIXA) THEN I2 = I1 + 1 F = A(I1) + ( (A(I2)-A(I1)) * (DIST-I1) ) ELSE IF (I1.EQ.NPIXA) THEN F = A(I1) ELSE F = RNULL ENDIF C(YOFF+NX) = F 400 CONTINUE C C points along center column ... DIST = ABS(NY-NPIXA) + 1. I1 = INT(DIST) I2 = I1 + 1 F = A(I1) + ( (A(I2)-A(I1)) * (DIST-I1) ) C(YOFF+NPIXA) = F C C if EVEN flag is set, repeat center column value IF (EVEN.EQ.1) C(YOFF+NPIXA+1) = F C C now the right half of each line YOFF2 = YOFF + NPIXA + EVEN DO 600, NX=1,NPIXA-1 C(YOFF2+NX) = C(YOFF2-NX-EVEN) 600 CONTINUE C YOFF = YOFF + NPIXC !update YOFF 1000 CONTINUE C C *** C build (upper) center line C *** C DO 1200, NX=1,NPIXA C(YOFF+NX) = A(NPIXA+1-NX) 1200 CONTINUE C C if EVEN flag is set, repeat center column value IF (EVEN.EQ.1) C(YOFF+NPIXA+1) = A(1) C C now the right half of this line YOFF2 = YOFF + NPIXA + EVEN DO 1400, NX=1,NPX C(YOFF2+NX) = C(YOFF2-NX-EVEN) 1400 CONTINUE C C if EVEN option - repeat center line IF (EVEN.EQ.1) THEN YOFF2 = NPIXA * NPIXC !point to line (NPIXA+1) YOFF = NPX * NPIXC !point to line NPIXA DO 1600, NX=1,NPIXC C(YOFF2+NX) = C(YOFF+NX) 1600 CONTINUE ENDIF C C *** C now repeat lower half of result frame C *** C YOFF2 = (NPIXA+EVEN) * NPIXC YOFF = (NPIXA-2) * NPIXC DO 2000, NY=1,NPX DO 1800, NX=1,NPIXC C(YOFF2+NX) = C(YOFF+NX) 1800 CONTINUE C YOFF2 = YOFF2 + NPIXC YOFF = YOFF - NPIXC 2000 CONTINUE C RETURN END SUBROUTINE SUBROT C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine SUBROT version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C pixels C C.PURPOSE C shift + rotate image around its start pixels by given angle C C.ALGORITHM C get the names of input + output frames from IN_A + OUT_B C get the new dimensions of output frame by rotating first the corner pixels C fill output frame with zeroes + rotate rest of input frame C C copy descriptor LHCUTS to output frame to keep same colours at same pixels! C C.INPUT/OUTPUT C the following keywords are used: C C IN_A/C/1/80 input frame C OUT_A/C/1/80 output frame C INPUTC rotation angle in degrees, C coords. of rotation point (may be = C), C scaling factors of stepsizes of input frame C or C R1,2,3 + P1,2,3 for ROTA/CLOCK + ROTA/COUNT C IN_B reference frame or "+" if no refframe C C.VERSIONS C 1.00 move to FORTRAN 77 + new interfaces C from version 2.40 as of 861029 C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,INTPOL,N,NAXISA INTEGER IMNOA,IMNOB,IMNOC INTEGER*8 PNTRA,PNTRC,PNTRD INTEGER NPIXA(3),NPIXB(3),STAT INTEGER N1,NVAL,STRIPE(2),MAPSIZ INTEGER NAXISC,ZOPT(2) INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*80 FRAMEA,FRAMEB,FRAMEC,LINE CHARACTER CUNITA*64,IDENTA*72,ACTION*2 C DOUBLE PRECISION STEPA(3),STARTA(3),STARTB(3),STEPB(3) DOUBLE PRECISION DVAL DOUBLE PRECISION DIN(4),DOUT(4) C REAL INPUTR(3),RNULL,RVAL REAL ANGDEG,ANGLE,AC,AS,AA,X1,X2,X3,Y1,Y2,Y3 REAL PIXROT(2),ROTOFF(2),ST(2),EN(2),RST(2),FST(2) REAL DIST,PIXSTP(2),QUOT,RQUOT,R1,R2,RROT(2) REAL CUTS(4),FACTO C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA FACTO /0.0174533/ ! Pi / 180. DATA INPUTR /1.0,1.0,0.0/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C ZOPT(1) = 0 ZOPT(2) = 0 INTPOL = 0 C C get input frame + map it CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NULO,STAT) CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT) CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXISA,UNI,NULO,STAT) CALL STDRDI(IMNOA,'NPIX',1,2,IAV,NPIXA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'START',1,2,IAV,STARTA,UNI,NULO,STAT) CALL STDRDD(IMNOA,'STEP',1,2,IAV,STEPA,UNI,NULO,STAT) CALL STDRDC(IMNOA,'CUNIT',1,1,64,IAV,CUNITA,UNI,NULO,STAT) CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENTA,UNI,NULO,STAT) IF (NAXISA.EQ.1) THEN CALL STETER(1,'input frame must be at least 2-dim frame...') ELSE IF (NAXISA.GT.2) THEN NAXISA = 2 !ignore higher dimensions... CALL STTPUT('only x-,y-plane will be processed...',STAT) ENDIF IF ((NPIXA(1).EQ.1) .OR. (NPIXA(2).EQ.1)) + CALL STETER(1,'We need a real 2-dim image ...') C C init frame pixls <-> world coords conversion CALL FPXWCO(0,IMNOA,DIN,DOUT,STAT) C C get name of output frame + User Null CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEB,UNI,NULO,STAT) CALL STKRDR('NULL',2,1,IAV,RNULL,UNI,NULO,STAT) C C check, that in, out_frame are different CALL GENEQF(FRAMEA,FRAMEB,IAV) IF (IAV.NE.0) + CALL STETER(9,'inframe and outframe are equal...') C C get string with rotation angle + rot_coords + scaling factors CALL STKRDC('INPUTC',1,1,60,IAV,LINE,UNI,NULO,STAT) C IF (LINE(1:1).EQ.'R') THEN CALL GENCNV(LINE(2:),1,1,NVAL,RVAL,DVAL,N) IF (NVAL.EQ.3) THEN !90 degrees counter clockwise INTPOL = 1 NPIXB(1) = NPIXA(2) NPIXB(2) = NPIXA(1) STEPB(1) = -STEPA(2) STEPB(2) = STEPA(1) DIN(1) = 1.0 DIN(2) = NPIXA(2) CALL FPXWCO(1,IMNOA,DIN,DOUT,STAT) STARTB(1) = DOUT(2) !STARTB(1) = ENDA(2) STARTB(2) = DOUT(1) !STARTB(2) = STARTA(1) C ELSE IF(NVAL.EQ.2) THEN !180 degrees INTPOL = 2 NPIXB(1) = NPIXA(1) NPIXB(2) = NPIXA(2) STEPB(1) = -STEPA(1) STEPB(2) = -STEPA(2) DIN(1) = NPIXA(1) DIN(2) = NPIXA(2) CALL FPXWCO(1,IMNOA,DIN,DOUT,STAT) STARTB(1) = DOUT(1) !STARTB(1) = ENDA(1) STARTB(2) = DOUT(2) !STARTB(2) = ENDA(2) ELSE !270 degrees INTPOL = 3 NPIXB(1) = NPIXA(2) NPIXB(2) = NPIXA(1) STEPB(1) = STEPA(2) STEPB(2) = -STEPA(1) DIN(1) = NPIXA(1) DIN(2) = 1 CALL FPXWCO(1,IMNOA,DIN,DOUT,STAT) STARTB(1) = DOUT(2) !STARTB(1) = STARTA(2) STARTB(2) = DOUT(1) !STARTB(2) = ENDA(1) ENDIF GOTO 3000 ENDIF C IAV = 1 CALL EXTRSS(LINE,',',IAV,FRAMEC,N) IF (N.LE.0) THEN CALL STETER(2,'Invalid rotation angle ...') ELSE CALL GENCNV(FRAMEC,2,1,NVAL,ANGDEG,DVAL,N) IF (N.LT.1) THEN CALL STETER(2,'Invalid rotation angle ...') ELSE ANGLE = ANGDEG * FACTO !convert from degrees to radians ENDIF ENDIF IF (LINE(IAV:IAV).EQ.'~') THEN N1 = (NPIXA(1) + 1) / 2 !use center pixel PIXROT(1) = N1 - 1 N1 = (NPIXA(2) + 1) / 2 !use center pixel PIXROT(2) = N1 - 1 GOTO 800 ENDIF C CALL EXTRSS(LINE,',',IAV,FRAMEC,N) IF (N.LE.0) THEN CALL STETER(3,'Invalid coords. of rotation point ...') ELSE PIXROT(1) = -9. IF ((FRAMEC(1:1).EQ.'C').OR.(FRAMEC(1:1).EQ.'c')) THEN N1 = (NPIXA(1) + 1) / 2 !use center pixel PIXROT(1) = N1 - 1 ELSE IF (FRAMEC(1:1).EQ.'<') THEN PIXROT(1) = 0. ELSE IF (FRAMEC(1:1).EQ.'>') THEN PIXROT(1) = NPIXA(1) - 1. ELSE IF (FRAMEC(1:1).EQ.'@') THEN CALL GENCNV(FRAMEC(2:),1,1,NVAL,RVAL,DVAL,N) PIXROT(1) = NVAL - 1. ELSE CALL GENCNV(FRAMEC,4,1,NVAL,RVAL,DVAL,N) PIXROT(1) = (DVAL - STARTA(1)) / STEPA(1) ENDIF ENDIF CALL EXTRSS(LINE,',',IAV,FRAMEC,N) IF (N.LE.0) THEN CALL STETER(3,'Invalid coords. of rotation point ...') ELSE PIXROT(2) = -9. IF ((FRAMEC(1:1).EQ.'C').OR.(FRAMEC(1:1).EQ.'c')) THEN N1 = (NPIXA(2) + 1) / 2 !use center pixel PIXROT(2) = N1 - 1 ELSE IF (FRAMEC(1:1).EQ.'<') THEN PIXROT(2) = 0. ELSE IF (FRAMEC(1:1).EQ.'>') THEN PIXROT(2) = NPIXA(2) - 1. ELSE IF (FRAMEC(1:1).EQ.'@') THEN CALL GENCNV(FRAMEC(2:),1,1,NVAL,RVAL,DVAL,N) PIXROT(2) = NVAL - 1. ELSE CALL GENCNV(FRAMEC,4,1,NVAL,RVAL,DVAL,N) PIXROT(2) = (DVAL - STARTA(2)) / STEPA(2) ENDIF ENDIF IF ((PIXROT(1).LT.-0.00001) .OR. + (PIXROT(2).LT.-0.00001)) + CALL STETER(3,'Invalid coords. of rotation point ...') C IF (LINE(IAV:IAV).NE.'~') THEN CALL GENCNV(LINE(IAV:),2,2,NVAL,INPUTR,DVAL,N) IF (N.LT.2) CALL STETER(4,'Invalid scale values ...') ENDIF C 800 AC = COS(ANGLE) AS = SIN(ANGLE) AA = MOD(ANGDEG,90.) !rotation by N*90 degrees ? IF (ABS(AA).LT.0.001) THEN IF (AS.LT.-0.1) THEN !yes. no interpolation needed... INTPOL = 3 !270 degree NPIXB(1) = NPIXA(2) NPIXB(2) = NPIXA(1) ELSE IF (AS.LT.0.1) THEN IF (AC.LT.-0.1) THEN INTPOL = 2 !180 degrees ELSE INTPOL = 4 !360 degrees ENDIF NPIXB(1) = NPIXA(1) NPIXB(2) = NPIXA(2) ELSE INTPOL = 1 !90 degrees NPIXB(1) = NPIXA(2) NPIXB(2) = NPIXA(1) ENDIF ENDIF ENDIF C C to get size of output frame, rotate corners around lower left corner C (assumed = 0,0) C R1 = NPIXA(1) - 1. R2 = NPIXA(2) - 1. CALL ROTO(R1,R2,AC,AS,X1,Y1) !(NX-1,NY-1) CALL ROTO(0.,R2,AC,AS,X2,Y2) !(0,NY-1) CALL ROTO(R1,0.,AC,AS,X3,Y3) !(NX-1,0) C C get result frame "corners" ST(1) = MIN(0.,X1,X2,X3) !rotated x-start EN(1) = MAX(0.,X1,X2,X3) !rotated x-end ST(2) = MIN(0.,Y1,Y2,Y3) !rotated y-start EN(2) = MAX(0.,Y1,Y2,Y3) !rotated y-end C C check refframe option CALL STKRDC('IN_B',1,1,80,IAV,FRAMEC,UNI,NULO,STAT) IF (FRAMEC(1:1).NE.'+') THEN ACTION(2:2) = 'F' CALL STFOPN(FRAMEC,D_OLD_FORMAT,0,F_IMA_TYPE,IMNOC,STAT) CALL STDRDI(IMNOC,'NAXIS',1,1,IAV,NAXISC,UNI,NULO,STAT) CALL STDRDI(IMNOC,'NPIX',1,3,IAV,NPIXB,UNI,NULO,STAT) IF (NAXISC.GT.2) NAXISC = 2 CALL STDRDD(IMNOC,'START',1,2,IAV,STARTB,UNI,NULO,STAT) CALL STDRDD(IMNOC,'STEP',1,2,IAV,STEPB,UNI,NULO,STAT) ELSE ACTION(2:2) = 'Q' ENDIF C C from coords. of rotation point determine shift after rotation QUOT = 1. - AC ROTOFF(1) = PIXROT(1)*QUOT + PIXROT(2)*AS ROTOFF(2) = PIXROT(2)*QUOT - PIXROT(1)*AS C C choose sample start in 0,0 pixel space according to option C IF (ACTION(2:2).EQ.'F') THEN DO 1000, N=1,2 PIXSTP(N) = STEPB(N) / STEPA(N) FST(N) = (STARTB(N) - STARTA(N)) / STEPA(N) ST(N) = FST(N) - ROTOFF(N) EN(N) = ST(N) + ( (NPIXB(N) - 1) * PIXSTP(N) ) 1000 CONTINUE ELSE DO 1500, N=1,2 PIXSTP(N) = INPUTR(N) STEPB(N) = STEPA(N) * PIXSTP(N) C C test, if we want to rotate around a "sample" point or not QUOT = PIXROT(N) / PIXSTP(N) RQUOT = FLOAT(INT(QUOT)) DIST = ABS(PIXROT(N) - RQUOT*PIXSTP(N)) IF (DIST.LE.(0.001*PIXSTP(N))) ZOPT(N) = 1 1500 CONTINUE C C choose corners such, that rotation point is hit ... N1 = ZOPT(1) * ZOPT(2) IF ( (INTPOL.EQ.0) .AND. (N1.EQ.1) ) THEN CALL ROTO(PIXROT(1),PIXROT(2),AC,AS,RROT(1),RROT(2)) DO 2000, N=1,2 DIST = !distance of fixed point to rotation point + RROT(N) - ST(N) QUOT = DIST / PIXSTP(N) !is forced to RQUOT = FLOAT(INT(QUOT)) + 1. !multiples of steps DIST = RQUOT * PIXSTP(N) ST(N) = RROT(N) - DIST DIST = EN(N) - ST(N) !distance of start to end point QUOT = DIST / PIXSTP(N) !is forced to RQUOT = FLOAT(INT(QUOT)) + 1. !multiples of steps EN(N) = ST(N) + (RQUOT*PIXSTP(N)) NPIXB(N) = NINT(RQUOT) + 1 !no. of points in new image 2000 CONTINUE ELSE NPIXB(1) = NINT( (EN(1) - ST(1))/PIXSTP(1) ) + 1 NPIXB(2) = NINT( (EN(2) - ST(2))/PIXSTP(2) ) + 1 ENDIF C C now update start of rotated image, so that rotation point stays fixed RST(1) = ST(1) + ROTOFF(1) RST(2) = ST(2) + ROTOFF(2) STARTB(1) = STARTA(1) + ( RST(1)*STEPA(1) ) !we work in 0,0 space STARTB(2) = STARTA(2) + ( RST(2)*STEPA(2) ) ENDIF C C and map it... 3000 N = NPIXB(1)*NPIXB(2) CALL STFCRE(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,N, + IMNOB,STAT) CALL STDWRI(IMNOB,'NAXIS',NAXISA,1,1,UNI,STAT) CALL STDWRI(IMNOB,'NPIX',NPIXB,1,NAXISA,UNI,STAT) CALL STDWRD(IMNOB,'START',STARTB,1,NAXISA,UNI,STAT) CALL STDWRD(IMNOB,'STEP',STEPB,1,NAXISA,UNI,STAT) CALL STDWRC(IMNOB,'IDENT',1,IDENTA,1,72,UNI,STAT) N = (NAXISA+1) * 16 CALL STDWRC(IMNOB,'CUNIT',1,CUNITA,1,N,UNI,STAT) C CALL STKRDI('MONITPAR',20,1,IAV,MAPSIZ,UNI,NULO,STAT) MAPSIZ = MAPSIZ * MAPSIZ !square frame ... IF (INTPOL.EQ.0) THEN N = NPIXA(1) * NPIXA(2) CALL STFMAP(IMNOA,F_I_MODE,1,N,IAV,PNTRA,STAT) C C for output buffer aim at size of the `mapsize frame' STRIPE(2) = MAPSIZ / NPIXB(1) IF (NPIXB(2).LT.STRIPE(2)) STRIPE(2) = NPIXB(2) N = NPIXB(1)*STRIPE(2) !size of output strip CALL STFXMP(N,D_R4_FORMAT,PNTRD,STAT) C INPUTR(1) = AC INPUTR(2) = AS INPUTR(3) = RNULL CALL ROTAIT(MADRID(PNTRA),MADRID(PNTRD),IMNOB,STRIPE(2), + NPIXA,NPIXB,ST,PIXSTP,INPUTR) ELSE !test, if direct disk I/O C C for each buffer aim at size of the `mapsize frame' STRIPE(1) = MAPSIZ / NPIXA(1) IF (NPIXA(2).LT.STRIPE(1)) STRIPE(1) = NPIXA(2) N = NPIXA(1)*STRIPE(1) !size of input strip CALL STFXMP(N,D_R4_FORMAT,PNTRC,STAT) STRIPE(2) = MAPSIZ / NPIXA(2) IF (NPIXA(1).LT.STRIPE(2)) STRIPE(2) = NPIXA(1) N = NPIXA(2)*STRIPE(2) !size of output strip CALL STFXMP(N,D_R4_FORMAT,PNTRD,STAT) C CALL ROTM90(INTPOL,MADRID(PNTRC),MADRID(PNTRD), + IMNOA,IMNOB,NPIXA,STRIPE) ENDIF C C copy descriptors + LHCUTS CALL DSCUPT(IMNOA,IMNOB,' ',STAT) CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT) CALL STDWRR(IMNOB,'LHCUTS',CUTS,1,4,UNI,STAT) C RETURN C END SUBROUTINE ROTO(X,Y,AC,AS,XOUT,YOUT) C IMPLICIT NONE C REAL X,Y,AC,AS,XOUT,YOUT C XOUT = X*AC - Y*AS YOUT = X*AS + Y*AC C RETURN END SUBROUTINE ROTAIT(A,B,IMNOB,STRIPB,NPIXA,NPIXB, + START,STEP,WORK) C C K. Banse 861029, 880902, 910204, 920306 C C do "backward" rotation, i.e. solve the system XROT = X*AC - Y*AS C YROT = X*AS + Y*AC for X,Y C to yield X = XROT*AC + YROT*AS C Y = YROT*AC - XROT*AS which takes AC**2 + AS**2 = 1 into account.... C IMPLICIT NONE C INTEGER IMNOB,STRIPB INTEGER NPIXA(2),NPIXB(2) INTEGER INDXB,KSTX,KX,KY INTEGER INOFF,DATOUT,KSIZE,YLIN INTEGER NO1,NO2,NO3,NO4,STAT INTEGER NX,NY,SIZEA,XDIMA,YDIMA C REAL A(*),B(*) REAL WORK(3),AC,AS,ZNULL REAL XT,YT,XFRAC,YFRAC,XLAST,YLAST REAL START(2),STEP(2) REAL XF,XSTEP,YSTEP,RX,RY,RYAS,RYAC C C init AC = WORK(1) AS = WORK(2) ZNULL = WORK(3) XF = START(1) XSTEP = STEP(1) YSTEP = STEP(2) XDIMA = NPIXA(1) YDIMA = NPIXA(2) XLAST = XDIMA - 1. YLAST = YDIMA - 1. SIZEA = XDIMA * YDIMA DATOUT = 1 YLIN = STRIPB INOFF = 0 RY = START(2) !init RY to start-y C C 400 IF ((INOFF+YLIN).GT.NPIXB(2)) YLIN = NPIXB(2) - INOFF INDXB = 0 C C loop through result lines DO 3000, NY=1,YLIN C RX = XF !init RX to start-x RYAS = RY * AS RYAC = RY * AC C DO 1000, NX=1,NPIXB(1) INDXB = INDXB + 1 !loop along result line XT = RX*AC + RYAS IF ((XT.LT.0.).OR.(XT.GT.XLAST)) THEN B(INDXB) = ZNULL GOTO 555 ENDIF YT = RYAC - RX*AS IF ((YT.LT.0.).OR.(YT.GT.YLAST)) THEN B(INDXB) = ZNULL GOTO 555 ENDIF C C interpolate... XT = XT + 1. !adjust 0,0 to 1,1 ... YT = YT + 1. KX = INT(XT) KY = INT(YT) XFRAC = XT - KX YFRAC = YT - KY KSTX = XDIMA*(KY-1) !get 1. pixel in this line - 1 NO1 = KX + KSTX !get pixel close to integer part of real pixel C NO2 = NO1 + 1 !get next pixel in x-direction IF (NO2-KSTX.GT.XDIMA) THEN !are we out of grid in x? IF (NO2.GT.SIZEA) THEN !also out in y? B(INDXB) = A(NO1) !yes. ELSE NO3 = NO1 + XDIMA !no. C> pixel = A(NO1)*(1.-YFRAC) + A(NO3)*YFRAC pixel = A(NO1)*(1.-XFRAC) + A(NO2)*YFRAC pixel = ( A(NO1)*(1.-XFRAC) + A(NO2)*XFRAC ) * (1.-YFRAC) + + ( A(NO3)*(1.-XFRAC) + A(NO4)*XFRAC ) * YFRAC