; $Id: dev_filter2.pro,v 1.1 2007/11/08 18:30:32 nathan Exp $
; NAME:
;	dev_filter2
;
; PURPOSE:
;	This function standardizes the contrast of an image
;	so that it is the same throughout.
;
; CATEGORY:
;	Image tool
;
; CALLING SEQUENCE:
;	result = DEV_FILTER(image, factor)
;
; INPUTS:
;	image:  The image that will be adjusted
;       factor: A constant that veries the strength of the filter.  Greater than one
;          causes low contrast areas to have a higher contrast than originally high contrast areas.
;          Less than one retains some of the relative contrast in the image.
;
;	
; KEYWORD PARAMETERS:
;	BOX:  Setting the keyword BOX causes the contrast to be standardized for rectangular
;		regions on the image.  The default is to standardize the contrast for each 
;		one pixel column.
;
; OUTPUTS:
;	This function returns the input image, but with the contrast standardized.
;
; RESTRICTIONS:
;	This procedure is specifically designed for height-time images one Carrington
;       rotation in length.  The ht-map should be 'tall,' not 'long,' i.e. as saved in the
;       data files in redwood (The dimensions should be something like [256, 5000]).
;
; PROCEDURE:
;	The pixel values are first centered around zero and then each box or column is divided
;       by its standard deviation.
;
; EXAMPLE:
;	result = DEV_FILTER(image, 1, /BOX)
;
; MODIFICATION HISTORY:
; 	Written by:	Jeff Walters, NRL Summer Student, July 1997
; $Log: dev_filter2.pro,v $
; Revision 1.1  2007/11/08 18:30:32  nathan
; moved from adam/jprogs
;
;
;	August 1997, 	Haleh Hadavand, NRL Summer Student, Added box method
;	Sept 1997, 	JHW, Added c2 compatibility, renamed dev_filter2
;       Aug 2001,       JWH, Added a target boxsize and minimized edge effects.	
;-

FUNCTION boxsize_error, size, box, target
   n = fix(size / box)
   e = (size - n*box)^2 + (box - target)^2
   boxsize_error = e
   return, boxsize_error
END

FUNCTION best_box_size, size, target, max_deviation
;*** This function returns the best box size to use
;*** for a given map length and a target box size.  The
;*** function minimizes an error metric that is defined
;*** in BOXSIZE_ERROR.  The error is a function of the 
;*** number of leftover pixles at the edge of a map and
;*** the distance from the preset target distance.

   n = 2*fix(max_deviation*2) ;make n even
   
   x = target - (n / 2) + indgen(n)
   
   e = fltarr(n)
   
   FOR i=0, n-1 DO BEGIN
      e(i) = boxsize_error(size, x(i), target)
   ENDFOR
   
   min = min(e, index)
   
   best_box_size = x(index)
   
   RETURN, best_box_size

END
  
FUNCTION dev_filter2, img, factor, BOX = box, C2 = c2
;***
;*** Modifiy the two constants below to change the box
;*** size.  

   SZ_RADIAL = 15
   SZ_TIME = 100
;___________________
   
   s = size(img)
   img = img - mean(img)
   cols = s(1)
   rows = s(2)
   nimg = img
   start = 30+38*KEYWORD_SET(c2)
   
   IF KEYWORD_SET(box) THEN BEGIN
      ;** 
      ;** Box Method
      ;** 
      
      ;** Calculate the box sizes.
      ;**
      ;** SZ_RADIAL and SZ_TIME will be adjusted slightly so that
      ;** a they are as close to the preset values as possible
      ;** but edge effects are minimized.  See best_box_size.pro
      ;** for details on how this is done
      
      stripx = best_box_size(cols-start, SZ_RADIAL, 1)
      stripy = best_box_size(rows, SZ_TIME, FIX(0.15*SZ_TIME))
      
      ;stripx=FIX(cols/10.0)
      ;stripy=FIX(rows/100.0)
      nnimg=dblarr(cols,rows)
   
      FOR i=start, cols - stripx,stripx DO BEGIN
         FOR j=0,rows-stripy,stripy DO BEGIN
            result = MOMENT(nimg(i:stripx-1+i, j:stripy-1+j))
            sd = result(1)
            sd =  (SQRT(sd))^(factor)
            IF (sd ge 0.001) THEN nnimg(i:stripx-1+i, j:stripy-1+j) = nimg(i:stripx-1+i, j:stripy-1+j) / (sd) ;- MEDIAN(nimg(i:stripx-1+i, j:stripy-1+j) / (sd))
         ENDFOR
      ENDFOR
   
      IF (i+1 LT cols) THEN BEGIN  
        FOR j=0,rows-stripy,stripy DO BEGIN
          result=MOMENT(nimg(i-stripx:cols-1,j:j+stripy-1))
          sd=result(1)
          sd=(SQRT(sd))^(factor)
          IF (sd ge 0.001) THEN nnimg(i-stripx:cols-1,j:j+stripy-1)=nimg(i-stripx:cols-1,j:j+stripy-1)/(sd) ;- median(nimg(i-stripx:cols-1,j:j+stripy-1)/(sd))
         ENDFOR
      ENDIF

      IF (j+1 LT rows) THEN BEGIN  
         FOR i=start, cols - stripx,stripx DO BEGIN
            result=MOMENT(nimg(i:i+stripx-1,j-stripy:rows-1))
            sd=result(1)
            sd=(SQRT(sd))^(factor)
            IF (sd ge 0.001) THEN nnimg(i:i+stripx-1,j-stripy:rows-1)=nimg(i:i+stripx-1,j-stripy:rows-1)/(sd) ;- median(nimg(i:i+stripx-1,j-stripy:rows-1)/(sd))
         ENDFOR
      ENDIF

      nnimg(0:start)=MEDIAN(nnimg)
      RETURN, nnimg
   ENDIF ELSE BEGIN
      ;** 
      ;** Column
      ;**
      FOR i=start, cols - 1 DO BEGIN
         result = MOMENT(nimg(i-1:i, *))
         sd = (SQRT(result(1)))^factor
         IF (sd ge 0.001) THEN nimg(i, *) = img(i, *) / (sd)
      ENDFOR
   ENDELSE
     nimg(0:start, *)=MEDIAN(nimg)
     help


   RETURN, nimg
END