Anisotropic smoothing operations

[Ben Tupper <btupper(at)bigelow.org>]

Hello,

How can I smooth a 2d image with a rectangular (rather than square)
smoothing window?  For example, how can I smooth an image with a moving
boxcar that is 9 columns by 5 rows in size.  The SMOOTH function accepts
a only single argument to define the width (and height in the 2d case)
of the smoothing window.

The CONVOL function seems to be the direction I want to go because I can
define arbitrary dimensions for the smoothing window.  Something like
the following:

	window = REPLICATE(1.0, 9,5)
	smoothed = CONVOL(original, window, TOTAL(window), /Edge_Truncate)

This works fine, unless the original data has any NANs.  SMOOTH and
MEDIAN permit NAN checking, but CONVOL does not.

Here's a routine with an example - the routine creates two windows, one
for NaN-less image smoothing and one for Nan-full image smoothing.  Each
window shows 4 images, from left to right they are: original, smoothed
with a [5,40] window using CONVOL, smoothed with a [40,5] window using
CONVOL, and finally, smoothed using SMOOTH and a square window.

Is there a way to make CONVOL treat missing data a SMOOTH does?  

Thanks,

Ben


;-----STARTHERE

PRO anisotropic

device, get_decomposed = decomp
tvlct, red, green,blue, /get

device, decomposed = 0
loadct, 2

n = 200
image = dist(n)
k1 = replicate(1.0, 5,40)
k2 = replicate(1.0, 40,5)
a = CONVOL(image, k1, total(k1), /Edge_Truncate)
b = CONVOL(image, k2, total(k2), /Edge_Truncate)
smoothed = Smooth(Image, 10, /Edge_truncate, /NAN)


Window, /free, title = 'Smoothed without NANs', xs = n*4+30, ys = n
tvscl, image
tvscl, a, n+10,0
tvscl, b, 2*n+20,0
tvscl, smoothed, 3*n +30,0

	
	;change some values to NAN
image[25:50, 25:50] = !Values.F_NAN
a = CONVOL(image, k1, total(k1), /Edge_Truncate)
b = CONVOL(image, k2, total(k2), /Edge_Truncate)
smoothed = Smooth(Image, 10, /Edge_truncate, /NAN)

!P.Multi = [0,2,2]
Window,/Free, title = 'Smoothed with NANs', xs = n*4+30, ys = n

index = where( FINITE(image) EQ 0,count)
if count gt 0 then image[index] = min(image,/nan)

index = where( FINITE(a) EQ 0,count)
if count gt 0 then a[index] = min(a,/nan)

index = where( FINITE(b) EQ 0,count)
if count gt 0 then b[index] = min(b,/nan)

index = where( FINITE(smoothed) EQ 0, count)
If Count gt 0 then smoothed[index] = min(smoothed,/nan)

tvscl, image
tvscl, a, n+10,0
tvscl, b, 2*n+20,0
tvscl, smoothed, 3*n +30,0

device, decomposed = decomp
tvlct,red,green,blue


END
;-----ENDHERE

-- 
Ben Tupper
Bigelow Laboratory for Ocean Sciences
180 McKown Point Rd.
W. Boothbay Harbor, ME 04575
btupper@bigelow.org