Kalman Stack filter

/////////////////////////////////////////////////////////////
///////////// Kalman filter for time series image stacks ///////
////the prediction bias determins the amount of filtering//////////////
//////the noise estimate has little effect on the outcome/////////////
//////////////////////////////////////////////////////////////////////////////
//Inserts a menu call in the Analysis menu, or can be called by typing 
// "kalman()" in the command line, it will only find 3 dimensional waves and outputs the filtered wave
// to a new 16bit wave with "_Kal" appended to the original file name
//Based on the imageJ code by Christopher Phillip Mauer

function kalman()

	//Initialize variables and input waves
	string inputwave
	Variable G = 0.8		//Filter gain
	Variable V = 0.06		// error estimate
	
	string list=wavelist("*",";","DIMS:3")
	prompt inputwave, "Wave Select", popup list
	prompt G, "Prediction Bias"
	prompt V, "Noise Estimate"
	doprompt "Kalman Filter", inputwave, G,V
			if(V_flag==1)
			abort
		endif
			
	wave input=$inputwave			// reference to wave being filtered
	redimension/s input			// convert to single floating point
	
	variable z = dimsize(input,2)	//counts the number of frames
	
	//Generate prediction seed
	Duplicate/FREE/o input, predicted
	redimension/N=(-1,-1,0) predicted
	multithread predicted=input[p][q][0]
	
	//Generate other variables
	duplicate/FREE/o predicted, one, observed
	one=1
	
	//Generate error seed
	Duplicate/Free/o predicted, Perror
	Perror=V
	Duplicate/FREE/o Perror, errorEst
	
	//Generate ouput wave
	Duplicate/FREE/o input, output
	
	variable i
	for(i=0;i<z;i+=1)		//Do filter
		multithread	observed=input[p][q][i]		//Get observed values
		
		matrixop/FREE/o Kalman=Perror/(Perror+errorEst)			//Calculate Kalman gain
		
		matrixop/FREE/o corrected= g*predicted+(1-g)*observed+kalman*(observed-predicted)		//calcuate corrected image
		
		matrixop/FREE/o correctedError = Perror*(one-kalman)					// calculate corrected noise estimate 
		
		matrixop/o Perror = correctedError										//Update predicted noise 
		Matrixop/o predicted = corrected											//Update prediction
		multithread output[][][i]=corrected[p][q]									//append corrected image to output stack	
		
	endfor
	
	redimension/w output
	String outputname=inputwave+"_Kal"	//Give filtered wave new name
	Duplicate/o output, $outputname
	
	newimage/k=1 $outputname 		//Display filtered wave with slider
	WMAppend3DImageSlider();

end

menu  "Analysis"
			"Kalman Filter", kalman()
End