;+
; PURPOSE:
;	Image editing: interpolation of image pixels, the coordinates
;	are given interactively. Set minimum or maximum values.
; CALLING SEQUENCE:
;	Edit, image [, xAxis, yAxis, [ xMin, xMax, yMin, yMax ] ]
; INPUTS:
;	image: 2D array
;	xAxis, yAxis: the corresponding axes
;	xMin, xMax, yMin, yMax: boundaries of the part displayed
;		of the image, in pixels.
; OUTPUT:
;	image: the edited image
; KEYWORD:
; 	status: 'Done' or 'Not Done'
; RESTRICTION:
;	Works only in interactive mode.
;	Works only in association with enrad
;
; MODIFICATION HISTORY:
;	Created in December 1991 by A. Csillaghy
;		Inst. of Astronomy, ETH Zurich
;	Interpolation modified in Jan 93 for non-regualr axes.
;	Addition of "Eliminate Pixel" option, Jan 93, A.Cs.
;	Addition of "DESPIKE" option: PSH, 2004/11/25 
;-

FUNCTION CoordToPix, coord, axis, n

  IF axis(0) GE axis( n-1 ) THEN $
	RETURN, Max( Where( axis GE coord )) $
  ELSE RETURN, Min( Where( axis GE coord ))

END ; Coord to pixel


FUNCTION InterpolZone, image, xMin, xMax, yMin, yMax, $
	 direction

  xDir = direction EQ 'X'
  yDir = direction EQ 'Y'
  IF xDir THEN im = image ELSE im = Transpose(image)
    
  nx = N_Elements( im(*,0) )
  ny = N_Elements( im(0,*) )
  left = (xMin*xDir + yMin*yDir) > 0
  right = (xMax*xDir + yMax*yDir) < (nx-1)
  low =(yMin*xDir + xMin*yDir) 
  high = (yMax*xDir + xMax*yDir) 
 
 low = low + 1

 WHILE low LT high DO BEGIN
    im( left:right, low ) = ReplLinInt( im( left:right, low ) )
    low = low + 1
  ENDWHILE

  IF direction EQ 'X' THEN RETURN, im $
  ELSE RETURN, Transpose( im )

END ; interpolate zone
  

PRO Edit, image, xAxis, yAxis, xMin, xMax, yMin, yMax, nx, ny, $
	status, oper

  status = 'Not Done'
  
  IF N_Params() LT 1 THEN BEGIN
    Message, 'Usage: Edit,image [,xAxis,yAxis ' + $
	'[,xMin,xMax,yMin,yMax ] ]', /INFO, /CONT
    RETURN
  ENDIF
  sizIm = Size( image )

  IF sizIm(0) EQ 0 THEN BEGIN
    Message, 'Image is not defined ... ', /INFO, /CONT
    RETURN
  END

  len = 16
  eMenu = StrArr( len )

  eMenu(0) = 'Interpolation in Y direction'
  eMenu(1) = 'Interactive definition of the interpolation zone'
  eMenu(2) = 'Point interpolation ... '
  eMenu(3) = 'Region interpolation ... '
  eMenu(4) = 'Interpolation below and/or above tresholds'
  eMenu(5) = 'Set minimum value ... '
  eMenu(6) = 'Set maximum value ... '
  emenu(7) = 'Set constant value in region...'
  eMenu(8) = 'Cut least value '
  eMenu(9) = 'Cut highest value '
  eMenu(10) = 'Eliminate wrong channel/scan (Y, resp X value)'
  eMenu(11) = 'Eliminate automatically wrong channel(s)'
  eMenu(12) = 'Despike (Light)'
  eMenu(13) = 'Despike (Strong)'
  eMenu(14) = 'Exit this menu'
  eMenu(15) = 'Redisplay this menu'

  IF N_Elements(nx) EQ 0 THEN nx = N_Elements( xAxis )
  IF N_Elements(ny) EQ 0 THEN ny = N_Elements( yAxis )
  IF N_Elements(xAxis) EQ 0 THEN xAxis = LIndGen( nx )
  IF N_Elements( yAxis ) EQ 0 THEN yAxis = LIndGen( ny )
  IF N_Elements( xMin ) EQ 0 THEN BEGIN
    xMin = 0L & xMax = nx-1
    yMin = 0L & yMax = ny-1
  ENDIF

  IF nx EQ 1 THEN stepPosX = 1 $
  ELSE stepPosX = (xAxis(1)-xAxis(0)) GE 0 
  IF ny EQ 1 THEN BEGIN
     stepPosY = 1 
     direction = 'X'
     xDir = 1 & yDir = 0
    eMenu(0) = 'Interpolation in X Direction'  
  ENDIF ELSE  BEGIN
    stepPosY = (yAxis(1)-yAxis(0)) GE 0 
    direction = 'Y' & yDir =1 & xDir = 0
  ENDELSE

  mode = 'Interactive'
  oper = ''

  REPEAT BEGIN

    choice = General_Menu( eMenu, 'Choose one Option or Exit', $
	OPAQUE = [ 0,0,0,0,0,0,1 ] )

    CASE choice OF
      1: BEGIN
            IF direction EQ 'Y'  THEN BEGIN
              direction = 'X'  & xDir = 1 & yDir = 0
            ENDIF ELSE BEGIN
              direction = 'Y'  & xDir = 0 & yDir = 1
            ENDELSE
            eMenu(0) = 'Interpolation in ' + direction + ' Direction'
          END
      2: BEGIN
            IF mode EQ 'Interactive' THEN mode = 'Static' $
            ELSE mode = 'Interactive'
            eMenu(1) = mode + ' Definition of the Interpolation Zone'
          END
      3: BEGIN
            IF mode EQ 'Static' THEN BEGIN
                Print, ''
                Print, 'Interpolate Pixel'
                Print, 'What is the coordinate of the pixel to interpolate ? '
                Read_Test, '     X Coordinate', Min( xAxis ), xCoord, $
                  HMS = GetHMS()
                Read_Test, '     Y Coordinate', Min( yAxis ), yCoord
                xPix = CoordToPix( xCoord, xAxis, nx )
                yPix = CoordToPix( yCoord, yAxis, ny )
            ENDIF ELSE BEGIN
                Pixel_Read, image( xMin: xMax, yMin: yMax ), $
                  xAxis( xMin: xMax ), yAxis( yMin: yMax ), dummy, $
                  xCoord, yCoord, /SIZES_ON, $
                  TEXT = 'Click on the pixel to interpolate', $
                  PIXELX = xPix, PIXELY = yPix
                xPix = xPix + xMin & yPix = yPix + yMin
            ENDELSE

            image = InterpolZone( image, xPix-1, xPix+1, $
                                  yPix-1, yPix+ 1, direction ) 
            status = 'Done' & oper = 'i'
        END                     ; point interpolation
        4: BEGIN
            IF mode EQ 'Static' THEN BEGIN
                Print, ''
                Print, 'Interpolate Region'
                Print, $
                  'What is the coordinate of the lower left pixel to interpol ? ' 
                Read_Test, '     X Coordinate', Min( xAxis ), xCoord1
                Read_Test, '     Y Coordinate', Min( yAxis ), yCoord1
                Print, $
                  'What is the coordinate of the upper right pixel to interpol ? ' 
                Read_Test, '     X Coordinate', Max( xAxis ), xCoord2
                Read_Test, '     Y Coordinate', Max( yAxis ), yCoord2
                xPix1 = CoordToPix( xCoord1, xAxis, nx )
                xPix2 = CoordToPix( xCoord2, xAxis, nx )
                yPix1 = CoordToPix( yCoord1, yAxis, ny )
                yPix2 = CoordToPix( yCoord2, yAxis, ny )
            ENDIF ELSE BEGIN
                yPix1 = 1 & yPix2 = 1 & xPix1 = 1 & xPix2 = 1
                DefineZone, image( xMin: xMax, yMin: yMax ), $
                  xAxis( xMin: xMax ), yAxis( yMin: yMax ), $
                  XPIXELS = xValues, YPIXELS = yValues, $
                  /SIZES_ON
                xPix1 = xValues(0) + xMin
                yPix1 = yValues(0) + yMin
                xPix2 = xValues(1) + xMin
                yPix2 = yValues(1) + yMin
            ENDELSE
            IF xPix2 LT xPix1 THEN BEGIN
                tmp = xPix1 & xPix1 = xPix2  & xPix2 = tmp
            ENDIF
            IF yPix2 LT yPix1 THEN BEGIN
              tmp = yPix1 & yPix1 = yPix2 & yPix2 =  tmp
            ENDIF
            image = InterpolZone( image, xPix1-1, xPix2+1, $
	yPix1-1, yPix2+1, direction ) 
            status = 'Done' & oper = 'i'
          END
          5: BEGIN
              Print, ''
              Print,'Interpolation Inside Tresholds'
              Read_Test,'Minimum value ', $
                StrTrim(MIN(image),2),minimum
              Read_Test,'Maximum value ', $
                StrTrim(Max(image),2),maximum
              image = TreshInterp( image, minimum, maximum, direction )
              status = 'Done'  & oper = 'int'
          END
          6: BEGIN
              Print,'Set Minimum Value'
              Read_Test,'New minimum: ', $
                StrTrim(MIN(image),2),minimum
              image = image > minimum
              status = 'Done'  & oper = 'smi'
          END
          7: BEGIN
              Print, ''
              Print,'Set Maximum Value'
              Read_Test,'New maximum: ', StrTrim(MAX(image),2), $
                maximum
              image = image < maximum 
              status = 'Done'  & oper = 'sma'
          END
          8: BEGIN
              yPix1 = 1 & yPix2 = 1 & xPix1 = 1 & xPix2 = 1
              DefineZone, image( xMin: xMax, yMin: yMax ), $
                xAxis( xMin: xMax ), yAxis( yMin: yMax ), $
                XPIXELS = xValues, YPIXELS = yValues, $
                /SIZES_ON
              xPix1 = xValues(0) + xMin
              yPix1 = yValues(0) + yMin
              xPix2 = xValues(1) + xMin
              yPix2 = yValues(1) + yMin
              IF xPix2 LT xPix1 THEN BEGIN
                  tmp = xPix1 & xPix1 = xPix2  & xPix2 = tmp
              ENDIF
              IF yPix2 LT yPix1 THEN BEGIN
                  tmp = yPix1 & yPix1 = yPix2 & yPix2 =  tmp
              ENDIF
              Read_Test,'New maximum: What is the constant value to set: ', $
                StrTrim(AVG(image(xPix1:xPix2,yPix1:yPix2)),2), $
                constVal
              image(xPix1:xPix2,yPix1:yPix2) = constVal 
              status = 'Done' & oper = 'cnst'
          END
     
      9: BEGIN
           image = CutMinMax( image, /MINIMUM )
           status = 'Done'   & oper = 'cmi'
         END
      10: BEGIN
            image = CutMinMax( image, /MAXIMUM )
            status = 'Done' & oper = 'cma'
          END
      11: BEGIN
               pixelX = 1 & pixelY = 1
               Pixel_Read, image( xMin: xMax, yMin: yMax ), $
	xAxis( xMin: xMax ), yAxis( yMin: yMax ), dummy, $
	PIXELX = pixelX, PIXELY = pixelY, /SIZES_ON, $
	TEXT = 'Click on the X or Y value to eliminate'
               pixelX = (pixelX + xMin) < (nx-2) > (1)
               pixelY = (pixelY + yMin) < (ny-2) > (1)
               IF direction EQ 'Y' THEN BEGIN
                 CASE  pixelY OF
	0: BEGIN
                          image = image(*,1:ny-1)
                          yAxis = yAxis(1: ny-1)
                        END
                    ny-1: BEGIN
                           image = image(*,0:ny-2)
                           yAxis = yAxis(0:ny-2)
                         END
                    ELSE: BEGIN
	    image = [ [image(*,0:pixelY-1)], $
		[image(*,pixelY+1:ny-1) ] ]
                        yAxis = [ yAxis(0:pixelY-1), yAxis(pixelY+1:ny-1) ]
                     END
                   ENDCASE
                   ny = ny-1
                   yMax = yMax < (ny-1)
               ENDIF ELSE BEGIN
                 CASE  pixelX OF
	0: BEGIN
                          image = image(1:nx-1,*)
                          xAxis = xAxis(1: nx-1,*)
                        END
                    nx-1: BEGIN
                           image = image(0:nx-2,*)
                           xAxis = xAxis(0:nx-2,*)
                         END
                    ELSE: BEGIN
	     image = [image(0:pixelX-1, *), $
		 image(pixelX+1:nx-1,*)]
                         xAxis =  [ xAxis(0:pixelX-1), xAxis(pixelX+1:nx-1) ]
                     END
                   ENDCASE
                   nx = nx-1
                   xMax = xMax < (nx-1)
               ENDELSE
               status = 'Done' & oper = 'elim'
             END
      12: BEGIN
              ElimWrongChannels, image, xAxis, yAxis
              nx = N_Elements( xAxis )
              ny = N_Elements( yAxis )
              xMax = xMax < (nx-1)  & yMax = yMax < (ny-1)
              status = 'Done' & oper = 'elima'
            END
      13: BEGIN
      	image=phoenix_despike(image)
      	;image=phoenix_despike(image,stat=6)
              status = 'Done' & oper = 'despikL'
      END
      14: BEGIN
      	image=phoenix_despike(image,stat=6)
              status = 'Done' & oper = 'despikS'
      END
      ELSE:
    ENDCASE

  END UNTIL (choice EQ len-1) OR (choice EQ 10)


END ; edit REgion user interface