C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**   LFN_GET_SPOTS      **
C**************************
C
CD-Routine: Get Refinement Nodals List - LFN_GET_SPOTS
C This routine is used to get a list of refinement nodal spots from the
C currently generated Laue Reflection List and determine their positions
C from the image. The spots found will be stored in the refinement
C observations list which will be automatically cleared on entry to
C the subroutine. The routine accesses a number of LDM parameters
C internally including image and spot related parameters and the
C spot search threshold. The background is normally calculated by taking a box
C three times the spot search box size and averaging the lowest 50%
C of the pixels.
CD-end
CD-Fortran:
      SUBROUTINE LFN_GET_SPOTS (NH, NK, NL, SPOTBOX, SORT, 
     +                          CANCEL, IMG, ITYPE, NF_OFF, 
     +                          MAXHALF, IDATA, IWORK, 
     +                          NSEARCH, NOBS, RMS, IERR, ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Get refinment nodals spots list
C
C Author:  John W. Campbell, September 1994
C
C Arguments:
C
      INTEGER NH, NK, NL
      LOGICAL SORT
      REAL SPOTBOX
      LOGICAL CANCEL
      EXTERNAL CANCEL
      INTEGER IMG(*)
      INTEGER ITYPE
      INTEGER NF_OFF
      INTEGER MAXHALF
      INTEGER IDATA(*)
      INTEGER IWORK(*)
      INTEGER NSEARCH
      INTEGER NOBS
      REAL RMS
      INTEGER IERR
      CHARACTER*(*) ERRSTR
C
CD-Parameters:
C NH          i  (R)  'h' nodal selection index (maximum absolute nodal 
C                     index for including spots in list).
C NK          i  (R)  'k' nodal selection index (maximum absolute nodal 
C                     index for including spots in list).
C NL          i  (R)  'l' nodal selection index (maximum absolute nodal 
C                     index for including spots in list).
C SPOTBOX     r  (R)  Spot box size for spot position determination
C                     (if 0.0 then use default box size if the 
C                     spot_width = 0.0 or radial masking option if
C                     the spot_width is > 0.0). If negative then then
C                     a box of size ABS(SPOTBOX) is used for both the
C                     spot box and background box e.g. for trying a large
C                     box when starting refinement from a rather poor match.
C SORT        l  (R)  Sort LRL so that spots in the observations list
C                     will be in an efficient order for accessing the
C                     in-memory image.
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 IMG(*)      i  (R)  Array containing the (possibly packed) image data
C                     for the current plate.
C ITYPE       i  (R)  Image data type in IMG =1 unsigned byte
C                                            =2 unsigned two byte
C                                            =3 signed (full) integer
C                                            =4 squashed i2 data
C NF_OFF      i  (R)  Fast offset between (slow) rasters of image (e.g. 
C                     1'st dimension of the 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                     as indicated by ITYPE.
C MAXHALF     i  (R)  Maximum half size in rasters allowed for background
C                     determining box which will be three times the spot 
C                     search box (each dimension) in size (unless a
C                     negative value was given for SPOTBOX (see above).
C IDATA       i  (*)  Work array to hold the spot data returned from the
C                     image. Must be dimensioned to at least
C                     (2*MAXHALF+1)*(2*MAXHALF+1)
C IWORK       i  (*)  Work array dimensioned to at least (2*MAXHALF+1).
C NSEARCH     i  (W)  No. of spot positions searched for.
C NOBS        i  (W)  Returns the number of spot position observations
C                     found and returned in in the refinement observations 
C                     list.
C RMS         r  (W)  Rms deviation in mm between predicted and observed
C                     spot positions.
C IERR        i  (W)  Error flag = 0   OK, spots found (but check IWARN 
C                                      also)
C                                = 1   No spots list generated or no 
C                                      spots in list
C                                = 2   Laue Reflection List needs to be
C                                      re-generated as LDM parameters
C                                      have been changed since it was 
C                                      created
C                                =-1   No. of 2'nd axis image rasters  
C                                      greater than 2*MAXHALF+1
C                                =-2   Required half box size greater
C                                      than MAXHALF
C                                =-3   Invalid image type requested
C                                =-4   Byte handling not available
C                                =-100 Cancelled via CANCEL routine
C ERRSTR      c  (W)  Error message string (max 80 chars) if IERR is
C                     non-zero
CD-end
C
C====== Parameter statements
C 

C
C====== GLOBALS:
C

C
C====== FUNCTION TYPES:
C

C
C====== LOCALS:
C
      LOGICAL OKNOD
      LOGICAL SPOV
      LOGICAL MTX
      LOGICAL GEN
      LOGICAL OVL
      LOGICAL USE_RADIAL
      INTEGER I
      INTEGER J
      INTEGER L
      INTEGER IR
      INTEGER IT
      INTEGER NBG
      INTEGER NBGPIX
      INTEGER IH
      INTEGER IK
      INTEGER IL
      INTEGER MULT
      INTEGER KERR
      INTEGER IORD
      INTEGER ISWAP
      INTEGER IFLAG
      INTEGER ND
      INTEGER ISPOT
      INTEGER NSPOTS
      INTEGER IPACK
      INTEGER IPLATE
      INTEGER NXRAST
      INTEGER NYRAST
      INTEGER IXSIZ
      INTEGER IYSIZ
      INTEGER KXSIZ
      INTEGER KYSIZ
      INTEGER IXC
      INTEGER IYC
      INTEGER ITHRESH
      INTEGER IBG
      INTEGER NOVER
      INTEGER N_H
      INTEGER N_K
      INTEGER N_L
      INTEGER MINHARM
      INTEGER MAXHARM
      INTEGER INCHARM
      REAL    XD
      REAL    YD
      REAL    ALAM
      REAL    DMTHR2
      REAL    ALNOD
      REAL    DTHR2NOD
      REAL    STHRESH
      REAL    SUMSQ
      REAL    X_CEN
      REAL    Y_CEN
      REAL    MM_RAST_X
      REAL    MM_RAST_Y
      REAL    SPOTW
      REAL    SPOTL
      REAL    SPOTB
      REAL    SPOTD
      REAL    SPOT_B
      REAL    SPOT_D
      REAL    TOO_FAR
      REAL    CGX
      REAL    CGY
      REAL    DIST2
      REAL    BGFACT
      REAL    EXPAND
      REAL    PERCNT
      CHARACTER*2 VALSTR
      CHARACTER*2 DUMSTR

C     LOGICAL OKNOD        !Non-spatially overlapping nodal flag
C     LOGICAL SPOV         !Spatial overlap flag
C     LOGICAL MTX          !Need to regenerate orientation matrix flag
C     LOGICAL GEN          !Need to re-predict LRL flag
C     LOCICAL OVL          !Need to re-calculate spatial overlaps flag
C     LOGICAL USE_RADIAL   !Use radial masking flag
C     INTEGER I            !Temp/loop variable
C     INTEGER J            !Temp/loop variable for sort
C     INTEGER L            !Temp/loop variable for sort
C     INTEGER IR           !Temp/loop variable for sort
C     INTEGER IT           !Temp/loop variable for sort
C     INTEGER NBG          !No. pixels in background box
C     INTEGER NBGPIX       !No. pixels to use as background
C     INTEGER IH           !h index
C     INTEGER IK           !k index
C     INTEGER IL           !l index
C     INTEGER MULT         !Multiplicity
C     INTEGER KERR         !Error flag
C     INTEGER IORD         !Axis order flag
C     INTEGER ISWAP        !Byte swap flag
C     INTEGER IFLAG        !LDM parameter status flag
C     INTEGER ND           !No. of decimal places fla
C     INTEGER ISPOT        !Spot index in spot list
C     INTEGER NSPOTS       !No. of spots in spot list
C     INTEGER IPACK        !Current pack number
C     INTEGER IPLATE       !Current plate number
C     INTEGER NXRAST       !No. x-rasters
C     INTEGER NYRAST       !No. y-raster
C     INTEGER IXSIZ        !Half box x size for finding spot
C     INTEGER IYSIZ        !Half box y size for finding spot
C     INTEGER KXSIZ        !Half box x size for finding background
C     INTEGER KYSIZ        !Half box y size for finding background
C     INTEGER IXC          !Nearest x pixel position of box centre
C     INTEGER IYC          !Nearest y pixel position of box centre
C     INTEGER ITHRESH      !Threshold above background for finding spots
C     INTEGER IBG          !Background value at spot position
C     INTEGER NOVER        !No. of pixels in box above threshold
C     INTEGER N_H          !Nodal 'h' index
C     INTEGER N_K          !Nodal 'k' index
C     INTEGER N_L          !Nodal 'l' index
C     INTEGER MINHARM      !Minimum harmonic
C     INTEGER MAXHARM      !Maximum harmonic
C     INTEGER INCHARM      !Harmonic increment
C     REAL    XD           !Distortion corrected 'xd' coordinate (rasters)
C     REAL    YD           !Distortion corrected 'yd' coordinate (rasters)
C     REAL    ALAM         !Wavelength
C     REAL    DMTHR2       !dmin_threshold2
C     REAL    ALNOD        !Wavelength for nodal 
C     REAL    DTHR2NOD     !dmin_threshold2 for nodal
C     REAL    STHRESH      !Spot search threshold
C     REAL    SUMSQ        !Sum of square of obs-calc position differences
C     REAL    X_CEN        !Corrected x centre in raster
C     REAL    Y_CEN        !Corrected y centre in rasters
C     REAL    MM_RAST_X    !mm to x-rasters conversion factor
C     REAL    MM_RAST_Y    !mm to y-rastrs conversion factor
C     REAL    SPOTW        !Spot width
C     REAL    SPOTL        !Spot length
C     REAL    SPOTB        !Spot border
C     REAL    SPOTD        !Spot delta
C     REAL    SPOT_B       !Spot box size to use
C     REAL    SPOT_D       !Spot diameter to use in closeness calculations
C     REAL    TOO_FAR      !Max shift poss. (mm) with all spot inside box
C     REAL    CGX          !Found spot c_of_g x position
C     REAL    CGY          !Found spot c_of_g y position
C     REAL    DIST2        !Square of obs-calc distance
C     REAL    BGFACT       !Background determining box is this factor * normal
C                          !spot box
C     REAL    EXPAND       !Expand factor for spot box on second try
C     REAL    PERCNT       !Percentage of pixels for background
C     CHARACTER*2 VALSTR   !LDM parameter value string (dummy)
C     CHARACTER*2 DUMSTR   !Dummy string
C
      DATA PERCNT /50.0/
C
C-------------------------------------------------------------------------------
C
C====== Initialise
C
      BGFACT = 3.0
      IF (SPOTBOX.LT.0.0) BGFACT = 1.0
      EXPAND = 1.5
      IF (SPOTBOX.LT.0.0) EXPAND = 1.0
      IERR = 0
      ERRSTR = ' '
      CALL LFN_CLROBS
      NOBS = 0
      NSEARCH = 0
      CALL LRL_NUMSPOTS (NSPOTS)
      IF (NSPOTS.LE.0) THEN
         IERR = 1
         ERRSTR = 'No spots list generated or no spots in list'
         RETURN
      ENDIF
      CALL LRL_GET_PLATE (IPACK, IPLATE)
      CALL LDM_CHK_GEN (IPACK, IPLATE, MTX, GEN, OVL)
      IF (GEN.OR.OVL) THEN
         IERR = 2
         ERRSTR = 'LRL needs to be regenerated with current parameters'
         RETURN
      END IF
      CALL LDM_GET_STHR (STHRESH, IFLAG, ND)
      SUMSQ = 0.0
      RMS = 0.0
      CALL LDM_MMRAST (IPACK, IPLATE, MM_RAST_X, MM_RAST_Y)
      CALL LRL_GET_XYCEN (X_CEN, Y_CEN)
      CALL LDM_GET_AXORD (VALSTR, IORD, ISWAP, IFLAG)
      CALL LDM_GET_NXRAS (NXRAST, IFLAG)
      CALL LDM_GET_NYRAS (NYRAST, IFLAG)
C
C====== Sort LRL if requested
C
      IF (SORT) THEN
         IF (IORD.LE.4) THEN
            CALL LRL_SORT (3)
         ELSE
            CALL LRL_SORT (4)
         END IF
      END IF
      IF (CANCEL(0)) THEN
         IERR = -100
         ERRSTR = 'Spot search cancelled'
         RETURN
      END IF
C
C====== Get spot and box parameters
C
      CALL LDM_GET_SPOTW (IPACK, IPLATE, SPOTW, IFLAG, ND, 
     +                    KERR, DUMSTR)
      USE_RADIAL = .FALSE.
      IF (SPOTBOX.EQ.0.0.AND.SPOTW.NE.0.0) USE_RADIAL = .TRUE.
C
C====== Set up box size and spot diameter for local use for normal or
C====== radially streaked spots
C
      IF (USE_RADIAL) THEN
         CALL LDM_GET_SPOTL (IPACK, IPLATE, SPOTL, IFLAG, ND, 
     +                       KERR, DUMSTR)
         CALL LDM_GET_SPOTB (IPACK, IPLATE, SPOTB, IFLAG, ND, 
     +                       KERR, DUMSTR)
         SPOT_B = SPOTL + (2.0 + SPOTB)/MM_RAST_X
         SPOT_D = SPOTW
      ELSE
         IF (SPOTBOX.EQ.0.0) THEN
            CALL LDM_GET_SPOTL (IPACK, IPLATE, SPOTL, IFLAG, ND, 
     +                          KERR, DUMSTR)
            CALL LDM_GET_SPOTD (IPACK, IPLATE, SPOTD, IFLAG, ND, 
     +                          KERR, DUMSTR)

            SPOT_B = 2.0*SPOTD - SPOTL
            IF (SPOT_B.LT.SPOTL) SPOT_B = SPOTL
            SPOT_D = SPOTL 
         ELSE
            SPOT_B = ABS(SPOTBOX)
            SPOT_D = SPOTL
         END IF
      END IF
C
C====== Get background related parameters
C
      KXSIZ = BGFACT*MM_RAST_X*SPOT_B/2.0
      KYSIZ = BGFACT*MM_RAST_Y*SPOT_B/2.0
      NBG = (2*KXSIZ+1)*(2*KYSIZ+1)
      NBGPIX = NBG*PERCNT/100.0
      IF (NBGPIX.LT.1) NBGPIX = 1
C
C====== Check that box size OK (note assumed that BGFACT>=1.5 i.e.
C====== maximum size of (expanded) spot search box)
C
      IF (KXSIZ.GT.MAXHALF.OR.KYSIZ.GT.MAXHALF) THEN
         IERR = -2
         ERRSTR = 'Program cannot handle required background box size'
         RETURN
      END IF
C
C====== Too_far (mm) is max. shift possible with all the spot inside box.
C
      TOO_FAR = (SPOT_B - SPOT_D)/2.0
C
C====== Get list of refinement spots
C
      DO 500 ISPOT = 1, NSPOTS
         IXSIZ = MM_RAST_X*SPOT_B/2.0
         IYSIZ = MM_RAST_Y*SPOT_B/2.0
         CALL LRL_GET (ISPOT, IH, IK, IL, XD, YD, ALAM, DMTHR2, MULT,
     +                 OKNOD, SPOV, KERR)
         IF (.NOT.OKNOD) GO TO 500
         IF (SPOV) GO TO 500
         CALL LRL_GET_HARM (ISPOT, IH, IK, IL, N_H, N_K, N_L, MULT,
     +                      MINHARM, MAXHARM, INCHARM, ALNOD, 
     +                      DTHR2NOD, KERR)
         IF (ABS(N_H).GT.NH) GO TO 500
         IF (ABS(N_K).GT.NK) GO TO 500
         IF (ABS(N_L).GT.NL) GO TO 500
         NSEARCH = NSEARCH + 1
         IXC = NINT(XD)
         IYC = NINT(YD)
C
C====== Get approximate background for use with spot position searches
C
         CALL LFN_IMG_RECT (IMG, ITYPE, IORD, NF_OFF, 
     +                      IXC-KXSIZ, IYC-KYSIZ,
     +                      2*KXSIZ+1, 2*KYSIZ+1, 
     +                      NXRAST, NYRAST,
     +                      0, IDATA, IWORK, KERR)
         IF (KERR.EQ.-1) THEN
            IERR = -3
            ERRSTR = 'Invalid image type requested'
            RETURN
         ELSE IF (KERR.EQ.-2) THEN
            IERR = -4
            ERRSTR = 'Byte handling not available'
            RETURN
         END IF
         IF (CANCEL(0)) THEN
            IERR = -100
            ERRSTR = 'Spot search cancelled'
            RETURN
         END IF
C
C====== Sort background pixels
C
         L = NBG/2 + 1
         IR = NBG
130      IF (L.GT.1) THEN
            L = L - 1
            IT = IDATA(L)
         ELSE
            IT = IDATA(IR)
            IDATA(IR) = IDATA(1)
            IR = IR - 1
            IF (IR.EQ.1) THEN
               IDATA(1) = IT
               GO TO 150
            END IF
         END IF
         I = L
         J = 2*L
140      IF (J.LE.IR) THEN
            IF (J.LT.IR) THEN
               IF (IDATA(J).LT.IDATA(J+1)) J = J + 1
            END IF
            IF (IT.LT.IDATA(J)) THEN
               IDATA(I) = IDATA(J)
               I = J
               J = J + J
            ELSE
               J = IR + 1
            END IF
            GO TO 140
         END IF
         IDATA(I) = IT
         GO TO 130
150      CONTINUE
         IBG = 0
         J = 0
         DO 160 I = 1, NBGPIX
            IF (IDATA(I).LE.0) GO TO 160
            IBG = IBG + IDATA(I)
            J = J + 1
160      CONTINUE
         IF (J.GT.0) IBG = IBG/J
C
C====== Search for spot
C
         ITHRESH = IBG + NINT(STHRESH)
         IF (USE_RADIAL) THEN
            CALL LFN_SPOT_RD (IPACK, IPLATE, IXC, IYC, 
     +                        IXSIZ, IYSIZ, IDATA, 
     +                        IWORK, ITHRESH, IMG, ITYPE, IORD, 
     +                        NF_OFF, NXRAST, NYRAST,
     +                        X_CEN, Y_CEN, NOVER, 
     +                        CGX, CGY, KERR)
         ELSE
            CALL LFN_SPOT_CG (IXC, IYC, IXSIZ, IYSIZ, IDATA, 
     +         IWORK, ITHRESH, IMG, ITYPE, IORD, NF_OFF,
     +         NXRAST, NYRAST, NOVER, 
     +         CGX, CGY, KERR)
         END IF
         IF (KERR.EQ.-1) THEN
            IERR = -3
            ERRSTR = 'Invalid image type requested'
            RETURN
         ELSE IF (KERR.EQ.-2) THEN
            IERR = -4
            ERRSTR = 'Byte handling not available'
            RETURN
         END IF
         IF (NOVER.GT.0) THEN
            DIST2 = (XD-CGX)**2 + (YD-CGY)**2
C
C====== Check each coord to see if some of spot is outside c of g box.
C
            IF (ABS(XD-CGX).GT.(TOO_FAR*MM_RAST_X)) THEN
C
C====== Because some of spot is outside box, expect true c of g to be
C====== further towards edge of box. An approximation of distance
C====== to move is the amount of spot that`s been cut off.
C====== ie. (cgx-xfd) + spot_diam/2 - box_size/2 (put into right units)
C

               IXC = CGX + (CGX-XD)
     +         + ((SPOT_D - SPOT_B)/2.*MM_RAST_X)
            ELSE
               IXC = CGX
            END IF
C
            IF (ABS(YD-CGY).GT.(TOO_FAR*MM_RAST_Y)) THEN
               IYC = CGY + (CGY-YD)
     +         + ((SPOT_D - SPOT_B)/2.*MM_RAST_Y)
            ELSE
               IYC = CGY
            END IF
C
C====== Re-centre if we didn`t get a true centre of gravity.
C
            IF (ABS(XD-CGX).GT.(TOO_FAR*MM_RAST_X) .OR. 
     +          ABS(YD-CGY).GT.(TOO_FAR*MM_RAST_Y)) THEN
C
C====== True c of g should be within a spot radius of here, but
C====== make box 1.5 times bigger just in case unless using fixed
C====== box.
C                 
               IXSIZ = SPOT_B/2.0*EXPAND*MM_RAST_X
               IYSIZ = SPOT_B/2.0*EXPAND*MM_RAST_Y
               IF (USE_RADIAL) THEN
                  CALL LFN_SPOT_RD (IPACK, IPLATE, IXC, IYC, 
     +                 IXSIZ, IYSIZ, 
     +                 IDATA, IWORK, ITHRESH, IMG, ITYPE,
     +                 IORD, NF_OFF, NXRAST, NYRAST,
     +                 X_CEN, Y_CEN, NOVER, 
     +                 CGX, CGY, KERR)
               ELSE
                  CALL LFN_SPOT_CG (IXC, IYC, IXSIZ, IYSIZ, 
     +                 IDATA, IWORK, ITHRESH, IMG, ITYPE, IORD, 
     +                 NF_OFF, NXRAST, NYRAST, NOVER, 
     +                 CGX, CGY, KERR)
               END IF
               DIST2 = (XD-CGX)**2 + (YD-CGY)**2
            END IF
            IF (NOVER.GT.0) THEN
               CALL LFN_ADDOBS (IH, IK, IL, CGX, CGY)
               SUMSQ = SUMSQ + DIST2
               NOBS = NOBS + 1
            END IF
         END IF
         IF (CANCEL(0)) THEN
            IERR = -100
            ERRSTR = 'Spot search cancelled'
            RETURN
         END IF
500   CONTINUE
      IF (NOBS.GT.0) RMS = SQRT(SUMSQ/NOBS)/MM_RAST_X
      RETURN
      END
