
CD-Section: In Memory Image Routines
C The routine described in this section is for extracting a section of 
C an image which has been stored in memory. Its use is associated with
C the development of the Laue Data Module (LDM) routines though it may
C be used independently of these.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LFN_IMG_RECT     **
C**************************
C
CD-Routine: Read Section from Image in Memory - LFN_IMG_RECT
C This routine is used to read in a rectangular section of an 'in memory' 
C image into an integer array. The data in the image itself may be stored in
C byte format, unsigned two-byte format, signed integer format or a
C 'squashed' two byte integer format.  In the 'squashed i2' format an
C intensity>32767 is as 65536-intensity/8.
CD-end
CD-Fortran:
      SUBROUTINE LFN_IMG_RECT (IMG, ITYPE, IORDER, NF_OFF, 
     +                         IAX1_ORIG, IAX2_ORIG, NAX1, NAX2,
     +                         NAX1_RASTS, NAX2_RASTS, IOUTLIM,
     +                         IDATA, IWORK, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Get a rectangular section from an image stored in memory
C
C Author:  John W. Campbell, October 1993
C
C Arguments:
C
      INTEGER IMG(*)
      INTEGER ITYPE
      INTEGER IORDER
      INTEGER NF_OFF
      INTEGER IAX1_ORIG
      INTEGER IAX2_ORIG
      INTEGER NAX1
      INTEGER NAX2
      INTEGER NAX1_RASTS
      INTEGER NAX2_RASTS
      INTEGER IOUTLIM
      INTEGER IDATA (NAX2,NAX1)
      INTEGER IWORK(*)
      INTEGER IERR
CD-Parameters:
C   IMG (*)    i (R)  Image data array (may hold signed integer data or
C                     unsigned byte data packed into an integer array or
C                     unsigned two byte data packed into an integer array
C                     or 'squashed' two byte data packed into an integer
C                     array.
C   ITYPE      i (R)  Image data type =1 unsigned byte
C                                     =2 unsigned two byte
C                                     =3 signed (full) integer
C                                     =4 'squashed i2' (if intensity>32767
C                                        then stored as 65536-intensity/8)
C   IORDER     i (R)  Order of data in image wrt. the two local axes
C                     =1 +ax1 +ax2
C                     =2 +ax1 -ax2
C                     =3 -ax1 +ax2
C                     =4 -ax1 -ax2
C                     =5 +ax2 +ax1
C                     =6 +ax2 -ax1
C                     =7 -ax2 +ax1
C                     =8 -ax2 -ax1
C   NF_OFF     i (R)  Fast offset between (slow) rasters of image (e.g. 
C                     first dimension of the image data if the data are 
C                     held in a two dimensional array). This value is the 
C                     number of elements for the image data type (bytes, 
C                     2-bytes, integers or 'squashed i2' as set by ITYPE)
C   IAX1_ORIG  i (R)  The origin pixel of the required section along the 
C                     first local axis
C   IAX2_ORIG  i (R)  The origin pixel of the required section along the 
C                     second local axis
C   NAX1       i (R)  Size of the section in pixels along the first local
C                     axis
C   NAX2       i (R)  Size of the section in pixels along the second local
C                     axis
C   NAX1_RASTS i (R)  No. of rasters in the image along the first 
C                     local axis    
C   NAX2_RASTS i (R)  No. of rasters in the image along the second 
C                     local axis
C   IOUTLIM    i (R)  Value to which returned data pixels are to be set if
C                     they lie outside the image
C   IDATA      i (W)  Array returning the image section data first local
C                     axis is the slow axis (2'nd dimension) and second
C                     local axis is the fast axis (1'st) dimension.
C   IWORK      i (*)  Work array dimensioned to at least NAX2 if IORDER
C                     = 1-4 or NAX1 if IORDER = 5-8; may be a dummy if
C                     ITYPE=3
C   IERR       i (W)  Error flag = 0 OK,
C                              =-1 Invalid type or order flag
C                              =-2 Byte handling not available
C                              = 1 Box partly outside image
C                              = 2 Box totally outside image
C
C              Note: pixels within the image are numbered from 1 upwards
C
C              Note: byte handling (via CCP4 library routines) must be
C                    available for data types 1 and 2
CD-end
C
C====== FUNCTION TYPES:
C
      LOGICAL CCPBYT
C
C====== LOCALS:
C
      INTEGER NBW
      INTEGER ISLOW_MIN
      INTEGER ISLOW_MAX
      INTEGER IFAST_MIN
      INTEGER IFAST_MAX
      INTEGER NS_RASTS
      INTEGER NF_RASTS
      INTEGER IAX1
      INTEGER IAX1_START
      INTEGER IAX1_END
      INTEGER IAX1_STEP
      INTEGER IAX2
      INTEGER IAX2_START
      INTEGER IAX2_END
      INTEGER IAX2_STEP
      INTEGER KS
      INTEGER KF
      INTEGER KOFF
      INTEGER MINEL
      INTEGER MAXEL
      INTEGER NELS
      INTEGER KTYP

C     INTEGER NBW               !No. bytes/word
C     INTEGER ISLOW_MIN         !Minimum pixel along image slow axis
C     INTEGER ISLOW_MAX         !Maximum pixel along image slow axis
C     INTEGER IFAST_MIN         !Minimum pixel along image fast axis
C     INTEGER IFAST_MAX         !Maximum pixel along image fast axis
C     INTEGER NS_RASTS          !No. of rasters in image along slow direction
C     INTEGER NF_RASTS          !No. of rasters in image along fast direction
C     INTEGER IAX1              !Loop through local axis 1
C     INTEGER IAX1_START        !Start of loop through local axis 1
C     INTEGER IAX1_END          !End of loop through local axis 1
C     INTEGER IAX1_STEP         !Loop step for local axis 1
C     INTEGER IAX2              !Loop through local axis 2
C     INTEGER IAX2_START        !Start of loop through local axis 2
C     INTEGER IAX2_END          !End of loop through local axis 2
C     INTEGER IAX2_STEP         !Loop step through local axis 2
C     INTEGER KS                !Counter for slow axis
C     INTEGER KF                !Counter for fast axis
C     INTEGER KOFF              !Offset to start of current fast axis raster
C     INTEGER MINEL             !Minimum element number for data translation
C     INTEGER MAXEL             !Maximum element number for data translation
C     INTEGER NELS              !No. of elements for data translation
C     INTEGER KTYP              !Data type flag for CCPCPI call
C
C-------------------------------------------------------------------------------
C
      IERR = 0 
C
C====== Check for valid flag values
C
      IF (ITYPE.LT.1.OR.ITYPE.GT.4) THEN
         IERR = -1
         RETURN
      END IF
      IF (IORDER.LT.1.OR.IORDER.GT.8) THEN
         IERR = -1
         RETURN
      END IF
C
C====== Check that byte handling available if requested
C
      IF (ITYPE.NE.3) THEN
         IF (.NOT.CCPBYT(NBW))THEN
            IERR = -2
            RETURN
         END IF
      END IF
C
C====== Get Box limits in image array and loop parameters for returning
C       data array
C
      IF (IORDER.EQ.1) THEN
         ISLOW_MIN = IAX1_ORIG
         IFAST_MIN = IAX2_ORIG
         IAX1_START = 1
         IAX1_END = NAX1
         IAX1_STEP = 1
         IAX2_START = 1
         IAX2_END = NAX2
         IAX2_STEP = 1
      ELSE IF (IORDER.EQ.2) THEN
         ISLOW_MIN = IAX1_ORIG
         IFAST_MIN = NAX2_RASTS - IAX2_ORIG - NAX2 + 2
         IAX1_START = 1
         IAX1_END = NAX1
         IAX1_STEP = 1
         IAX2_START = NAX2
         IAX2_END = 1
         IAX2_STEP = -1
      ELSE IF (IORDER.EQ.3) THEN
         ISLOW_MIN = NAX1_RASTS - IAX1_ORIG - NAX1 + 2
         IFAST_MIN = IAX2_ORIG
         IAX1_START = NAX1
         IAX1_END = 1
         IAX1_STEP = -1
         IAX2_START = 1
         IAX2_END = NAX2
         IAX2_STEP = 1
      ELSE IF (IORDER.EQ.4) THEN
         ISLOW_MIN = NAX1_RASTS - IAX1_ORIG - NAX1 + 2
         IFAST_MIN = NAX2_RASTS - IAX2_ORIG - NAX2 + 2
         IAX1_START = NAX1
         IAX1_END = 1
         IAX1_STEP = -1
         IAX2_START = NAX2
         IAX2_END = 1
         IAX2_STEP = -1
      ELSE IF (IORDER.EQ.5) THEN
         ISLOW_MIN = IAX2_ORIG
         IFAST_MIN = IAX1_ORIG
         IAX2_START = 1
         IAX2_END = NAX2
         IAX2_STEP = 1
         IAX1_START = 1
         IAX1_END = NAX1
         IAX1_STEP = 1
      ELSE IF (IORDER.EQ.6) THEN
         ISLOW_MIN = IAX2_ORIG
         IFAST_MIN = NAX1_RASTS - IAX1_ORIG - NAX1 + 2
         IAX2_START = 1
         IAX2_END = NAX2
         IAX2_STEP = 1
         IAX1_START = NAX1
         IAX1_END = 1
         IAX1_STEP = -1
      ELSE IF (IORDER.EQ.7) THEN
         ISLOW_MIN = NAX2_RASTS - IAX2_ORIG - NAX2 + 2
         IFAST_MIN = IAX1_ORIG
         IAX2_START = NAX2
         IAX2_END = 1
         IAX2_STEP = -1
         IAX1_START = 1
         IAX1_END = NAX1
         IAX1_STEP = 1
      ELSE IF (IORDER.EQ.8) THEN
         ISLOW_MIN = NAX2_RASTS - IAX2_ORIG - NAX2 + 2
         IFAST_MIN = NAX1_RASTS - IAX1_ORIG - NAX1 + 2
         IAX2_START = NAX2
         IAX2_END = 1
         IAX2_STEP = -1
         IAX1_START = NAX1
         IAX1_END = 1
         IAX1_STEP = -1
      END IF
      IF (IORDER.LT.5) THEN
         ISLOW_MAX = ISLOW_MIN + NAX1 - 1
         IFAST_MAX = IFAST_MIN + NAX2 - 1
         NS_RASTS = NAX1_RASTS
         NF_RASTS = NAX2_RASTS
      ELSE
         ISLOW_MAX = ISLOW_MIN + NAX2 - 1
         IFAST_MAX = IFAST_MIN + NAX1 - 1
         NS_RASTS = NAX2_RASTS
         NF_RASTS = NAX1_RASTS
      END IF
C
C====== If some or all of requested section lies outside image then
C       initialise returned data array to the given out of limits
C       value
C
      IF (ISLOW_MIN.LT.1.OR.ISLOW_MAX.GT.NS_RASTS.
     +   OR.IFAST_MIN.LT.1.OR.IFAST_MAX.GT.NF_RASTS) THEN
         IERR = 1
         DO 50 IAX1 = 1, NAX1
            DO 40 IAX2 = 1, NAX2
               IDATA(IAX2,IAX1) = IOUTLIM
40          CONTINUE
50       CONTINUE
         IF (ISLOW_MIN.GT.NS_RASTS.OR.ISLOW_MAX.LT.1) THEN
            IERR = 2
            RETURN
         END IF
         IF (IFAST_MIN.GT.NF_RASTS.OR.IFAST_MAX.LT.1) THEN
            IERR = 2
            RETURN
         END IF
      END IF
C
C======  Signed (full) integer data
C
      IF (ITYPE.EQ.3) THEN
         IF (IORDER.LT.5) THEN
            KS = ISLOW_MIN - 1 
            KOFF = (KS-1)*NF_OFF
            DO 120 IAX1 = IAX1_START, IAX1_END, IAX1_STEP
               KS = KS + 1
               KOFF = KOFF + NF_OFF
               IF (KS.LT.1.OR.KS.GT.NS_RASTS) GO TO 120
               KF = IFAST_MIN - 1
               DO 100 IAX2 = IAX2_START, IAX2_END, IAX2_STEP
                  KF = KF + 1
                  IF (KF.LT.1.OR.KF.GT.NF_RASTS) GO TO 100
                  IDATA(IAX2,IAX1) = IMG(KOFF+KF)
100            CONTINUE
120         CONTINUE
         ELSE
            KS = ISLOW_MIN - 1 
            KOFF = (KS-1)*NF_OFF
            DO 220 IAX2 = IAX2_START, IAX2_END, IAX2_STEP
               KS = KS + 1
               KOFF = KOFF + NF_OFF
               IF (KS.LT.1.OR.KS.GT.NS_RASTS) GO TO 220
               KF = IFAST_MIN - 1
               DO 200 IAX1 = IAX1_START, IAX1_END, IAX1_STEP
                  KF = KF + 1
                  IF (KF.LT.1.OR.KF.GT.NF_RASTS) GO TO 200
                  IDATA(IAX2,IAX1) = IMG(KOFF+KF)
200            CONTINUE
220         CONTINUE
         END IF
         RETURN
      ELSE
C
C====== Unsigned byte or unsigned integer*2 data
C       
         IF (ITYPE.EQ.1) THEN
            KTYP = 1
         ELSE 
            KTYP = 3
         END IF
         MINEL = IFAST_MIN
         IF (MINEL.LT.1) MINEL = 1
         MAXEL = IFAST_MAX
         IF (MAXEL.GT.NF_RASTS) MAXEL = NF_RASTS
         NELS = MAXEL - MINEL + 1
         IF (IORDER.LT.5) THEN
            KS = ISLOW_MIN - 1
            KOFF = (KS-1)*NF_OFF
            DO 320 IAX1 = IAX1_START, IAX1_END, IAX1_STEP
               KS = KS + 1
               KOFF = KOFF + NF_OFF
               IF (KS.LT.1.OR.KS.GT.NS_RASTS) GO TO 320
               CALL CCPCPI (IWORK, IMG, KOFF+MINEL, KOFF+MAXEL, KTYP)
               KF = IFAST_MIN - MINEL
               DO 300 IAX2 = IAX2_START, IAX2_END, IAX2_STEP
                  KF = KF + 1
                  IF (KF.LT.1.OR.KF.GT.NELS) GO TO 300
                  IF (ITYPE.EQ.4.AND.IWORK(KF).GT.32767)
     +               IWORK(KF) = 8*(65536-IWORK(KF))
                  IDATA(IAX2,IAX1) = IWORK(KF)
300            CONTINUE
320         CONTINUE
         ELSE
            KS = ISLOW_MIN - 1
            KOFF = (KS-1)*NF_OFF
            DO 420 IAX2 = IAX2_START, IAX2_END, IAX2_STEP
               KS = KS + 1
               KOFF = KOFF + NF_OFF
               IF (KS.LT.1.OR.KS.GT.NS_RASTS) GO TO 420
               CALL CCPCPI (IWORK, IMG, KOFF+MINEL, KOFF+MAXEL, KTYP)
               KF = IFAST_MIN - MINEL
               DO 400 IAX1 = IAX1_START, IAX1_END, IAX1_STEP
                  KF = KF + 1
                  IF (KF.LT.1.OR.KF.GT.NELS) GO TO 400
                  IF (ITYPE.EQ.4.AND.IWORK(KF).GT.32767)
     +               IWORK(KF) = 8*(65536-IWORK(KF))
                  IDATA(IAX2,IAX1) = IWORK(KF)
400            CONTINUE
420         CONTINUE
         END IF
         RETURN
      END IF
      END
