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
    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
    //Generate error seed
    Duplicate/Free/o predicted, Perror
    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   
    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


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




Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More