% Author: Mihai Ivanovici, MIV Imaging Venture
% Affiliation: Department of Electronics and Computers, Transilvania University of Brasov, Romania
%
% Date: July 2008
% Last update: June 2010
%
% This Matlab function implements the extension to the colour domain of the Voss probabilistic algorithm for 
% the estimation of fractal dimension and computes the lacunarity curve for a given colour image.
% Input parameters: the name of the file containing the colour image
% Outpur parameters: the figure representing the colour lacunarity curve
% The run of this Matlab script may take up to several minutes for a 500 x 500 image.
%
% You may use the code as you wish for academic purposes, as long as you acknowledge our work by referencing the following article: 
% M. Ivanovici, N. Richard, "The Lacunarity of Colour Fractal Images", ICIP’09, the International Conference on Image Processing, Cairo, Egypt, November 7-11, 2009.

function l = color_lacunarity( imgfilename )
	LMAX = 41; 			% the maximum size of the hyper-cubes
	local_dims = 3 : 2 : LMAX; 	% the sizes of the hyper-cubes
	mmax = power( LMAX, 2 );	% the maximum number m of pixels

	P = zeros( mmax, length( local_dims ) );	% the probability matrix
	M = zeros( 1, length( local_dims ) );		% the first order moment
	M2 = zeros( 1, length( local_dims ) );		% the second order moment
	Lambda = zeros( 1, length( local_dims ) );	% the lacunarity

	img = imread( imgfilename );
	[ n, m, o ] = size( img );

	skip = floor( LMAX / 2 ) + 1;
	for i = skip : n - skip
	       for j = skip : m - skip
	       		for k = 1 : length( local_dims )
				cube_size = local_dims( k );
				
				%count how many pixels fall inside the hyper-cube
				cpix = img( i, j, 1 : 3 );
				cs = floor( cube_size / 2 );
				region = img( ( i - cs ) : ( i + cs ), ( j - cs ) : ( j + cs ), 1 : 3 );
				npix = length( find( ...	
						region( :, :, 1 ) >= cpix( 1 ) - cs & region( :, :, 1 ) <= cpix( 1 ) + cs &	...
						region( :, :, 2 ) >= cpix( 2 ) - cs & region( :, :, 2 ) <= cpix( 2 ) + cs &	...
						region( :, :, 3 ) >= cpix( 3 ) - cs & region( :, :, 3 ) <= cpix( 3 ) + cs ...
					) );

				P( npix, k ) = P( npix, k ) + 1;
			end
	       end
	end

	%normalize the probability matrix
	P = P / ( (n - LMAX) * (m - LMAX) );
	
	%compute the statistical moments and lacunarity \Lamba(\delta)
	for i = 1 : length( local_dims )
		for m = 1 : mmax;
			M( i ) = M( i ) + m * P( m, i );
			M2( i ) = M2( i ) + m * m * P( m, i );
		end
		Lambda( i ) = ( M2( i ) - M( i ) * M( i ) ) / ( M( i ) * M( i ) ); 
	end

	% plot and save the colour lacunarity curve
	handle = figure;
	plot( local_dims, Lambda, 'b*-', 'LineWidth', 2 );
	set( gca, 'FontSize', 14 );
	xlabel( '\delta' );
	ylabel( '\Lambda(\delta)' );
	axis( [ 0 44 0 3 ] );
	grid on;
	
	pos = strfind( imgfilename, '.' );
	output_name = imgfilename( 1 : pos-1 );
	saveas( gcf, [output_name, '_colour_lacunarity.fig'], 'fig' );
	print( handle, '-dpng', [output_name, '_colour_lacunarity.png' ] );
    
    l = Lambda;