FUNCTION DENOISE, image, coeff, softthreshold
;+
; $Id: denoise.pro,v 1.1 2008/03/31 14:36:39 nathan Exp $
; NAME:
;	DENOISE
;
; PURPOSE:
;	Denoise an image using wavelet soft thresholding.
;
; CATEGORY:
;	IMAGE PROCESSING
;
; CALLING SEQUENCE:
;	Result = DENOISE(Image, Coeff, Softtheshold)
;
; INPUTS:
;	Image:	The image to be denoised.
;	Coeff:	May be 4, 12, or 20.  The different values specify
;		different wavelet bases.
;	Softthreshold:  Wavelet coefficients that are below an automatically
;		set threshold are scaled by this amount.  If set to 1, then
;		the filter does nothing.  The lower the value (to zero), the
;		stronger the denoising.
;
;
; OUTPUTS:
;	The denoised image. 
;
; PROCEDURE:
;	Uses WTN for the discrete wavelet transform.  WTN can only handle images
;	that have dimensions that are powers of 2, so the image is first 
;	congrided up to the closest power of 2 dimensions.  After the 
;	denoising, the image is congrided back to its original size.  The 
;	threshold used to determine which wavelet coefficients are 
;	insignificant is SIGMA * SQRT(2*LOG(N)), where SIGMA is the
;	stardard deviation of the image and N is the number of points in 
;	the image.
;
; EXAMPLE:
;	NewImage = DENOISE(Image, 4, .7)
;
; MODIFICATION HISTORY:
; 	Written by:	Jeff Walters, Aug 2001
;-
; $Log: denoise.pro,v $
; Revision 1.1  2008/03/31 14:36:39  nathan
; moved from adam/programs, was dated Nov 7, 2007
;



;*** A simple routine to apply wavelet coefficient 
;*** thresholding

;*** Image must have dims that are a power of 2
   s = size(image)
   nx = s(1)
   ny = s(2)
   
   newx = 2.^(ceil(alog(nx)/alog(2)))
   newy = 2.^(ceil(alog(ny)/alog(2)))
  
   newimage = CONGRID(image, newx, newy, cubic=-0.5)
   
   n = newx*newy
   d = stddev(newimage)
   
   thresh = d*SQRT(2*alog(n))
   
   trans = wtn(newimage, coeff)
   
   trans = trans - mean(trans)
   
   q = where(abs(trans) le thresh)
   
   trans(q) = softthreshold*trans(q)
   
   newimage = WTN(trans, coeff, /INVERSE)
   
   newimage = CONGRID(newimage, nx, ny, cubic = -0.5)
   
   RETURN, newimage
END