//VolumeFractionOfBone
//ImageJ macro for the determination of bone volume fraction in
//3D image sets 
//Copyright (C) 2007 Michael Doube
//m.doube at qmul.ac.uk
//Last modified 2007-02-28 17:53 UTC

//This program is free software; you can redistribute it and/or
//modify it under the terms of the GNU General Public License
//as published by the Free Software Foundation (version 2).
//http://www.gnu.org/copyleft/gpl.html

//This program is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//GNU General Public License for more details.

//You should have received a copy of the GNU General Public License
//along with this program; if not, write to the Free Software
//Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.


var startslice;
var endslice;

macro "Start Slice [j]" {
startslice = getSliceNumber();
}

macro "End Slice [l]" {
endslice = getSliceNumber();
}

macro "Volume Fraction of Bone [q]" {
requires("1.34h");
n = getSliceNumber();
row = nResults;
vroihisto = newArray(256);
getThreshold(lower, upper);
if (lower==-1) exit("Use Image>Adjust>Threshold to set the threshold");
setBatchMode(true);
vmax = 0;
if (startslice < 1) startslice = 1;
if (endslice < 1) endslice = nSlices;
if (startslice > endslice) exit ("Start slice must be a lower number than end slice!");
print("For slice "+startslice+" to slice "+endslice);
for (i=startslice; i<=endslice; i++) {
	setSlice(i);
	getRawStatistics(count, mean, min, max, std, histogram);
	vcount = vcount + count;
	if (max > vmax) vmax = max;
	if (max < upper) top = max;
	else top = upper;
	for (m=0; m<=top;m++){
		vroihisto[m] = vroihisto[m] + histogram[m];
	}
       	for (k=0; k<lower; k++) {
            	background = background + histogram[k];
       	}
}
vfraction = (vcount - background) / vcount;

if (vmax < upper) top = vmax;
else top = upper;

av=0;
cumulative = 0;
for (q = lower; q<=top; q++){
	cumulative = cumulative + vroihisto[q];
   	av = av + q * vroihisto[q];
}
vbmmd = av / cumulative;

for (q = 0; q<=vmax; q++){
   	roiav = roiav + q * vroihisto[q];
}
roiav = roiav / vcount;

if (vfraction == 0) mode = 0/0;
else {
	mode = lower;
	for (q=lower; q<=top;q++){
		if (vroihisto[q] > vroihisto[mode]) mode = q;
	}
}
name = getTitle();
setResult("Label", row, name);
setResult("Vfb", row, vfraction);
setResult("Mean Matrix", row, vbmmd);
setResult("Mode Matrix", row, mode);
setResult("Mean ROI", row, roiav);
updateResults();
setSlice(n);
setBatchMode(false);
}