 
CD-Section: Determine Spot Positions and Sizes
C A set of routines is available to determine the scanner coordinates of
C of the spots on a Laue image, to determine the spot lengths and widths
C and to analyse the results. A routine to output the results is also
C available.
CD-end
C
C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************
C**     LFN_FINDSPOTS      **
C****************************
C
CD-Routine: Find Spot Positions in an Image - LFN_FINDSPOTS
C This routine is used to find the Laue spot positions in an image
C stored in memory. A pre-calculated background image is passed to the
C routine and the routine searches for pixels at a certain threshols
C above background to determine those pixels which are presumed to belong
C to spots. The routine continuously updates the spot position
C by updating the calculation of its centre of gravity as the spot is
C built up from contiguous above threshold pixels. At the same time as the
C spot is being built up, a rotated image of the spot is also being 
C built up along the horizontal axis so that the spot length is just
C determined as the horizontal extent of the finished spot and the spot
C width as the vertical extent of the rotated spot. A list of spot details
C is stored in an in-memory spot list managed by the 'DML' set of routines.
CD-end
CD-Fortran:
      SUBROUTINE LFN_FINDSPOTS (IMG, ITYPE, IORDER, NF_OFF, 
     +                          NAX1_RASTS, NAX2_RASTS,
     +                          IMG_BG, NF_BGOFF, NBGCMP, 
     +                          IDATA, IWORK,
     +                          XCEN, YCEN, RMIN, RMAX, 
     +                          XLOW, XHIGH, YLOW, YHIGH, REXCLUDE, 
     +                          THRESHFAC, PROGRESS, CANCEL,
     +                          MINDX_SPOT, NSPOT, IERR, 
     +                          ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C PURPOSE: Find spot positions for a Laue image
C
C AUTHOR: D. Nguti
C
C ARGUMENTS:
C
      INTEGER IMG(*)
      INTEGER ITYPE
      INTEGER IORDER
      INTEGER NF_OFF
      INTEGER NAX1_RASTS
      INTEGER NAX2_RASTS
      INTEGER IMG_BG(*)
      INTEGER NF_BGOFF
      INTEGER NBGCMP
      INTEGER IDATA(*)
      INTEGER IWORK(*)
      REAL XCEN
      REAL YCEN
      REAL RMIN
      REAL RMAX
      REAL XLOW
      REAL XHIGH
      REAL YLOW
      REAL YHIGH
      REAL REXCLUDE
      REAL THRESHFAC
      EXTERNAL PROGRESS
      LOGICAL CANCEL
      EXTERNAL CANCEL
      INTEGER MINDX_SPOT
      INTEGER NSPOT
      INTEGER IERR
      INTEGER KERR
      CHARACTER*(*) ERRSTR
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                     A fourth 'squashed' i2 format is also available.
C   ITYPE      i (R)  Image data type =1 unsigned byte
C                                     =2 unsigned two byte
C                                     =3 signed (full) integer
C                                     =4 'squashed i2'
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
C                     (e.g. first dimension of the image data if the data
C                     are held in a two dimensional array). This value is
C                     the number of elements for the image data type
C                     (bytes, 2-bytes or integers as set by ITYPE)
C   IMG_BG     i (R)  Background image (same data format as IMG)
C   NF_BGOFF   i (R)  Fast offset between (slow) rasters of background 
C                     image (e.g. first dimension of the background image 
C                     data if the data are held in a two dimensional array).
C                     This value is the number of elements for the image 
C                     data type (bytes, 2-bytes or integers as set by 
C                     ITYPE)
C   NBGCMP     i (R)  No. of pixels in each direction which were compressed
C                     into 1 for the background calculation.
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   IDATA      i (*)  Work array to hold the spot data returned from the
C                     image. Must be dimensioned to at least the maximum
C                     no. of rasters in either of the two dimensions of
C                     the full image.
C   IWORK      i (*)  Work array dimensioned to at least the maximum
C                     no. of rasters in either of the two dimensions of
C                     the full image.
C   XCEN       r (R)  'x' centre of pattern image in rasters (corrected)
C   YCEN       r (R)  'y' centre of pattern image in rasters (corrected)
C   RMIN       r (R)  Minimum radius of active area (rasters)
C   RMAX       r (R)  Maximum radius of active area (rasters)
C   XLOW       r (R)  Low limit on 'x' of area to include (0.0=image min.)
C                     (rasters)
C   XHIGH      r (R)  High limit on 'x' of area to include (0.0=image max.)
C                     (rasters)
C   YLOW       r (R)  Low limit on 'y' of area to include (0.0=image min.)
C                     (rasters)
C   YHIGH      r (R)  High limit on 'y' of area to include (0.0=image max.)
C                     (rasters)
C   REXCLUDE   r (R)  The routine uses RMIN+REXCLUDE and RMAX-REXCLUDE as 
C                     actual limits for min and max radii to avoid edge of
C                     detector problems with background calculation 
C                     (rasters). Also uses XLOW+REXCLUDE, XHIGH-REXCLUDE, 
C                     YLOW+REXCLUDE and YHIGH-REXCLUDE as limits if 
C                     XLOW, XHIGH etc. are set. A suggestion is to use 
C                     a value value of (dimension in pixels of background
C                     calculation box)/sqrt(2) if XDLF_BG_CALC was used
C                     to calculate the background image.
C   THRESHFAC  i (R)  Uses pixels more than THRESHFAC*background above 
C                     background as being pixels which belong to spots.
C   PROGRESS     ( )  Subroutine which will be called at intervals during
C                     the spots find which may be used to display its 
C                     progress if required (See below for its arguments).
C                     Note that the routine must be declared as EXTERNAL 
C                     in the calling program.
C   CANCEL     f (R)  Logical function used to interrupt procedure.
C                     It has one dummy argument and must return a
C                     logical value of .TRUE. to cancel the procedure
C                     or .FALSE. to let it proceed.
C                     The actual function passed must be declared as
C                     both LOGICAL and EXTERNAL in the calling program.
C                     (The dummy function LFN_DUMMY_CNCL may be used)
C   MINDX_SPOT i (W)  'DML' record storage routines index.
C                     Spots list returned. Seven words/spot
C                     Word 1  SX    Spot 'x' position (rasters)   (real)
C                     Word 2  SY    Spot 'y' position (rsaters)   (real)
C                     Word 3  SUMI  Spot intensity                (real)
C                     Word 4  NPIX  No. of pixels in spot         (integer)
C                     Word 5  SPOT_LENGTH Spot length (pixels)    (integer)
C                     Word 6  SPOT_WIDTH  Spot width (pixels)     (integer)
C                     Word 7  SPOT_ANGLE  Spot angle in degrees   (real)
C   NSPOT      i (W)  No. of spots found.
C   IERR       i (W)  Error flag =  0, OK
C                                =  1, Memory allocation error for spots 
C                                      list
C                                =  2, Too many active spots
C                                =100, Cancelled
C   ERRSTR     c (W)  Error message string (max 80 chars) set if IERR
C                     is non-zero.
C
C                    Subroutine PROGRESS (IFLAG, NSTRIPS, ISTRIP)
C
C                    This routine called by the spot find routine 
C                    with the following parameters:-
C
C                        IFLAG  i (R) Integration pass = 0 start
C                                                      = 1 spot search in
C                                                          progress
C                                                      = 2 finished
C                        NSTRIP i (R) No. strips for spot search
C                        ISTRIP i (W) Current strip
C
CD-end
C
C====== Parameter statements
C
 
C
C====== GLOBALS:
C
      INCLUDE 'lfn_findspots.finc'
C
C====== FUNCTION TYPES:
C
 
C
C====== LOCALS:
C
      INTEGER I, J, L
      INTEGER IBGVAL(2)
      INTEGER ITHRESH
      INTEGER IBL
      INTEGER IEL
      LOGICAL TEST1
      LOGICAL LINE
      INTEGER IBUF(7)
      INTEGER MXCENT
      INTEGER IRMIN
      INTEGER IRMINSQ
      INTEGER IRMAX
      INTEGER IXCEN
      INTEGER IYCEN
      INTEGER IXL
      INTEGER IXH
      INTEGER IYL
      INTEGER IYH
      INTEGER IXRAST
      INTEGER IDX
      INTEGER IDELX
      INTEGER IXBG
      INTEGER IYPT
      INTEGER IYBG
      INTEGER ILEN
      INTEGER ICENT
      INTEGER NCENT
      INTEGER IY
      INTEGER NSTRIPS
      INTEGER ISTRIP
      INTEGER KXMIN
      INTEGER KXMAX
      INTEGER KYMIN
      INTEGER KYMAX
      REAL PI
      REAL DEG
      REAL XSQ
      REAL FXVAL
      REAL SUMYOD
      REAL SUMOD
      REAL X_ORIGIN
      REAL Y_ORIGIN
      REAL X_ROT
      REAL Y_ROT
      REAL R1
      REAL SX
      REAL SY
      REAL SPOT_LENGTH
      REAL SPOT_WIDTH
      REAL SPOT_ANGLE
      REAL SUMINT
      INTEGER NPIX
      EQUIVALENCE (SX,IBUF(1)), (SY,IBUF(2)), (SUMINT,IBUF(3))
      EQUIVALENCE (NPIX,IBUF(4)), (SPOT_LENGTH,IBUF(5))
      EQUIVALENCE (SPOT_WIDTH,IBUF(6)), (SPOT_ANGLE,IBUF(7))
C
C-------------------------------------------------------------------------------
C
C====== Initialisations
C
      NCENT = 0
      NSPOT = 0
      IERR = 0
      ERRSTR = ' '
      MXCENT = MAX_ACT_SPOTS
      PI = ACOS(-1.0)
      DEG = 180.0/PI
      KXMIN = 1
      IF (XLOW.GT.0.0) KXMIN = NINT(XLOW)
      KXMAX = NAX1_RASTS
      IF (XHIGH.GT.0.0) KXMAX = NINT(XHIGH)
      KYMIN = 1
      IF (YLOW.GT.0.0) KYMIN = NINT(YLOW)
      KYMAX = NAX2_RASTS
      IF (YHIGH.GT.0.0) KYMAX = NINT(YHIGH)
      KXMIN = NINT(KXMIN+REXCLUDE)
      KXMAX = NINT(KXMAX-REXCLUDE)
      KYMIN = NINT(KYMIN+REXCLUDE)
      KYMAX = NINT(KYMAX-REXCLUDE)
C
C====== Create output spots list
C
      CALL DMLF_INIT(500,500,7,MINDX_SPOT)
      IF (MINDX_SPOT.LT.0) THEN
         IERR = 1
         ERRSTR = 'Cannot allocate memory for spots list'
         RETURN
      END IF
C
      IRMAX = NINT(RMAX-REXCLUDE)
      IRMIN = NINT(RMIN+REXCLUDE)
      IRMINSQ = IRMIN*IRMIN
C
C====== Calculate limits along first local axis
C
      IXCEN = NINT(XCEN)
      IYCEN = NINT(YCEN)
      IXL = MAX(1,IXCEN-IRMAX,KXMIN)
      IXH = MIN(NAX1_RASTS,IXCEN+IRMAX,KXMAX)
      NSTRIPS = IXH - IXL + 1
      ISTRIP = 0
      CALL PROGRESS (0, NSTRIPS, 0)
C
C====== Loop over stripes
C
      DO 130 IXRAST = IXL , IXH
         IDX = IXCEN - IXRAST
         XSQ = IDX*IDX
C
C====== Determine limits on second local axis
C
         IDELX = NINT(SQRT(REAL(IRMAX*IRMAX-XSQ)))
         IYL = MAX(1,IYCEN-IDELX,KYMIN)
         IYH = MIN(NAX2_RASTS,IYCEN+IDELX,KYMAX)

         CALL LFN_IMG_RECT(IMG,ITYPE,IORDER,NF_OFF,IXRAST,1,1,
     +                     NAX2_RASTS,NAX1_RASTS,NAX2_RASTS,0,IDATA,
     +                     IWORK,KERR)
C
C====== Radius within which spots are picked
C
C====== Loop over pixels to find lines of continuous threshold density
C
         FXVAL = REAL(IXRAST)
         LINE = .FALSE.
         IXBG = (IXRAST-1)/NBGCMP + 1
         DO 70 IYPT = IYL , IYH
            IYBG = (IYPT-1)/NBGCMP + 1
 
            CALL LFN_IMG_RECT(IMG_BG,ITYPE,IORDER,NF_BGOFF,IXBG,IYBG,1,
     +                        1,NAX1_RASTS/NBGCMP,NAX2_RASTS/NBGCMP,0,
     +                        IBGVAL,IWORK,KERR)
 
C            ITHRESH = IDATA(IYPT) - IBGVAL(1) - KTHRESH
            ITHRESH = IDATA(IYPT) - THRESHFAC*IBGVAL(1)
            IF (ABS(IXCEN-IXRAST).LE.IRMIN) THEN
               IF ((IYPT-IYCEN)*(IYPT-IYCEN)
     +            +(IXRAST-IXCEN)*(IXRAST-IXCEN).LE.IRMINSQ)
     +            ITHRESH = -1
            END IF
C
            IF (.NOT.LINE) THEN
               IF (ITHRESH.LE.0) GO TO 70
               LINE = .TRUE.
               IBL = IYPT
               SUMYOD = REAL(ITHRESH*IYPT)
               SUMOD = REAL(ITHRESH)

            ELSE IF (ITHRESH.LE.0 .OR. IYPT.EQ.IYH) THEN
               LINE = .FALSE.
               IEL = IYPT - 1
               ILEN = IEL - IBL + 1
C
C====== Inner loop to decide which spots are ovelapping
C
               DO 30 ICENT = 1 , NCENT
                  J = IPT(ICENT)
C
C====== Skipping condition
C
                  TEST1 = IEL.LT.IYBEG(J) .OR. IBL.GT.IYEND(J)
                  IF (TEST1) GO TO 30

                  SUMXI(J) = SUMXI(J) + SUMOD*FXVAL
                  SUMYI(J) = SUMYI(J) + SUMYOD
                  SUMI(J) = SUMI(J) + SUMOD
                  IXCENT(J) = NINT(SUMXI(J)/SUMI(J))
                  IYCENT(J) = NINT(SUMYI(J)/SUMI(J))
                  NPIXCENT(J) = NPIXCENT(J) + ILEN
                  IYBEG(J) = IBL
                  IYEND(J) = IEL
                  IXMAX(J) = IXRAST
C
C====== Rotate pixels in strip by azimuthal angle
C
                  DO 20 IY = IBL , IEL
                     X_ORIGIN = IXRAST - IXCEN
                     Y_ORIGIN = IY - IYCEN
                     CALL LFN_ROTATE_XY(AZIMUTHAL_ANGLE(J),X_ORIGIN,
     +                              Y_ORIGIN,X_ROT,Y_ROT)

                     SPOT_LENGTH_MIN(J) = MIN(X_ROT,SPOT_LENGTH_MIN(J))
                     SPOT_LENGTH_MAX(J) = MAX(X_ROT,SPOT_LENGTH_MAX(J))

                     SPOT_WIDTH_MIN(J) = MIN(Y_ROT,SPOT_WIDTH_MIN(J))
                     SPOT_WIDTH_MAX(J) = MAX(Y_ROT,SPOT_WIDTH_MAX(J))
 20               CONTINUE
C
C====== Flag spot if line overlaps search area limits
C
                  IF (IBL.EQ.IYL .OR. IYPT.EQ.IYH) IFLAG(J) = 1
                  GO TO 70
 30            CONTINUE
C
C====== Count this active spot
C
               NCENT = NCENT + 1
               IF (NCENT.GT.MXCENT) THEN
                  IERR = 2
                  ERRSTR = 'Too many active spots, maximum is'
                  CALL STR_ADDI(ERRSTR,MXCENT,1)
                  GO TO 140
               END IF
C
C====== Find empty slot in which to store info about this spot
C
               DO 40 I = 1 , MXCENT
                  IF (NPIXCENT(I).EQ.0) THEN
                     IPT(NCENT) = I
                     GO TO 50
                  END IF
 40            CONTINUE



 50            L = IPT(NCENT)
               NPIXCENT(L) = ILEN
               SUMXI(L) = SUMOD*FXVAL
               SUMYI(L) = SUMYOD
               SUMI(L) = SUMOD
               IYBEG(L) = IBL
               IYEND(L) = IEL
               IXCENT(L) = NINT(SUMXI(L)/SUMI(L))
               IYCENT(L) = NINT(SUMYI(L)/SUMI(L))
               IXMAX(L) = IXRAST

               X_ORIGIN = IXRAST - IXCEN
               Y_ORIGIN = IBL - IYCEN

               AZIMUTHAL_ANGLE(L) = -ATAN2(Y_ORIGIN,X_ORIGIN)
C
C====== Initialise min/max rotated variables
C
               SPOT_LENGTH_MIN(L) = 10000.0
               SPOT_LENGTH_MAX(L) = -10000.0

               SPOT_WIDTH_MIN(L) = 10000.0
               SPOT_WIDTH_MAX(L) = -10000.0
C
C====== Rotate pixels in strip by azimuthal angle
C
               DO 60 IY = IBL , IEL
                  X_ORIGIN = IXRAST - IXCEN
                  Y_ORIGIN = IY - IYCEN
                  CALL LFN_ROTATE_XY(AZIMUTHAL_ANGLE(L),X_ORIGIN,
     +                           Y_ORIGIN,X_ROT,Y_ROT)

                  SPOT_LENGTH_MIN(L) = MIN(X_ROT,SPOT_LENGTH_MIN(L))
                  SPOT_LENGTH_MAX(L) = MAX(X_ROT,SPOT_LENGTH_MAX(L))

                  SPOT_WIDTH_MIN(L) = MIN(Y_ROT,SPOT_WIDTH_MIN(L))
                  SPOT_WIDTH_MAX(L) = MAX(Y_ROT,SPOT_WIDTH_MAX(L))

 60            CONTINUE
C
C====== Flag spots if srip overlaps search area
C
               IF (IBL.EQ.IYL .OR. IXRAST.EQ.IXL) THEN
                  IFLAG(L) = 1
               ELSE
                  IFLAG(L) = 0
               END IF
C
C====== Ongoing line
C
            ELSE
               SUMYOD = SUMYOD + REAL(IYPT*ITHRESH)
               SUMOD = SUMOD + REAL(ITHRESH)

            END IF
 70      CONTINUE
C
C====== End of loop over pixels
C
C====== Loop over active spots and store those that have finished
C
C====== Skip if no active spots
C
         IF (NCENT.EQ.0) GO TO 120
         ICENT = 0
 80      ICENT = ICENT + 1
         J = IPT(ICENT)
C
C====== Spots considered finished if no pixels for spot found on last stripe
C
         IF (IXRAST.EQ.IXMAX(J)) GO TO 110
C
C====== Reject spot if it overlaps limits of search area
C
         IF (IFLAG(J).EQ.1) GO TO 90
C
C====== Spots coodrinates as centre of gravity
C
         SX = SUMXI(J)/SUMI(J)
         SY = SUMYI(J)/SUMI(J)
C
C====== Store spots if centroid lies within RMM of search box centre
C
         SPOT_WIDTH = ABS(SPOT_WIDTH_MIN(J)-SPOT_WIDTH_MAX(J)) + 1.0
         SPOT_LENGTH = ABS(SPOT_LENGTH_MIN(J)-SPOT_LENGTH_MAX(J)) + 1.0
         R1 = SQRT((IXCEN-SX)**2+(IYCEN-SY)**2)
         IF (R1.GT.IRMAX) GO TO 90
         NSPOT = NSPOT + 1
         SUMINT = SUMI(J)
         NPIX = NPIXCENT(J)
         SPOT_ANGLE = ATAN2(SY-IYCEN,SX-IXCEN)*DEG
         CALL DMLF_ADDREC(MINDX_SPOT,IBUF,IERR)
         IF (IERR.NE.0) THEN
            IERR = 1
            ERRSTR = 'Memory allocation error when adding spot to list'
            GO TO 140
         END IF

C
C====== Spot stored or rejected
C====== Set flag and shuffle IPT down to fill gap in array
C
 90      NPIXCENT(J) = 0
         NCENT = NCENT - 1
C
C====== Skip if last spot
C
         IF (ICENT.GT.NCENT) GO TO 120


         DO 100 I = ICENT , NCENT
            IPT(I) = IPT(I+1)
 100     CONTINUE


 110     IF (ICENT.LT.NCENT) GO TO 80
 120     IF (CANCEL(0)) THEN
            IERR = 100
            ERRSTR = 'Find spots routine interrupted'
            GO TO 140
         END IF
         ISTRIP = IXRAST-IXL+1
         CALL PROGRESS(1,NSTRIPS,ISTRIP)
 130  CONTINUE
 140  CALL PROGRESS(2,NSTRIPS,ISTRIP)
      RETURN
      END
