forward_function MSEdenoise

pro minMSEdenoise, stats=stats, sigStats=sigStats
common NoisyImage, in, noisy, sigma
stats = dblarr(8,10,8)
;OPENW, unit, 'stats.txt', /GET_LUN
;printf,unit, 'Sigma Noisy MSE Bayes Threshold Bayes MSE Opt Threshold ', $
;	'Opt MSE %Error Thresh %Error MSE'

images = ['pentagon.tiff', 'boat.tiff', 'testImg.tiff', 'girl.tiff', $
	'u2.tiff','stream.tiff', 'couple.tiff', 'airport.tiff']
sigStats = dblarr(8,2)

for j=0, 7 do begin
	in = read_tiff(images[j])
	in = in[0:511,0:511]
	sigStats[j,0] = mean(in)
	sigStats[j,1] = stddev(in)
	for i=0, 9 DO BEGIN
		print, i
		sigma = i*0.5+0.1
		stats[0,i,j] = sigma
		noise = randomn(1024,512,512,/double)*sigma
		noisy = in+noise
		stats[1,i,j] = total((noisy-in)^2)

		;; Get the Bayesian threshold and result
		clean = DIcal_denoise(noisy, thresh=thresh, sigma=sigma)
		MSE = total((in-clean)^2)
		stats[2,i,j] = 0;thresh
		stats[3,i,j] = MSE

		;; Minimize the MSE
		bestThresh=[10];[thresh]
		constrained_min, bestThresh,[0,1.0e30],[0,1.0e30], 0, $
			'MSEdenoise', stat, EPSTOP=1.0e-4, REPORT='err.txt'

		stats[4,i,j] = bestThresh
		stats[5,i,j] = MSEdenoise(bestThresh)
	
		stats[6,i,j] = (stats[2,i,j]-stats[4,i,j])/stats[4,i,j] *100
		stats[7,i,j] = (stats[3,i,j]-stats[5,i,j])/stats[5,i,j] *100
		thresh = 0
	ENDFOR
ENDFOR
print, 'done'
END

function MSEdenoise, thresh
	common NoisyImage, in, noisy, sigma
	clean = DIcal_denoise(noisy, sigma=sigma, thresh=thresh)
	return, [total((in-clean)^2)]
END