' This version uses mxGaussianBlur as the blurrer (PSF)

Option Explicit

Dim Iters
Dim Radius
Dim StdDev

Const mxErosionKF = 9
Const mxUserKF = 10
Const mxGaussianBlur = 11

Dim mxapp	' the MaximDL application object
Dim ImageNum

''''''''''''
'Const imagename = "C:\d\sxv\05-06-26\m16HaRGB.tif"
'Dim imagename
'imagename=inputbox("Image name: ")
Iters=CInt(inputbox("How many iterations? "))
Radius=CDbl(inputbox("Enter FWHM (in pixels): "))
StdDev=CDbl(inputbox("Standard Deviation of Background: "))
''''''''''''''
Radius = Radius/2.35482

Dim image
Dim numerator
Dim denominator
Dim result
Dim tmp_result
Dim image_high_vals
Dim image_low_vals
Dim background


Set mxapp = CreateObject("Maxim.Application")

Set image = mxapp.CurrentDocument

Set background = image.Duplicate
background.DisplayName = "Background"
background.KernelFilter mxErosionKF, 7   ' 7x7 erosion
background.KernelFilter mxGaussianBlur, 5  ' Gaussian blur 3 pixel radius
background.AddConstant(5*StdDev)

Set image_high_vals = image.Duplicate
image_high_vals.DisplayName = "Image high_vals"
image_high_vals.Subtract(background)

Set image_low_vals = image_high_vals.Duplicate
image_low_vals.DisplayName = "Image low_vals"
image_low_vals.Add(background)
image_low_vals.Subtract(image)

Set result = image_high_vals.Duplicate
result.DisplayName = "Deconvolved " + image.DisplayName
Set tmp_result = image_high_vals.Duplicate
image.DoNotClose = TRUE
result.DoNotClose = TRUE

for ImageNum=1 to Iters
	Set numerator = image_high_vals.Duplicate
	Set denominator = result.Duplicate
	numerator.DisplayName = "Numerator"
	denominator.DisplayName = "Denominator"
	tmp_result.BringToTop
	denominator.KernelFilter mxGaussianBlur,  Radius 
	numerator.Divide(denominator)
	numerator.KernelFilter mxGaussianBlur,  Radius
	result.Multiply(numerator)
	numerator.Close
	denominator.Close
	tmp_result.Close
	Set tmp_result = result.Duplicate
	tmp_result.Add(background)
	tmp_result.Subtract(image_low_vals)
	tmp_result.DisplayName = "Richardson-Lucy iteration "+cstr(ImageNum)
next

tmp_result.Close
result.Add(background)
result.Subtract(image_low_vals)

image_low_vals.Close
image_high_vals.Close
background.Close

Dim flag
flag = result.SetFITSKey("HISTORY", "Richardson-Lucy/Vanderbei. Iters="+CStr(Iters)+", Radius="+CStr(2*Radius)+", StdDev="+CStr(StdDev))
'flag = result.setFITSKey("HISTORY","Richardson-Lucy (Vanderbei) "+CStr(Iters)+" iters, "+CStr(Radius)+" pixels, "+CStr(StdDev)+" std dev")