CiVol.m

function [Volumes, N] = CiVol(mask, D)
% CIVOL Intrinsic volumes of the search region in pixel.
%	CiVol(mask, D) returns the intrinsic Volumes (N=Volumes(D))
%	If mask is a scalar, it specifies the number of pixels in the search region
% 		E.g., [Volumes, N] = CiVol(pi*128^2, 2);
%	If mask is a binary matrix, where the 1's define the search region
% 		E.g.,
% 		mask = double(imread('faceMask.tif'));
% 		mask = (mask - min(mask(:))) / (max(mask(:)) - min(mask(:)));
% 		[Volumes, N] = CiVol(mask, length(size(mask)));
%	D (D<=3) is the dimensionality of the search region
% 
% From K. J. Worsley's (keith.worsley@mcgill.ca) stat_threshold help section:
%	Volumes is a vector of D+1 components, the intrinsic volumes of the search 
%	region = [Euler characteristic, 2*caliper diameter, 1/2 surface area, volume],
%	e.g. [1, 4*r, 2*pi*r^2, 4/3*pi*r^3] for a sphere of radius r in 3D. For a 2D 
%	search region, it is [1, 1/2 perimeter length, area]. The random field theory 
%	threshold is based on the assumption that the search region is a sphere, which 
%	is a very tight lower bound for any non-spherical region.
% 
% See also DEMOSTAT4CI
% 
% The Stat4Ci toolbox is free (http://mapageweb.umontreal.ca/gosselif/stat4ci.html); if you use 
% it in your research, please, cite us:
%	Chauvin, A., Worsley, K. J., Schyns, P. G., Arguin, M. & Gosselin, F. (2004).  A sensitive 
%	statistical test for smooth classification images.
% 
% Alan Chauvin & FrŽdŽric Gosselin (frederic.gosselin@umontreal.ca), 20/08/2004

if length(mask) == 1,
	
	N = mask;

else

	% the mask have to be 1 in the region of interest and 0 elsewhere
	minMask = min(mask(:));
	maxMask = max(mask(:));
	
	N = sum( (mask(:) - minMask)./  (maxMask - minMask) );
	
end;

EC = 1;

switch D
case 1

	Volumes(1) = EC; 			% Euler characteristic
		
	Volumes(2) = N;			% length
	
case 2
	
	% N = pi*r^2
	r = sqrt(N/pi);

	Volumes(1) = EC; 			% Euler characteristic
	
	Volumes(2) = pi * sqrt(N/pi)	% caliper diameter
	Volumes(2) = pi * r	% caliper diameter
	
	Volumes(3) = N;			% surface area
	
case 3

	% N = (4/3)*pi*r^3
	r = (3/4/pi*N)^(1/3);	% r is the radius of the sphere
	
	Volumes(1) = EC; 			% Euler characteristic
	
	Volumes(2) =  4 * r;		% 2 * caliper diameter
	
	Volumes(3) = 2 * pi * r;	% 1/2 surface area
		
	Volumes(4) = N;				% volume

otherwise
	error('Not implemented for D > 3.')
end;