; @ Chris Jacobsen, Stony Brook 
;+ 
;NAME: 
;	EDGEGAUSS 
; 
;LAST CHANGED: ----------------------------------- 25-feb-00 
; 
;PURPOSE: 
;	This procedure uses gaussian smoothing to edge enhance images 
; as part of the stack_align procedure 
; 
;CATEGORY: 
;	STAND ALONE: image processing 
; 
;CALLING SEQUENCE: 
; for procedures: 
;	EDGEGAUSS, image_array, sigma, dc_value, $ 
;    zero_edge=zero_edge, help=help 
; 
;CALLED FROM AXIS: 
;	Stack_align 
; 
;INPUTS: 
;	IMAGE_ARRAY - image 
;	SIGMA - width of Gaussian smooth 
;	DC_VALUE - offset to be removed prior to smoothing 
; 
;KEYWORDS: 
;	ZERO_EDGE - if set, force edges of image to zero 
;	HELP - print information message 
; 
;OUTPUTS: 
;	The edge-sharpened image is returned as image_array 
; 
;COMMON BLOCKS: none 
; 
;PROCEDURE: 
; Revised June 30, 1998 to let one zero out the trace along the 
; gaussian border.  If you set /zero_edge, the DC value from 
; 1 sigma to 2 sigma from the border is calculated and subtracted 
; from the image, and it is also returned to you 
; 
;MODIFICATION HISTORY: 
; (30-Jun-98 cjj) add zero_edge keyword 
; (25-feb-00 aph) AXIS standard documentation 
;- 
 
PRO edgegauss, image_array, sigma, dc_value, $ 
    zero_edge=zero_edge, help=help 
 
IF (n_params() eq 0) OR keyword_set(help) THEN BEGIN 
    print,'edgegauss, image_array, sigma' 
    return 
ENDIF 
 
svec = size(image_array) 
IF (svec(0) NE 2) THEN BEGIN 
    print,'image_array must be 2D' 
    return 
ENDIF 
nx = svec(1) 
ny = svec(2) 
 
n_sigma = -alog(1.e-6) 
n_roll = fix(1+sigma*n_sigma)>2 
exparg = findgen(n_roll)/float(sigma) 
rolloff = 1.-exp(-0.5*exparg*exparg) 
 
IF keyword_set(zero_edge) THEN BEGIN 
    edge_total = 0. 
    edge_count = long(0) 
 
    ; Bottom edge 
    y1 = sigma 
    y2 = y1+sigma-1 
    edge_total = edge_total+$ 
      total(image_array(0:(nx-1),y1:y2)) 
 
    ; Top edge 
    y3 = ny-2*sigma 
    y4 = y3+sigma-1 
    edge_total = edge_total+$ 
      total(image_array(0:(nx-1),y3:y4)) 
 
    ; Left edge 
    x1 = sigma 
    x2 = x1+sigma-1 
    edge_total = edge_total+$ 
      total(image_array(x1:x2,y2:y3)) 
 
    ; Right edge 
    x3 = nx-2*sigma 
    x4 = x3+sigma-1 
    edge_total = edge_total+$ 
      total(image_array(x3:x4,y2:y3)) 
 
    n_pixels = $ 
        long(y2-y1+1)*long(nx)+long(y4-y3+1)*long(nx)+$ 
        long(x2-x1+1)*long(y3-y2+1)+long(x4-x3+1)*long(y3-y2+1) 
    dc_value = edge_total/float(n_pixels) 
    image_array = temporary(image_array) - dc_value 
ENDIF 
 
; Bottom edge 
xstart = 0 
xstop = nx-1 
iy = 0 
FOR i_roll = 0, (n_roll-1) DO BEGIN 
    image_array(xstart:xstop,iy) = $ 
      image_array(xstart:xstop,iy) * rolloff(iy) 
    xstart = (xstart+1)<(nx/2-1) 
    xstop = (xstop-1)>(nx/2) 
    iy = (iy+1)<(ny-1) 
ENDFOR 
 
; Top edge 
xstart = 0 
xstop = nx-1 
iy = ny-1 
FOR i_roll = 0, (n_roll-1) DO BEGIN 
    image_array(xstart:xstop,iy) = $ 
      image_array(xstart:xstop,iy) * rolloff(ny-1-iy) 
    xstart = (xstart+1)<(nx/2-1) 
    xstop = (xstop-1)>(nx/2) 
    iy = (iy-1)>0 
ENDFOR 
 
; Left edge 
ystart = 1 
ystop = ny-2 
ix = 0 
FOR i_roll = 0, (n_roll-1) DO BEGIN 
    image_array(ix,ystart:ystop) = $ 
      image_array(ix,ystart:ystop) * rolloff(ix) 
    ystart = (ystart+1)<(ny/2-1) 
    ystop = (ystop-1)>(ny/2) 
    ix = (ix+1)<(nx-1) 
ENDFOR 
 
; Right edge 
ystart = 1 
ystop = ny-2 
ix = nx-1 
FOR i_roll = 0, (n_roll-1) DO BEGIN 
    image_array(ix,ystart:ystop) = $ 
      image_array(ix,ystart:ystop) * rolloff(nx-1-ix) 
    ystart = (ystart+1)<(ny/2-1) 
    ystop = (ystop-1)>(ny/2) 
    ix = (ix-1)>0 
ENDFOR 
 
return 
end 
 
