C @(#)rectimag.for	19.1 (ESO-DMD) 02/25/03 13:31:23
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
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION: RECTIMAG
C.AUTHOR:         D. Baade
C.VERSION:        850809  version 1.2
C.VERSION;        850820  K. Banse  
C.KEYWORDS:       Geometrical rectification of 2-D images, rebinning, 
C                 flux conservation
C.PURPOSE:        Geometrically rectificy distorted 2-D images, rebin but 
C                 roughly maintain the original step sizes, take care that 
C                 flux is conserved
C.LANGUAGE:       F77+ESOext
C.ALGORITHM:      a) rectification: 
C                    The parameters of a bipolynomial regression of 
C		     degree (3,3) have to be calculated first by 
C	             REGRESSION/POLY on TABLE data (e.g. with positions of 
C                    reseaux marks).
C                    They are stored in the keyword COEFFX (1,...,20) and
C                    COEFFY (1,...,20).
C                 b) rebinning:
C                    Each pixel is routinely subdivided into 9 subpixels
C                    (unless SUBFAC=1.). Assuming that
C                    the coordinates of the input pixel refer to its center
C                    and that its contents is the average flux integrated over
C                    its area (no dead space), a simple linear interpolation
C                    with nearest neighboring original pixels is done for each
C                    of the 9 subpixels. After this interpolation the flux
C                    summed over the subpixels is scaled to the flux of the
C                    original pixel so that flux conservation is achieved on a
C                    pixel basis. Finally, each subpixel is imaged to the
C                    rectified coordinate frame and its contents added to the
C                    appropriate pixel(s) (of about the same size as in the 
C                    input frame). Optionally (NREP.GT.1) each subpixel can 
C                    first be subdivided into still smaller units. Sub-subpixels
C                    originating from the same parent-subpixel all contain the
C                    same flux, i.e. no further interpolation is done.
C.RESTRICTIONS:      a) works in pixel space only, i.e STEP(x,y) = (+1,+1)
C                       in BOTH INPUT and OUTPUT frame
C                    b) it is implicitly assumed that the positions of the
C                       reseaux marks, too, are meassured in these same units
C.INPUT/OUTPUT:      The following keywords are used:
C                    IN_A/C/1/60     name of input frame
C                    OUT_A/C/1/60    name of output frame
C                    COEFFX/D/1/20   regression coefficients for X-coordinate
C                    COEFFY/D/1/20   regression coefficients for Y-coordinate
C                    INPUTI/I/1/1    repetition factor (for the further 
C                                    subdivision of subpixe
C                    INPUTR/R/1/1    substepping factor
C.VERSION:           871123   Rein H. Warmels  ESO-FORTRAN Conversion
C 
C 010625		last modif
C 
C -------------------------------------------------------------------
      PROGRAM      RCTFRA
C
      IMPLICIT     NONE
C
      INTEGER      MADRID(1)
      INTEGER      I,IAV
      INTEGER      KUN,KNUL
      INTEGER      AIMF, CIMF, DIMF
      INTEGER      NAXIS,NPIX(2),ISTAT
      INTEGER      NREP,OUTPIX(2),NPIXIP(2)
C 
      INTEGER*8    IPNTRA,IPNTRC,IPNTRD
C
      REAL         OUTSTA(2),SUBFAC
C 
      DOUBLE PRECISION STEP(2),START(2)
      DOUBLE PRECISION COEFFX(20),COEFFY(20)
C 
      CHARACTER    FRAMEA*60,FRAMEC*60
      CHARACTER    IDENT*72,CUNIT*64
C
      COMMON        /VMR/MADRID
C 
      INCLUDE       'MID_INCLUDE:ST_DEF.INC'
C 
      DATA          IDENT /' '/, CUNIT /' '/
C 
      INCLUDE       'MID_INCLUDE:ST_DAT.INC'
C
C *** begin code
C
      CALL STSPRO('RECTIMAG')                       ! get into midas environment
      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUN,KNUL,ISTAT) ! name of input frame
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2,
     2            NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTRA,
     +            AIMF,ISTAT)                                       ! map frame
C
C *** insist that frame be two-dimensional
      IF (NAXIS.NE.2) 
     +   CALL STETER(1,'*** FATAL: Only 2-d frames allowed ... ')
C
C *** get regression coefficients, repetition factor and substepping flag
      CALL STKRDD('COEFFX',1,20,IAV,COEFFX,KUN,KNUL,ISTAT)
      CALL STKRDD('COEFFY',1,20,IAV,COEFFY,KUN,KNUL,ISTAT)
      CALL STKRDI('INPUTI',1,1,IAV,NREP,KUN,KNUL,ISTAT)
      CALL STKRDR('INPUTR',1,1,IAV,SUBFAC,KUN,KNUL,ISTAT)
C
C *** limit nrep to 5
      IF (NREP.GT.5) THEN
         NREP   = 5
         CALL STTPUT('*** WARNING: Substepping factor set to 5!',
     2               ISTAT)
      END IF
C
C *** calculate values of outsta and outpix for output frame
      CALL OUTDIM(START,STEP,NPIX,COEFFX,COEFFY,OUTSTA,OUTPIX)
C
C *** get name of output frame, map frame
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,KUN,KNUL,ISTAT)
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     2            NAXIS,OUTPIX,OUTSTA,STEP,IDENT,CUNIT,IPNTRC,
     3            CIMF,ISTAT)
C
C *** and initialize its contents to zero
      CALL INTNUL(OUTPIX(1)*OUTPIX(2),MADRID(IPNTRC))
C
C *** copy descriptor lhcuts because the cuts shouldn't change much ...
      CALL STDCOP(AIMF,CIMF,4,'LHCUTS',ISTAT)
C
C *** let subroutines rectinterpol, rectscale, and rectmap do the work
      IF (SUBFAC.GT.1.1) THEN                        ! deconvolution to be done
C
C *** create dummy file to hold intermediate results after deconvolution
         DO 10,I=1,2
            NPIXIP(I) = NPIX(I)*SUBFAC
   10    CONTINUE
         CALL STIPUT('IPOLDUM',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     2            NAXIS,NPIXIP,START,STEP,IDENT,CUNIT,IPNTRD,
     3            DIMF,ISTAT)
         CALL INTNUL(NPIXIP(1)*NPIXIP(2),MADRID(IPNTRD))
         CALL RCTINT(MADRID(IPNTRA),NPIX(1),NPIX(2),MADRID(IPNTRD))
         CALL RCTMAP(MADRID(IPNTRC),OUTPIX(1),OUTPIX(2),START,
     2               MADRID(IPNTRD),NPIXIP(1),NPIXIP(2),COEFFX,COEFFY,
     3               OUTSTA,SUBFAC,NREP)
C  
      ELSE                                                  ! (no deconvolution)
         CALL RCTMAP(MADRID(IPNTRC),OUTPIX(1),OUTPIX(2),START,
     2               MADRID(IPNTRA),NPIX(1),NPIX(2),COEFFX,COEFFY,
     3               OUTSTA,SUBFAC,NREP)
      END IF
C
      CALL STSEPI
      END
