https://github.com/quantixed/VolumeFinder
Tip revision: 708b6d9a96330da323d07708274e68b0a8bf8615 authored by Stephen Royle on 09 March 2017, 11:37:32 UTC
SNR added
SNR added
Tip revision: 708b6d9
VolumeFinder.ipf
#pragma TextEncoding = "MacRoman" // For details execute DisplayHelpTopic "The TextEncoding Pragma"
#pragma rtGlobals=3 // Use modern global access method and strict wave access.
// Menu item for easy execution
Menu "Macros"
Submenu "Volume Finder"
"Find Volumes", VolumeFinder()
End
End
// Use this function to load and process TIFF stacks in one folder
// TIFFs must be binarized versions of amira segmentation files
Function VolumeFinder()
String expDiskFolderName, expDataFileName
String FileList, ThisFile, wName
Variable FileLoop, nWaves, i
NewPath/O/Q/M="Please find disk folder" ExpDiskFolder
if (V_flag!=0)
DoAlert 0, "Disk folder error"
Return -1
endif
PathInfo /S ExpDiskFolder
ExpDiskFolderName=S_path
FileList=IndexedFile(expDiskFolder,-1,".tif")
Variable nFiles=ItemsInList(FileList)
Variable /G vol
Variable /G nI
Make/O/D/N=(nFiles) volWave, volHWave
Make/O/T/N=(nFiles) fileWave
Make/O/D/N=(nFiles) nPointWave
for (FileLoop = 0; FileLoop < nFiles; FileLoop += 1)
ThisFile = StringFromList(FileLoop, FileList)
expDataFileName = ReplaceString(".tif",ThisFile,"") // get rid of .tif
expDataFileName = ReplaceString(".labels",expDataFileName,"") // get rid of .labels
ImageLoad/O/T=tiff/C=-1/LR3D/Q/P=expDiskFolder ThisFile
VolumeCalc($ThisFile)
fileWave[FileLoop] = expDataFileName
volWave[FileLoop] = vol
nPointWave[FileLoop] = nI
Wave kCloud // path of the 3D convex hull
wName = expDataFileName + "_kCloud"
Rename kCloud $wName
Wave HullWave // 1D wave of polygonarea for each 2D hull per layer
volHWave[FileLoop] = sum(HullWave)
wName = expDataFileName + "_HullWave"
Rename HullWave $wName
Wave pCloud // point cloud of all pixels over threshold
wName = expDataFileName + "_pCloud"
Rename pCloud $wName
Wave LayerWave // number of pixels over threshold per layer
wName = expDataFileName + "_LayerWave"
Rename LayerWave $wName
wName = expDataFileName + "_LDens"
MatrixOp/O $wName = LayerWave / HullWave
WaveTransform zapnans $wName
KillWaves /Z $ThisFile //should already be killed
endfor
Duplicate nPointWave densityWave
densityWave /= volWave
//Display results
DoWindow/K Results
Edit/N=Results fileWave,nPointWave,volWave,densityWave,volHWave
DoWindow/K MTvol
Display/N=MTvol nPointWave vs fileWave
ModifyGraph swapXY=1
SetAxis/A/R left
SetAxis/A/E=1/N=1 bottom
Label bottom "Point Volume"
DoWindow/K Spindlevol
Display/N=spindlevol volWave vs fileWave
ModifyGraph swapXY=1
SetAxis/A/R left
SetAxis/A/E=1/N=1 bottom
Label bottom "Hull Volume"
DoWindow/K Density
Display/N=densityvol densityWave vs fileWave
ModifyGraph swapXY=1
SetAxis/A/R left
SetAxis/A/E=1/N=1 bottom
Label bottom "Density"
DoWindow/K summaryLayout
NewLayout/N=summaryLayout
AppendLayoutObject/W=summaryLayout graph MTvol
AppendLayoutObject/W=summaryLayout graph spindlevol
AppendLayoutObject/W=summaryLayout graph densityvol
#If igorversion()>=7
LayoutPageAction size(-1)=(595, 842), margins(-1)=(18, 18, 18, 18)
#EndIf
ModifyLayout units=0
ModifyLayout frame=0,trans=1
Execute /Q "Tile"
ScaleIt(12,12,60)
End
// This function does the calculations for each dataset
/// @param m0 matrix for processing
Function VolumeCalc(m0)
Wave m0
Variable timer = startmstimer
NVAR /Z vol
NVAR /Z nI //global variables
nI = Dimsize(m0,0)
Variable nJ = Dimsize(m0,1)
Variable nK = Dimsize(m0,2)
Make/O/N=((nI*nJ*nK),3) m1 = NaN //slows code but worth it
Variable l = 0 // rows in result wave
vol = 0
Variable tempvar, V_value
Variable i, j, k
For (k = 0; k < nK; k += 1) // layers
For (j = 0; j < nJ; j += 1) // columns
For (i = 0; i < nI; i += 1) // rows
If (m0[i][j][k] != 0)
m1[l][0] = i
m1[l][1] = j
m1[l][2] = k
l += 1
EndIf
EndFor
EndFor
EndFor
MatrixOp/O w0 = sumRows(m1)
WaveTransform ZapNans w0
nI = numpnts(w0)
Duplicate/O/R=(0,nI-1) m1, pCloud
KillWaves w0, m1, m0 // also gets rid of image to free memory
Duplicate/O/R=[*][2] pCloud, pCloudZ
Redimension/N=-1 pCloudZ
Make/O/N=(nK) HullWave = 0, LayerWave = 0
// HullWave will hold the area of each convex hull per layer
// LayerWave will hold the number of pixels that are over threshold in that layer
String wName, wList
For (k = 0; k < nK; k += 1) // layers
FindValue/V=(k) /Z pCloudZ
If(V_Value > -1) // test to see if it's worth calculating this layer
Duplicate/O/R=[*][0] pCloud, pCloudX
Duplicate/O/R=[*][1] pCloud, pCloudY
Redimension/N=-1 pCloudX, pCloudY
pCloudX = (pCloudZ == k) ? pCloudX : NaN
pCloudY = (pCloudZ == k) ? pCloudY : NaN
WaveTransform zapnans pCloudX
WaveTransform zapnans pCloudY
LayerWave[k] = numpnts(pCloudX)
Convexhull /C pCloudX, pCloudY
Wave W_XHull, W_YHull
HullWave[k] = PolygonArea(W_XHull,W_YHull)
Duplicate/O W_XHull, W_ZHull
W_ZHull = k
wName = "kCloud_" + num2str(k)
Concatenate/O/KILL {W_XHull,W_YHull,W_ZHull}, $wName
EndIf
EndFor
wList=WaveList("kCloud_*",";","")
Concatenate/O/KILL/NP=0 wList, kCloud
Triangulate3D/VOL kCloud
vol = V_Value
KillWaves pCloudX, pCloudY, pCloudZ
printf "%g\r", stopmstimer(timer)/1e6
End
/// @param xnm voxel size, x dimension in nm, i.e. 12
/// @param ynm voxel size, y dimension in nm, i.e. 12
/// @param znm voxel size, z dimension in nm, i.e. 60
Function ScaleIt(xnm,ynm,znm)
Variable xnm, ynm, znm
//This will scale the points to real world values
Variable scale = (xnm * ynm * znm) / (10^9) //in µm^3
//need to scale MTs in a different way
If(xnm != ynm)
Print "xnm and ynm are not equal. Please check"
EndIf
Variable MTscale = ((1/3) * (xnm * ((PI*12.5)^2))) / (10^9)
// assumes MTs are 3 px wide and are single MTs
// scale to µm^3
Wave nPointWave, volWave
nPointWave *=MTscale
volWave *=scale
Label /W=MTvol bottom, "Point Volume (µm\S3\M)"
Label /W=spindlevol bottom, "Hull Volume (µm\S3\M)"
End