Standard deviation in ROI of image
sjr51
I am calculating signal to noise ratio of spots in a bunch of images. I am specifying the ROI for the signal and an ROI for the noise using a mask wave where 0 is outside the ROI and 1 is inside.
Calculating the mean can be done using a nice one-liner from AG that I found in another thread. I'm curious if there is a similar one line MatrixOp for getting the standard deviation.
Here is my current code:
// snip
m0 = 0 // m0 is mask wave for signal ROI
MakeMaskMean(m0,xx,yy) // external function
MatrixOP/O aa = sum(m0 * imgMat) / sum(m0) // calculates the mean
qW[i][0] = aa[0]
m1 = 0 // m1 is mask wave for noise ROI
MakeMaskSD(m1,xx,yy) // external function
MatrixOP/O bb = sum(m1 * imgMat) / sum(m1) // get mean
Duplicate/O/FREE imgMat,m2
m2 -=bb[0] // subtract mean from image
m2[][] *=m1[p][q] // apply the mask
MatrixOP/O cc = sqrt(sum(m2 * m2) / sum(m1)) // calculate SD
qW[i][1] = cc[0]
// snip
m0 = 0 // m0 is mask wave for signal ROI
MakeMaskMean(m0,xx,yy) // external function
MatrixOP/O aa = sum(m0 * imgMat) / sum(m0) // calculates the mean
qW[i][0] = aa[0]
m1 = 0 // m1 is mask wave for noise ROI
MakeMaskSD(m1,xx,yy) // external function
MatrixOP/O bb = sum(m1 * imgMat) / sum(m1) // get mean
Duplicate/O/FREE imgMat,m2
m2 -=bb[0] // subtract mean from image
m2[][] *=m1[p][q] // apply the mask
MatrixOP/O cc = sqrt(sum(m2 * m2) / sum(m1)) // calculate SD
qW[i][1] = cc[0]
// snip
Getting SD is 5 lines and I wonder if it would be faster to push the values in the ROI to a 1D wave and do WaveStats, or if there is some MatrixOP magic I am missing.
Thanks for any help.
Imagestats returns average and standard deviation and takes an ROI mask all in one line.
Andy
March 29, 2017 at 08:43 am - Permalink
March 29, 2017 at 08:55 am - Permalink
MakeMaskMean(m0,xx,yy) // external function
MatrixOP/O aa = sum(m0 * imgMat) / sum(m0) // calculates the mean
qW[i][0] = aa[0]
m1 = 0 // m1 is mask wave for noise ROI
MakeMaskSD(m1,xx,yy) // external function
MatrixOP/O bb = sum(m1 * imgMat) / sum(m1) // get mean
Variable bb0=bb[0]
// There is no need to waste time on this duplication; MatrixOP takes care of it:
// Duplicate/O/FREE imgMat,m2
// m2 -=bb[0] // subtract mean from image
// m2[][] *=m1[p][q] // apply the mask
// MatrixOP/O cc = sqrt(sum(m2 * m2) / sum(m1)) // calculate SD
MatrixOP/O cc=sqrt(sum(magsqr((imgMat-bb0)*m1))/sum(m1))
qW[i][1] = cc[0]
In practice, depending on the size of your images and ROI's, ImageStats may be more efficient because in IP7 ImageStats is automatically multithreaded.
Also note that when you work with ROI in MatrixOP, included pixels correspond to a value of 1 in the mask while ImageStats requires that included pixels are marked by 0 in the ROI wave.
A.G.
WaveMetrics, Inc.
March 30, 2017 at 10:40 am - Permalink
Thanks!
March 30, 2017 at 01:08 pm - Permalink