C+
C NAME:
C	smei_frm_measle
C PURPOSE:
C	Measle detector
C CATEGORY:
C CALLING SEQUENCE:
	subroutine smei_frm_measle(kill,icam,nx,ny,frame,count,count_big,mask,indx)
C INPUTS:
C	kill		logical		.TRUE. : set measles to bad value
C					.FALSE.: set measles to neighbours average
C	icam		integer		camera number (1,2,3)
C					not used yet
C	nx		integer		horizontal frame size (1272 for mode 0)
C	ny		integer		vertical frame size (256 for mode 0)
C	frame(nx,ny)	real		SMEI frame after pedestal and dark current have been subtracted
C	indx(nx,ny)	integer		scratch array
C OUTPUTS:
C	frame(nx,ny)	real		SMEI frame with value of measles replace by a
C					an neighbours average.
C					Only the uncovered pixels are tested. Pedestal and dark current
C					areas are not modified.
C	mask(nx,ny)	integer		1 if pixel is OK
C					0 if pixel is flagged as measle
C					mask will be 1 for all pedestal and dark current pixels.
C	count		integer		# measles found
C	count_big	integer		# big measles (with at most 3 of the neighbour pixels
C					modified too
C INCLUDE:
	include		'smei_frm_layout.h'
C CALLS:
C	IndexR4, iArrR4ValuePresent, ArrR4GetMinMax, ArrR4Constant
C	smei_FitPlane, BadR4, smei_frm_mode, Say
C PROCEDURE:
C	Uses the ADUs departure of
C	a pixel from its neighbors-average z0 to trigger its replacement/elimination.
C	The response must be above a fixed value "val_cut" and be at least double the
C	8-pixel average response z0. Also, the departure from the neighbors average
C	is compared to the difference of neighboring pixels on opposing sides of the
C	center pixel.
C MODIFICATION HISTORY:
C	JUN-2004, Andy Buffington (UCSD/CASS; abuffington@ucsd.edu), Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    logical	kill
	    integer	icam
	    integer	nx
	    integer	ny
	    real	frame(nx,ny)
	    integer	count
	    integer	count_big
	    integer	mask (nx,ny)
	    integer	indx (*)

	character	cSay*11				/'FindMeasles'/
	integer		smei_frm_mode

	logical		bfix(-1:1,-1:1)
	logical		bfixpix
	real		z   (-1:1,-1:1)

	! Mode 0 numbers are used by cal_pattern (closed shutter calibration pattern)
	! Mode 1 is used by frm_pattern for camera 3 (usually run in 2x2 binning, mode 1)
	! Mode 2 is used by frm_pattern for cameras 1 and 2 (usually in 4x4 binning, mode 2)

	real		val_cut_mode(0:2)		/27.76,55.52, 13.88/	!55.52 is a 1000-electron cut: we have 500 here.
	real		dev_cut_mode(0:2)		/25.00,25.00,400.00/	!5 standard deviation cut, since this is sigma**2

	bad  = BadR4()
	mode = smei_frm_mode(nx, nbin)

	!if (icam .eq. 3) then
	!    if (mode .eq. 2) call Say(cSay,'E','mode', 'wrong mode for camera 3')
	!else
	!    if (mode .eq. 1) call Say(cSay,'E','mode', 'wrong mode for camera 1 or 2')
	!end if

	! Sort frame
	! If there are any bad values in frame these end up at the bottom of the
	! sorted array (BadR4 is a very big negative value).

	ntot = nx*ny			! Total # pixels
	call IndexR4(1,ntot,1,ntot,frame,indx)

	val_cut = val_cut_mode(mode)
	dev_cut = dev_cut_mode(mode)

	ileft  = (SMEI__FRM_LEFT -1)/nbin+1	! Column with partially covered pixels
	iright =  SMEI__FRM_RIGHT   /nbin	! Column with partially covered pixels
	jtop   =  SMEI__FRM_NY      /nbin

	count     = 0
	count_big = 0
	call ArrI4Constant(ntot,1,mask)

	! Process frame from the maximum down, i.e. indx is accessed in reverse order.

	do n=ntot,1,-1

	    ip = indx(n)-1		! Zero-based index into frame

	    i  = mod(ip,nx)		! Zero-based row
	    j  = ip/nx			! Zero-based column

	    i  = i+1			! One-based row
	    j  = j+1			! One-based column

	    ! Loop over all pixels in the open area of the CCD.
	    ! Stay away from the edges (we are looking at neighbours)

	    if (i .gt. ileft .and. i .lt. iright .and. j .gt. 1 .and. j .lt. jtop) then

		val = frame(i,j)

		! Quit if any one pixel is below threshold

		!if (val .ne. bad .and. val .ne. -10.0 .and. val .lt. val_cut) return
		if (val .ne. bad .and. val .lt. val_cut) return

		! Copy 3x3 neighbourhood of i,j into z

		do jj=-1,1
		    do ii=-1,1
			z(ii,jj) = frame(i+ii,j+jj)
		    end do
		end do

		! Bail if 3x3 has any empty pixel

		if (iArrR4ValuePresent(9,z,Bad) .eq. 0) then

		    call smei_FitPlane(z,z,z0,rcx,rcy,sumzz)
		    call ArrR4GetMinMax(9,z,zmin,zmax)

		    ! z now contains residuals (after subtracting planar fit).
		    ! Bail if not above threshold

		    if (z(0,0) .ge. max(val_cut,zmax) .and. z(0,0)*z(0,0) .ge. dev_cut*sumzz) then

			! Center pixel is too high: count as measle

			count = count+1				! Increase measle counter

			if (kill) z0 = bad

			z0fit = z0				! Planar fit at center pixel
			nfix = 0				! Counter for pixels that need fixing

			do jj=-1,1
			    do ii=-1,1

				val = z(ii,jj)			! Center pixel always gets fixed
				bfixpix = ii .eq. 0 .and. jj .eq. 0

				if (.not. bfixpix) then		! Not the center pixel
				    bfixpix = val .gt. val_cut
								! Take away neighbor's contribution
				    if (bfixpix) then		! Neighbour is too high
					nfix = nfix+1
					if (z0 .ne. bad) z0 = z0-val/8.0
				    end if
				end if

				bfix(ii,jj) = bfixpix

			    end do
			end do

			if (nfix .eq. 0) then			! Only center needs fixing

			    frame(i,j) = z0
			    mask (i,j) = 0

			else if (nfix .le. 3) then		! Center, plus 1 to 3 neighbours need fixing

			    count_big = count_big+1

			    do jj=-1,1
				do ii=-1,1
				    if (bfix(ii,jj)) then
					frame(i+ii,j+jj) = z0 	! Put in average of low neighbours
					mask (i+ii,j+jj) = 0
				    end if
				end do
			    end do

			else					! Too many neighbors to fix
								! only fix center pixel
			    frame(i,j) = z0fit
			    mask (i,j) = 0

			end if

		    end if

		end if

	    end if

	end do

	return
	end