matlab_cc.m

%
%	A MATLAB program for automated analysis of vertical canopy images	
%
% 	Authors: Jaakko Heikkinen & Lauri Korhonen 2006-2009
%	Citation: Korhonen, L. & Heikkinen, J. 2008. Automated analysis of in situ canopy images for the estimation of forest canopy cover. 
%		Forest Science 55(4): 323-334.
%	Thresholding algorithm: Nobis, M. & Hunziker, U. 2005. Automatic thresholding for hemispherical canopy photographs based on 
%		edge detection. Agricultural and forest meteorology 128: 243-250.
%	This version analyzes 4:3 area of the image determined by hor_angle parameter.
%	If the thresholding fails, you can try adjusting min- or max_thresh and sky parameters to make binarization work.

%Constants

filetype='jpg';	%Which type of files are read
distance=8;		%Distance used in defining whether the maxima are real
sky=11; 		%If the maximum brightness difference between sky and canopy is smaller than this value, the image is supposed to be pure sky!
medf_dist=7; 	%Distance used in the median filtering of the maxima
min_thresh=60;
max_thresh=220;	%These define what range of 8-bit values (0-255) is analyzed. Speeds up the measurement and eliminates some extra peaks.
elem_size=10;	%Diameter of the circular structuring element used in painting; 10 OR for 640x480 resolution images
true_focal=6.5;	%Focal length of the camera, mm
ccd_hor=6.16;	%Horizontal size of the ccd sensor
ccd_ver=4.62;	%Vertical size of the ccd sensor
hor_angle=90 	%Horizontal angle of view that should be used, degrees. Larger than cameras AOV -> whole image is analyzed.
thresh=180;		%Initial guess for the threshold

files=dir(['*.', filetype]); 	%Read all files of the specified type

results=cell(length(files),4);	%Create array for the results
results{1,1}='Image';
results{1,2}='Threshold'; 
results{1,3}='C_Closure'; 
results{1,4}='C_Cover'; 

disp('Analyzing image:');

%%%%%%%%%%%%%%% START LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for h=1:length(files) 	%Loop goes through all images h in directory
    
	fileid=files(h).name(1:length(files(h).name)-(length(filetype)+1));
	disp(fileid);
	
	rgb=imread(files(h).name); 		%Read RGB image
	[rows,cols,channels]=size(rgb);	%Read image dimensions
	if cols ~=640					%Image is scaled to 640x480 so that size of the structuring element is correct
		rgb = imresize(rgb,[480 640],'nearest');
	end
	
	
	subplot(221) 
	image(rgb); %Show RGB-image in upper left subplot
	title(files(h).name,'FontSize',14,'Interpreter','none');
	axis image
	hold on;
	[rows,cols,channels]=size(rgb); %Save image dimesions
	
	orig_angle=2*atan(ccd_hor/(2*true_focal))*180/pi; %Calculate horizontal angle of view in degrees
	
	blue=rgb(:,:,3);	%Extract blue channel
	blue=medfilt2(blue); %Median filtering to decrease noise
	
	
	%%%%%%%%%%%%%%%  THRESHOLDING  %%%%%%%%%%%%%%%%%%%%
	
	blue11=blue(1:(rows-1),1:(cols-1));	%Left top corner subpixel matrix
	blue12=blue(1:(rows-1),2:cols);		%Right top corner subpixel matrix
	blue21=blue(2:rows,1:(cols-1));		%Left bottom corner subpixel matrix
	blue22=blue(2:rows,2:cols);			%Right bottom corner subpixel matrix

	bin=false(size(blue));   			%Make empty 0-array
	difftable=zeros(256,1);          	%Empty vector for later use
	
	for k=min_thresh:max_thresh	%Loop through 8-bit blue channel value range
		
		bin=(blue<=k);	%Binarization with the current threshold

		bin11=bin(1:(rows-1),1:(cols-1));	%Subimages for comparing the pixels
		bin12=bin(1:(rows-1),2:cols);
		bin21=bin(2:rows,1:(cols-1));
		bin22=bin(2:rows,2:cols);
		
		%Pixels 11 vs 12
		is_difference=bin11~=bin12;
		dif11_12=abs(blue11(is_difference)-blue12(is_difference));
		
		%Pixels 11 vs 21
		is_difference=bin11~=bin21;
		dif11_21=abs(blue11(is_difference)-blue21(is_difference));
		
		%Pixels 11 vs 22
		is_difference=bin11~=bin22;
		dif11_22=abs(blue11(is_difference)-blue22(is_difference));
		
		%Pixels 12 vs 21
		is_difference=bin12~=bin21;
		dif12_21=abs(blue12(is_difference)-blue21(is_difference));
		
		%Calculate denominator
		div=(length(dif11_12)+length(dif11_21)+length(dif11_22)+length(dif12_21));
		if div ~=0
			%Calculate mean brightness diffenrence between the threshold pixels
			difftable(k)=( sum(dif11_12)+sum(dif11_21)+sum(dif11_22)+sum(dif12_21))/div;
		else
			difftable(k)=0; %No differences 
		end
		
	end 

	%Median filtering to smooth the difference histogram
	midmed=medfilt1(difftable, medf_dist);
	
	%%%%%%%%%%%% REDUCTION OF AOV FOR FURTHER ANALYSIS %%%%%%%%%%%%%%% 
	%Reduction is done in 4:3 aspect ratio, angle defined in horizontal direction
	wholeblue=rgb(:,:,3);	%Extract the blue channel for whole image for later use
	if orig_angle > hor_angle
	
		%Calculate the dimensions of the area to be analyzed
		new_cols=round(hor_angle*cols/orig_angle);
		new_rows=round(rows*new_cols/cols);	
		cmin=round(cols/2-new_cols/2);
		cmax=round(cols/2+new_cols/2);
		rmin=round(rows/2-new_rows/2);
		rmax=round(rows/2+new_rows/2);	
		
		blue=rgb(rmin:rmax,cmin:cmax,3);	%Extract the blue channel from the area to be analyzed
		[rows,cols]=size(blue); 			%New image dimensions
		
		%Delineate the area to be analyzed in the image
		frame=[cmin,rmin;cmax,rmin;cmax,rmax;cmin,rmax;cmin,rmin];
		plot(frame(:,1),frame(:,2),'LineWidth',1,'Color','r');
		hold off;
	end
		

	%%%%%%%%%% DRAW THRESHOLD PLOT %%%%%%%%%%%%%%%%
	
	subplot(222)
	bar(0:255,difftable,'FaceColor',[0.5 0.5 0.5],'EdgeColor',[0.5 0.5 0.5]);
	hold on

	plot(0:255,midmed, 'LineWidth',2,'Color','black');
		
	%The maximum difference is accepted only if the values at the distance of 8 and 16 are smaller than the peak value   
	is_maximum=false(256,1);

	%Find the possible maxima
	for top=min_thresh:max_thresh
		midleft=midmed(top-distance);
		midright=midmed(top+distance);
		extremeleft=midmed(top-2*distance);
		extremeright=midmed(top+2*distance);
		is_maximum(top)= (extremeleft<midleft) & (midleft<top)...
			& (midright<top) & (extremeright<midright);
	end
	possible_maxima=find(is_maximum==true);	
	
	if ~ isempty(possible_maxima) 
		%If possible maxima were found, find the largest one
		[medimax,max_index]=max(midmed(possible_maxima));

		if medimax>sky	%The brightness difference indicates the image is not only sky
		
			thresh=possible_maxima(max_index);	%Final threshold is determined
			
			%Finish the threshold plot
			plot(thresh,midmed(thresh),'Marker','o','MarkerSize',12,...
				'LineWidth',2,'MarkerFaceColor','white','MarkerEdgeColor','black');
			text(thresh-15,1.12*midmed(thresh),num2str(thresh),'FontSize',12);
			set(gca,'XLim',[0 255],'YLim',[0 1.3*midmed(thresh)]);
			title(strcat('Threshold = ',num2str(thresh)),'FontSize',13);
		else
			set(gca,'XLim',[0 255],'YLim',[0 1.3*max(midmed)]);
			disp('Empty image')
		end
	else
		set(gca,'XLim',[0 255],'YLim',[0 1.3*max(midmed)]);
	end
	hold off

	
	%%%%%%%%%%%%%%% CANOPY CLOSURE %%%%%%%%%%%%%%%%%%%%%%%%%
   
	binim=ones(size(blue));	%Initialize result image with white
	wholebinim=ones(size(wholeblue));	%Initialize result image with white
	%If a maximum was found and the image is not just sky, create a binary image
	if (~isempty(possible_maxima) && medimax>sky)
		binim(blue<thresh)=0;
		wholebinim(wholeblue<thresh)=0;
	end 
	%Otherwise use threshold from previous image
	
	%Calculate canopy closure
	whites=sum(sum(binim));
	cclosure=(1-whites/(rows*cols))*100;	
	%Print binary image
	colormap(gray(2));	%Set colors	
	subplot(223)
	set(gca,'YDir','normal');   
	image(binim*2);	%Plot image, multiplication needed to see the color differences
	axis image
	title(strcat('Canopy closure = ',num2str(round(cclosure)),'%'),'FontSize',13);
	
	
	%%%%%%%%%%%%%%%%%  CANOPY COVER %%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	%Paint the crowns black using morphological opening and closing
	elem=strel('disk',elem_size);
	opeclo=imopen(binim,elem);	%Because objects are black, do the closing with opening operation
	opeclo=imclose(opeclo,elem);

	%Calculate pixel sum
	whites=sum(sum(opeclo));
	ccover=(1-whites/(rows*cols))*100;
	
	%Plot painted image
	subplot(224)
	image(opeclo*2);	%Plot image, multiplication needed to see the color differences
	axis image
	title(strcat('Canopy cover = ',num2str(round(ccover)),'%'),'FontSize',13);
	orient portrait
	
	%%%%%%%%%%%%%% SAVE BIN AND CHECK IMAGES%%%%%%%%%%%%%%%%%%%%%
    
	%Save check image as jpg 
	print('-djpeg','-r100',['check_',files(h).name])
	
	%Save thresholded and painted images as bmp
	binim=binim>0; %Convert to 2-bit
	imwrite(binim, ['abin_',fileid,'_',num2str(hor_angle),'.bmp'],'bmp')	%Only analyzed part, binary
	opeclo=opeclo>0;
	imwrite(opeclo, ['painted_',fileid,'_',num2str(hor_angle),'.bmp'],'bmp')	%Only analyzed part, painted
	wholebinim=wholebinim>0;
	imwrite(wholebinim, ['bin_',fileid,'.bmp'],'bmp')	%Whole image, binary
	close all
	
	%Add results to array
	results{h+1,1}=files(h).name;
	results{h+1,2}=thresh;
	results{h+1,3}=cclosure;
	results{h+1,4}=ccover;
	
end
%Show and save results
results
xlswrite(['Matlab_CC',num2str(hor_angle),'.xls'],results);

clear all %Clear memory