C @(#)align.for 17.1.1.1 (ESO-IPG) 01/25/02 17:39: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 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 ALIGN C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C program ALIGN version 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS: C images, rotation C C.PURPOSE: C given 2 tables with coordinates of the same objects in different images, C calculate the coefficients of the transformation function C C.ALGORITHM: C Use least squares approximation. C C.INPUT/OUTPUT: C C The following keywords are used: C ACTION/C/1/1 A for aligning, R for rotation C IN_A/C/1/80 input_table1 C IN_B/C/1/80 input_table2 C INPUTC/C/1/1 method for transformation C P4/C/1/20 overlay flag, optional intensity C C TRANSFRM/D/1/9 this key will be filled with the C transformation values as in IHAP C C.VERSIONS C 1.00 taken from version 3.10 as of 871123 of original ALIGN C which was written by P. Grosbol C C 010110 last modif C C-------------------------------------------------------------------------- C IMPLICIT NONE C INTEGER COLNM1(3),COLNM2(3),COLNMR(3) INTEGER BOX(4),XFIG(5),YFIG(5) INTEGER OVINTS INTEGER IAV,ICOUNT INTEGER N,NN,NOPNTS INTEGER NCNT1,NCNT2,NROW1,NROW2 INTEGER STAT,LL,IDNO(2000) INTEGER TID1,TID2,IMNO INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*80 TABLE1,TABLE2,FRAME CHARACTER CBUF*80,METHOD*1,RESIDU*1,OVFLAG*20 CHARACTER LABEL1(3)*16,LABEL2(3)*16,IDENT1*8 CHARACTER IDENT2*8,ID(2000)*8 C REAL R1,R2 REAL XX1(2,2000),XX2(2,2000),TEMP1(3),TEMP2(3) REAL A,A1,A2,A3,A4,CX,CY,R,X,Y REAL RBUF(17),PBUF(6) REAL RINF(8) C LOGICAL TABNUL(3) LOGICAL ISEL1, ISEL2 C DOUBLE PRECISION TR(6),F1,F2,ANG,DBUF(9) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C INCLUDE 'MID_INCLUDE:IDIMEM.INC' INCLUDE 'MID_INCLUDE:IDIDEV.INC' C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up MIDAS environment + enable automatic error abort C CALL STSPRO('ALIGN') ICOUNT = 0 CALL STKRDC('ACTION',1,1,1,IAV,METHOD,UNI,NULO,STAT) IF (METHOD.EQ.'R') THEN CALL ALIGNR CALL STSEPI ENDIF C C get input tables and columns CALL STKRDC('IN_A',1,1,80,IAV,TABLE1,UNI,NULO,STAT) CALL STKRDC('IN_B',1,1,80,IAV,TABLE2,UNI,NULO,STAT) C C look for column labels LABEL1(3) = 'IDENT ' LABEL2(3) = 'IDENT ' LL = INDEX(TABLE1,',') IF (LL.LE.0) THEN LABEL1(1) = 'XCEN ' LABEL1(2) = 'YCEN ' ELSE CBUF(1:) = TABLE1(LL+1:)//' ' TABLE1(LL:) = ' ' LL = INDEX(CBUF,',') IF (LL.LE.0) + CALL STETER(9,'invalid syntax for 1. input table') LABEL1(1) = CBUF(2:LL-1) !omit the leading ':' LABEL1(2) = CBUF(LL+2:) ENDIF C LL = INDEX(TABLE2,',') IF (LL.LE.0) THEN LABEL2(1) = 'XCEN ' LABEL2(2) = 'YCEN ' ELSE CBUF(1:) = TABLE2(LL+1:)//' ' TABLE2(LL:) = ' ' LL = INDEX(CBUF,',') IF (LL.LE.0) + CALL STETER(9,'invalid syntax for 2. input table') LABEL2(1) = CBUF(2:LL-1) !omit the leading ':' LABEL2(2) = CBUF(LL+2:) ENDIF C IF (TABLE1.EQ.TABLE2) + CALL STETER(1,'input and output table have to be different...') C get method CALL STKRDC('INPUTC',1,1,2,IAV,CBUF,UNI,NULO,STAT) CALL UPCAS(CBUF,CBUF) !make sure we have upper case METHOD = CBUF(1:1) RESIDU = CBUF(2:2) C C test, if required columns are there CALL TBTOPN(TABLE1,F_IO_MODE,TID1,STAT) CALL TBTOPN(TABLE2,F_I_MODE,TID2,STAT) DO 100,N=1,3 CALL TBLSER(TID1,LABEL1(N),COLNM1(N),STAT) IF (COLNM1(N).EQ.-1) THEN CBUF(1:) = 'column '//LABEL1(N)// + 'not found in table '//TABLE1//' ' CALL STETER(2,CBUF) ENDIF CALL TBLSER(TID2,LABEL2(N),COLNM2(N),STAT) IF (COLNM2(N).EQ.-1) THEN CBUF(1:) = 'column '//LABEL2(N)// + 'not found in table '//TABLE2//' ' CALL STETER(2,CBUF) ENDIF 100 CONTINUE C C check, that IDENT columns are of type character CALL TBFGET(TID1,COLNM1(3),CBUF,IAV,NN,STAT) IF (NN.NE.D_C_FORMAT) THEN CALL STETER(22, + '1. Table: Column :IDENT not character type ...') ENDIF CALL TBFGET(TID2,COLNM2(3),CBUF,IAV,NN,STAT) IF (NN.NE.D_C_FORMAT) THEN CALL STETER(22, + '2. Table: Column :IDENT not character type ...') ENDIF C C if RESIDU is set, look for residual columns IF (RESIDU.EQ.'Y') THEN CBUF(1:) = ' ' CALL TBLSER(TID1,'XRESIDUAL',COLNMR(1),STAT) IF (COLNMR(1).EQ.-1) THEN CALL TBCINI(TID1,D_R4_FORMAT,1,'G12.6',CBUF(1:16), + 'XRESIDUAL',COLNMR(1),STAT) ENDIF CALL TBLSER(TID1,'YRESIDUAL',COLNMR(2),STAT) IF (COLNMR(2).EQ.-1) THEN CALL TBCINI(TID1,D_R4_FORMAT,1,'G12.6',CBUF(1:16), + 'YRESIDUAL',COLNMR(2),STAT) ENDIF CALL TBLSER(TID1,'RESIDUALS',COLNMR(3),STAT) IF (COLNMR(3).EQ.-1) THEN CALL TBCINI(TID1,D_R4_FORMAT,1,'G12.6',CBUF(1:16), + 'RESIDUALS',COLNMR(3),STAT) ENDIF ENDIF C C get values from tables CALL TBIGET(TID1,N,NROW1,N,N,N,STAT) !get total no. of rows CALL TBIGET(TID2,N,NROW2,N,N,N,STAT) !get total no. of rows C C here we loop and look for equal identifiers DO 1000,NCNT1=1,NROW1 CALL TBSGET(TID1,NCNT1,ISEL1,STAT) IF (ISEL1) THEN CALL TBRRDR(TID1,NCNT1,2,COLNM1,TEMP1, !get next row of table1 + TABNUL,STAT) CALL TBERDC(TID1,NCNT1,COLNM1(3), + IDENT1,TABNUL,STAT) DO 400,NCNT2=1,NROW2 CALL TBSGET(TID2,NCNT2,ISEL2,STAT) IF (ISEL2) THEN CALL TBERDC(TID2,NCNT2, + COLNM2(3),IDENT2,TABNUL,STAT) IF (IDENT1.EQ.IDENT2) THEN !match, read also x,y values CALL TBRRDR(TID2,NCNT2,2,COLNM2,TEMP2, + TABNUL,STAT) GOTO 800 ENDIF ENDIF 400 CONTINUE ENDIF C C no matching identifier in TABLE2 so skip this row GOTO 1000 C C matching identifiers 800 ICOUNT = ICOUNT + 1 IF (ICOUNT.GT.2000) THEN ICOUNT = 2000 GOTO 1500 !we only take max. 2000 points ELSE ID(ICOUNT) = IDENT1 IDNO(ICOUNT) = NCNT1 !save row no. XX1(1,ICOUNT) = TEMP1(1) XX1(2,ICOUNT) = TEMP1(2) XX2(1,ICOUNT) = TEMP2(1) XX2(2,ICOUNT) = TEMP2(2) ENDIF 1000 CONTINUE C C check if errors C 1500 IF (ICOUNT.EQ.0) + CALL STETER(3,'no identifiers match') C C now calculate the transformation C CALL TRANS(XX2,XX1,ICOUNT,METHOD,ANG,F1,F2,TR,STAT) IF (STAT.NE.0) + CALL STETER(4,'problems in calculating the transformation') C C fill keyword TRANSFRM C ANG = 57.29578D0 * ANG !transform to degrees DBUF(1) = ANG DBUF(2) = F1 DBUF(3) = F2 DO 1800,N=1,6 DBUF(N+3) = TR(N) 1800 CONTINUE CALL STKWRD('TRANSFRM',DBUF,1,9,UNI,STAT) C C display results WRITE(CBUF,10001) DBUF(1) CALL STTPUT(CBUF,STAT) WRITE(CBUF,10002) DBUF(2),DBUF(3) CALL STTPUT(CBUF,STAT) CBUF = ' ' CALL STTPUT(CBUF,STAT) C C if desired, draw computed new position into overlay plane CALL STKRDC('P4',1,1,20,IAV,OVFLAG,UNI,NULO,STAT) IF (OVFLAG(1:1).EQ.'o') OVFLAG(1:1) = 'O' IF (OVFLAG(1:1).EQ.'O') THEN CALL DTOPEN(1,STAT) NN = INDEX(OVFLAG,',') IF (NN.LE.0) THEN OVINTS = 255 !default intensity is 255 ELSE CALL GENCNV(OVFLAG(NN+1:),1,1,OVINTS,OVINTS,OVINTS,N) IF (N.LE.0) OVINTS = 255 ENDIF C C get info of displayed image CALL DTGICH(QDSPNO,QIMCH,FRAME,RINF,STAT) CALL STFOPN(FRAME,D_OLD_FORMAT,0,F_IMA_TYPE,IMNO,STAT) CALL PIXXCV('INIT',IMNO,RBUF,STAT) !set up PIXXCV ENDIF C C compute residuals and mean error A = 0. A1 = TR(1) A2 = TR(2) A3 = TR(3) A4 = TR(4) CX = TR(5) CY = TR(6) C WRITE(CBUF,10003) CALL STTPUT(CBUF,STAT) DO 2200,N=1,ICOUNT X = A1*XX1(1,N) + A2*XX1(2,N) + CX Y = A3*XX1(1,N) + A4*XX1(2,N) + CY R1 = X - XX2(1,N) R2 = Y - XX2(2,N) R = SQRT( R1*R1 + R2*R2 ) WRITE(CBUF,10004) ID(N),X,Y,R1,R2,R IF (RESIDU.EQ.'Y') THEN !save residuals in input table TEMP1(1) = R1 TEMP1(2) = R2 TEMP1(3) = R CALL TBRWRR(TID1,IDNO(N),3,COLNMR,TEMP1,STAT) ENDIF C CALL STTPUT(CBUF,STAT) A = A + R !sum up... IF (OVFLAG(1:1).EQ.'O') THEN PBUF(1) = X PBUF(2) = Y CALL PIXXCV('WRS',0,PBUF,STAT) !world c. -> screen c. BOX(1) = PBUF(5) - 5 BOX(2) = PBUF(6) - 5 BOX(3) = PBUF(5) + 5 BOX(4) = PBUF(6) + 5 TEMP1(1) = -1.0 TEMP1(2) = -1.0 CALL BLDGRA('REC',BOX,TEMP1,XFIG,YFIG,5,NOPNTS) CALL IIGPLY(QDSPNO,QOVCH,XFIG,YFIG,NOPNTS,OVINTS,1,STAT) ENDIF 2200 CONTINUE C A = A/ICOUNT !normalize... WRITE(CBUF,10000) A CALL STTPUT(CBUF,STAT) IF (OVFLAG(1:1).EQ.'O') CALL DTCLOS(QDSPNO) C C that's it folks... C CALL STSEPI C C Formats 10000 FORMAT('mean error = ',E15.6) 10001 FORMAT('rotation angle = ',F8.3,' degrees') 10002 FORMAT('scaling factors in x,y :',E15.6,', ',E15.6) 10003 FORMAT('ID X Y x-residual ', + ' y-residual residuals') 10004 FORMAT(A,5E14.6) C END SUBROUTINE TRANS(X1,X2,N,METH,ANG,F1,F2,TR,STAT) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C subroutine TRANS version 2.10 850812 C K. Banse ESO - Garching C C.KEYWORDS: C transformation, coordinate sets C C.PURPOSE: C Calculate the coefficients of a linear transformation between the C two point sets X1 and X2 of N points each. C The transformation is assumed to be of the form: C C X1(1,k) = TR(1)*X2(1,k) + TR(2)*X2(2,k) + TR(5) C X1(2,k) = TR(3)*X2(1,k) + TR(4)*X2(2,k) + TR(6) C C or alternatively C C X1(1,k) = F1*cos(ANG)*X2(1,k) + F2*sin(ANG)*X2(2,k) + TR(5) C X1(2,k) = -F1*sin(ANG)*X2(1,k) + F2*cos(ANG)*X2(2,k) + TR(6) C C.ALGORITHM: C Use least square algorithm. C C.INPUT/OUTPUT: C call as TRANS(X1,X2,N,METH,ANG,F1,F2,TR,STAT) C C input par: C X1: R*4 array 1. set of points C X2: R*4 array 2. set of points C N: I*4 no. of points in X1 and X2 C METH: char.exp. method to use, C F = free paramaters, E = F1 and F2 are equal C U = F1 = F2 = 1.0 C S = shift only: F1 = F2 = 1.0, angle = 0.0 C C output par: C ANG: R*8 rotation angle in radians C F1: R*8 scaling factor in x-direction C F2: R*8 scaling factor in y-direction C TR: R*8 array array of coefficients C ANG, F1, F2 are derived from TR ... C STAT: I*4 return status C C-------------------------------------------------------------------------- C IMPLICIT NONE C CHARACTER*(*) METH C INTEGER I,N,STAT C REAL X1(2,N),X2(2,N) C DOUBLE PRECISION TR(6),F1,F2,ANG,AC DOUBLE PRECISION RN,X,Y,XI,YI,XX,XY,YY,XXI,XYI,XIY,YYI DOUBLE PRECISION A,B,C,D,DINV DOUBLE PRECISION DATAN2,DCOS,DSIGN,DSIN,DSQRT !system functions... C C check input IF ( ((METH(1:1).EQ.'F').AND.(N.LT.3)) .OR. + ((METH(1:1).EQ.'E').AND.(N.LT.2)) .OR. + ((METH(1:1).EQ.'U').AND.(N.LT.2)) .OR. + ((METH(1:1).EQ.'S').AND.(N.LT.1)) ) THEN CALL STTPUT('not enough points in tables...',STAT) STAT = -1 RETURN ELSE STAT= 0 ENDIF C C init sums RN = 1.D0 / DBLE (FLOAT(N) ) X = 0.D0 Y = 0.D0 XI = 0.D0 YI = 0.D0 XX = 0.D0 XY = 0.D0 YY = 0.D0 XXI = 0.D0 XYI = 0.D0 XIY = 0.D0 YYI = 0.D0 C C compute the sums DO 1000, I=1,N A = DBLE (X2(1,I)) B = DBLE (X2(2,I)) C = DBLE (X1(1,I)) D = DBLE (X1(2,I)) C X = X + A !sum of x-coords of X2 Y = Y + B !sum of y-coords of X2 XX = XX + A*A XY = XY + A*B YY = YY + B*B XI = XI + C YI = YI + D XXI = XXI + A*C XYI = XYI + A*D XIY = XIY + C*B YYI = YYI + B*D 1000 CONTINUE C C normalize coefficients XX = XX - X*X*RN XY = XY - X*Y*RN YY = YY - Y*Y*RN XXI = XXI - X*XI*RN XYI = XYI - X*YI*RN XIY = XIY - XI*Y*RN YYI = YYI - Y*YI*RN C C branch according to method IF (METH(1:1).EQ.'E') GOTO 2000 IF (METH(1:1).EQ.'U') GOTO 3000 IF (METH(1:1).EQ.'S') GOTO 4000 C C fall through to method = FREE (all paramters are free) D = XX*YY - XY*XY IF (DABS(D).LT.1.D-20) THEN CALL STTPUT('points not well chosen...',STAT) STAT = 1 RETURN ENDIF C A = (YY*XXI - XY*XIY) / D B = (XX*XIY - XY*XXI) / D C = (YY*XYI - XY*YYI) / D D = (XX*YYI - XY*XYI) / D C TR(1) = A TR(2) = B TR(3) = C TR(4) = D TR(5) = (XI - A*X - B*Y)*RN TR(6) = (YI - C*X - D*Y)*RN C ANG = DATAN2(B,D) F1 = DSQRT(A*A + C*C) F2 = DSQRT(B*B + D*D) AC = DCOS(ANG) IF (DABS(AC).GT.0.5D0) THEN F1 = DSIGN( F1,TR(1)*AC) F2 = DSIGN( F2,TR(4)*AC) ELSE AC = DSIN(ANG) F1 = DSIGN( F1,-TR(3)*AC) F2 = DSIGN( F2,TR(2)*AC) ENDIF RETURN C C method = EQUAL (F1 = F2) 2000 D = XX + YY IF (DABS(D).LT.1.D-20) THEN CALL STTPUT('points not well chosen...',STAT) STAT = 1 RETURN ELSE DINV = 1.D0 / D ENDIF C A = (XXI + YYI) * DINV B = (XIY - XYI) * DINV C TR(1) = A TR(2) = B TR(3) = -TR(2) TR(4) = TR(1) TR(5) = (XI - A*X - B*Y)*RN TR(6) = (YI + B*X - A*Y)*RN C ANG = DATAN2(B,A) F1 = DSQRT(A*A+B*B) F2 = F1 RETURN C C method = UNIT (F1 = F2 = 1.0) 3000 A = DATAN2( (XIY-XYI),(XXI+YYI) ) ANG = A TR(1) = DCOS(ANG) TR(2) = DSIN(ANG) TR(3) = -TR(2) TR(4) = TR(1) TR(5) = (XI - TR(1)*X - TR(2)*Y)*RN TR(6) = (YI + TR(2)*X - TR(1)*Y)*RN C F1 = 1.0D0 F2 = 1.0D0 RETURN C C method = SHIFT (F1 = F2 = 1.0, angle = 0.0) 4000 ANG = 0.0D0 TR(1) = 1.0D0 TR(2) = 0.0D0 TR(3) = 0.0D0 TR(4) = 1.0D0 TR(5) = (XI - X) * RN TR(6) = (YI - Y) * RN C F1 = 1.0D0 F2 = 1.0D0 RETURN C END SUBROUTINE ALIGNR C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine ALIGNR version 1.00 861028 C K. Banse ESO - Garching C 1.10 910205 C C.KEYWORDS C pixels C C.PURPOSE C scale + rotate image according to output from ALIGN/IMAGE 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 TRANSFRM/D/1/9 rotation specs + transfromation matrix C IN_B/C/1/80 reference frame, if not "+" C NULL?R/2/1 User Nullvalue C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,N,NAXISA INTEGER*8 PNTRA,PNTRB INTEGER IMNOA,IMNOB,IMNOR,STAT INTEGER NPIXA(3),NPIXB(3) INTEGER NAXISC INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*80 FRAMEA,FRAMEB,FRAMEC CHARACTER CUNITA*64,IDENTA*72,ACTION*2 C REAL RNULL C DOUBLE PRECISION STEPA(3),STARTA(3) DOUBLE PRECISION STARTB(3),STEPB(3),STARTC(3),END(2) DOUBLE PRECISION X1,X2,X3,X4,Y1,Y2,Y3,Y4 DOUBLE PRECISION ST(2),EN(2),F1,F2,A1,A2,A3,A4 DOUBLE PRECISION CX,CY,TR(6),CUTS(4) DOUBLE PRECISION TOT(9),TRR(6) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA IDENTA /' '/ DATA CUNITA /' '/ C C get input frame + map it C 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) IF (NAXISA.EQ.1) THEN CALL STETER(1,'We need a 2-dim image ...') ELSE IF (NAXISA.GT.2) THEN NAXISA = 2 CALL STTPUT('We only work on the 1. plane ...',STAT) ENDIF IF ((NPIXA(1).EQ.1).OR.(NPIXA(2).EQ.1)) + CALL STETER(1,'We need a 2-dim image ...') C C get output frame CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEB,UNI,NULO,STAT) C C get transfomation coefficients + user Null CALL STKRDD('TRANSFRM',1,9,IAV,TOT,UNI,NULO,STAT) F1 = TOT(2) F2 = TOT(3) A1 = TOT(4) A2 = TOT(5) A3 = TOT(6) A4 = TOT(7) CX = TOT(8) CY = TOT(9) CALL STKRDR('NULL',2,1,IAV,RNULL,UNI,NULO,STAT) C C to get size of output frame, rotate corners END(1) = STARTA(1) + (NPIXA(1)-1)*STEPA(1) END(2) = STARTA(2) + (NPIXA(2)-1)*STEPA(2) C X1 = A1*STARTA(1) + A2*STARTA(2) + CX Y1 = A3*STARTA(1) + A4*STARTA(2) + CY X2 = A1*STARTA(1) + A2*END(2) + CX Y2 = A3*STARTA(1) + A4*END(2) + CY X3 = A1*END(1) + A2*STARTA(2) + CX Y3 = A3*END(1) + A4*STARTA(2) + CY X4 = A1*END(1) + A2*END(2) + CX Y4 = A3*END(1) + A4*END(2) + CY C C get result frame "corners" ST(1) = MIN(X1,X2,X3,X4) !rotated x-start EN(1) = MAX(X1,X2,X3,X4) !rotated x-end ST(2) = MIN(Y1,Y2,Y3,Y4) !rotated y-start EN(2) = MAX(Y1,Y2,Y3,Y4) !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,1,IMNOR,STAT) CALL STDRDI(IMNOR,'NAXIS',1,1,IAV,NAXISC,UNI,NULO,STAT) CALL STDRDI(IMNOR,'NPIX',1,3,IAV,NPIXB,UNI,NULO,STAT) IF (NAXISC.EQ.1) THEN CALL STETER(10,'We need a 2-dim reference image ...') ELSE IF (NAXISC.GT.2) THEN NAXISC = 2 ENDIF IF ((NPIXB(1).EQ.1).OR.(NPIXB(2).EQ.1)) + CALL STETER(10,'We need a 2-dim reference image ...') CALL STDRDD(IMNOR,'START',1,2,IAV,STARTC,UNI,NULO,STAT) CALL STDRDD(IMNOR,'STEP',1,2,IAV,STEPB,UNI,NULO,STAT) ELSE ACTION(2:2) = 'Q' STEPB(1) = STEPA(1) * F1 STEPB(2) = STEPA(2) * F2 ENDIF C IF (STEPB(1).LT.0.) THEN STARTB(1) = EN(1) ELSE STARTB(1) = ST(1) ENDIF IF (STEPB(2).LT.0.) THEN STARTB(2) = EN(2) ELSE STARTB(2) = ST(2) ENDIF C NPIXB(1) = NINT( (EN(1) - ST(1))/ABS(STEPB(1)) ) + 1 NPIXB(2) = NINT( (EN(2) - ST(2))/ABS(STEPB(2)) ) + 1 C C I think the following code section was wrong all the time... C so take it out (000426) CC IF (ACTION(2:2).EQ.'F') THEN CC STARTB(1) = CC + STEPB(1) * FLOAT( INT( (STARTB(1)-STARTC(1))/STEPB(1) ) ) CC + + STARTC(1) CC STARTB(2) = CC + STEPB(2) * FLOAT( INT( (STARTB(2)-STARTC(2))/STEPB(2) ) ) CC + + STARTC(2) CC ENDIF C C invert transformation matrix CALL INVERS(TOT(4),TRR,STAT) DO 1000, N=1,6 TR(N) = TRR(N) !R*8 -> R*4 1000 CONTINUE C C and map it... CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,1, + NAXISA,NPIXB,STARTB,STEPB, + IDENTA,CUNITA,PNTRB,IMNOB,STAT) C C and do the rotation CALL MATROT(MADRID(PNTRA),MADRID(PNTRB),NPIXA,NPIXB, + STARTA,END,STEPA,STARTB,STEPB,TR,RNULL) 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 INVERS(TRIN,TROUT,STAT) C IMPLICIT NONE C INTEGER STAT C DOUBLE PRECISION TRIN(6),TROUT(6) DOUBLE PRECISION DD C DD = TRIN(1)*TRIN(4) - TRIN(2)*TRIN(3) IF (DABS(DD).LT.1.D-20) THEN STAT = 1 RETURN ELSE STAT = 0 ENDIF C TROUT(1) = TRIN(4) / DD TROUT(2) = -TRIN(2) / DD TROUT(3) = -TRIN(3) / DD TROUT(4) = TRIN(1) / DD TROUT(5) = -( TROUT(1)*TRIN(5) + TROUT(2)*TRIN(6) ) TROUT(6) = -( TROUT(3)*TRIN(5) + TROUT(4)*TRIN(6) ) RETURN END SUBROUTINE MATROT(A,B,NPIXA,NPIXB,START,END,STEP, + STARTB,STEPB,TR,ZNULL) C C K. Banse 881028, 910205, 000419 C C do "matrix" rotation, i.e. compute: C C X = XROT*A1 + YROT*A2 + XC C Y = YROT*A3 + XROT*A4 + YC C C in real world coordinates C C IMPLICIT NONE C INTEGER NPIXA(2),NPIXB(2) INTEGER NX,NY,INDX INTEGER NO1,NO2,NO3,NO4,XDIMA,SIZEA INTEGER KX,KY,KSTX C REAL A(*),B(*) REAL ZNULL C DOUBLE PRECISION TR(6) DOUBLE PRECISION Y1,X,Y,DX,DY DOUBLE PRECISION B1,B2,B3,B4,XC,YC,BB1,BB3 DOUBLE PRECISION XT,YT,XFRAC,YFRAC,XST,YST DOUBLE PRECISION STARTB(2),STEPB(2),START(2),END(2),STEP(2) DOUBLE PRECISION MINI(2), MAXI(2) C C init B1 = TR(1) B2 = TR(2) B3 = TR(3) B4 = TR(4) XC = TR(5) YC = TR(6) DX = B1 * STEPB(1) DY = B3 * STEPB(1) XDIMA = NPIXA(1) SIZEA = NPIXA(1)*NPIXA(2) IF (START(1).LT.END(1)) THEN MINI(1) = START(1) MAXI(1) = END(1) ELSE MINI(1) = END(1) MAXI(1) = START(1) ENDIF C IF (START(2).LT.END(2)) THEN MINI(2) = START(2) MAXI(2) = END(2) ELSE MINI(2) = END(2) MAXI(2) = START(2) ENDIF C XST = 1.D0/STEP(1) YST = 1.D0/STEP(2) BB1 = B1 * STARTB(1) BB3 = B3 * STARTB(1) INDX = 0 C C loop through result lines Y1 = STARTB(2) !init Y1 to start-y DO 2000,NY=1,NPIXB(2) C X = BB1 + B2*Y1 + XC Y = BB3 + B4*Y1 + YC C DO 1000,NX=1,NPIXB(1) INDX = INDX + 1 !loop along result line IF ((X.LT.MINI(1)) .OR. (Y.LT.MINI(2)) .OR. + (X.GT.MAXI(1)) .OR. (Y.GT.MAXI(2))) THEN B(INDX) = ZNULL GOTO 500 ENDIF C C move to pixel space XT = ( (X-START(1)) * XST ) + 1.D0 !also adjust 0,0 to 1,1 YT = ( (Y-START(2)) * YST ) + 1.D0 C C interpolate... 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 int. part of real B(INDX) 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(INDX) = A(NO1) !yes. ELSE NO3 = NO1 + XDIMA !no. B(INDX) = A(NO1)*(1.D0-YFRAC) + A(NO3)*YFRAC ENDIF GOTO 500 ENDIF C NO3 = NO1 + XDIMA !get low pixel in next line up IF (NO3.GT.SIZEA) THEN !are we out of grid in y? B(INDX) = A(NO1)*(1.D0-XFRAC) + A(NO2)*XFRAC GOTO 500 ENDIF C NO4 = NO3 + 1 B(INDX) = ( A(NO1)*(1.D0-XFRAC) + A(NO2)*XFRAC ) + * (1.D0-YFRAC) + + + ( A(NO3)*(1.D0-XFRAC) + A(NO4)*XFRAC ) * YFRAC CC B(INDX) = A(NO1) + CC + XFRAC * (A(NO2) - A(NO1)) + CC + YFRAC * (A(NO3) - A(NO1)) + CC + XFRAC * YFRAC * (A(NO1) - A(NO2) - A(NO3) + A(NO4)) C 500 X = X + DX Y = Y + DY 1000 CONTINUE C Y1 = Y1 + STEPB(2) !move to next y-coord. 2000 CONTINUE C RETURN END