bmanalyze( varargin )

function h = bmanalyze( varargin )
% BMANALYZE Run analysis routines on BrainMaps.org high resolution tiled brain image data.  
%   H = BMANALYZE(url) applies granulometry analysis on a per image
%   tile basis to brain image at url at the highest level of resolution.
%
%   H = BMANALYZE(url,function) applies 'function' on a per image
%   tile basis to brain image at url at the highest level of resolution.
%   Current functions are 'granulometry', 'mean', and 'var'.  It is
%   straightforward to add your own functions by editing the code below.
%
%   H = BMANALYZE(url,function,level) applies 'function' on a per image
%   tile basis to brain image at url at the given level of resolution.
%
%   Notes
%   -------
%   BMANALYZE may take several hours or even days to run.  If you want to see the type
%   of results BMANALYZE produces, try BMANALYZEPLOT.
%
%   To obtain image url's to use with BMANALYZE, use BMANALYZEGETPATHS, or
%   alternatively, find a brain dataset you
%   want to use at http://brainmaps.org, click on the 'Metadata' link, and
%   then click on 'Generate Paths' to generate a list of all url's for the
%   selected brain dataset.
%
%   Example
%   -------
%   Run Granulometry analysis on coronal mouse brain Nissl section.
%
%       bmanalyze('http://brainmaps.org/HBP2/mus.musculus/cor/002-nissl/m08a/', 'granulometry');
%       bmanalyzePlot;
%
%   See also BMANALYZEPLOT, BMANALYZEINFO, BMANALYZEGETPATHS.
%
%   Copyright 2006 BrainMaps.org.

try
    url = varargin{1}
catch
	error('try "help bmanalyze"')
end

a_pre=imread(strcat(url,'TileGroup0/0-0-0.jpg'));
figure,imagesc(a_pre),axis equal

tic

url
url2=strcat(url,'ImageProperties.xml');

 

theStruct = parseChildNodes(xmlread(url2));

full_width=str2num(theStruct.Attributes(6).Value)
full_height=str2num(theStruct.Attributes(1).Value)
numtiles=str2num(theStruct.Attributes(3).Value)
tile_size=256;

try
    num_resolutions = varargin{3}
catch
	num_resolutions = ceil(log(ceil(max(full_width,full_height)/tile_size))/log(2.0));
end

num_resolutions_max = ceil(log(ceil(max(full_width,full_height)/tile_size))/log(2.0))

 

%This next loop enumerates over all image tiles and saves to the array
%'tilehold'

tilehold=single([]);
tile_count=0;

for res=1:num_resolutions
    res
    for col=0:floor(full_height/(256*(2^(num_resolutions_max-res)) )) 
        for row=0:floor(full_width/(256*(2^(num_resolutions_max-res)) )) 
           tilehold(tile_count+1,1:4)=single([floor((tile_count+1)/256), res, row , col]);
           tile_count=tile_count+1;                    
        end
    end
end
[x y]=size(tilehold);

out=zeros(floor(full_width/256),floor(full_height/256));
specount=0;
spetotal=[];


        
try
    analysisType = varargin{2}
catch
    analysisType = 'granulometry'
end

%This next loop performs the image analysis on a per tile basis.

for tileindex=1:x

    if (tilehold(tileindex,2)==num_resolutions)
        url3=strcat(url, 'TileGroup', num2str(tilehold(tileindex,1)),'/',num2str(tilehold(tileindex,2)),'-',num2str(tilehold(tileindex,3)),'-',num2str(tilehold(tileindex,4)),'.jpg')
        
        try
            a=imread(url3);
        catch
            url3=strcat(url, 'TileGroup', num2str(tilehold(tileindex,1)-1),'/',num2str(tilehold(tileindex,2)),'-',num2str(tilehold(tileindex,3)),'-',num2str(tilehold(tileindex,4)),'.jpg')
            try
                a=imread(url3);
            catch
                url3=strcat(url, 'TileGroup', num2str(tilehold(tileindex,1)+1),'/',num2str(tilehold(tileindex,2)),'-',num2str(tilehold(tileindex,3)),'-',num2str(tilehold(tileindex,4)),'.jpg')
                a=imread(url3);
            end
        end
        

        
        
        switch analysisType
            case 'mean'
                out( tilehold(tileindex,3)+1, tilehold(tileindex,4)+1 ) = mean(mean(rgb2gray(a))'); 
            case 'variance'
                out( tilehold(tileindex,3)+1, tilehold(tileindex,4)+1 ) = var(mean(rgb2gray(a))');
            case 'granulometry'
                

   
                I=-im2double(rgb2gray(a))+1;
                for counter = 0:36
                    remain = imopen(I, strel('disk', counter));
                    intensity_area(counter + 1) = sum(remain(:));
                end
                intensity_area_prime= diff(intensity_area);
%               figure,plot(intensity_area_prime, 'm - *'), grid on;
%               title('Granulometry (Size Distribution) of Cell Bodies');
%               set(gca, 'xtick', [0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36]);
%               xlabel('radius of cell bodes (pixels)');
%               ylabel('Sum of pixel values in cell bodies as a function of radius');
%               xlim([0 36])

                spetotal(specount+1,1:36)=intensity_area_prime;
                specount=specount+1;
        
        end  %switch statement
        
        
        
    end
end

%Saving results for use with BMANALYZEPLOT
save BMANALYZE



%Plotting some results.  Use BMANALYZEPLOT for more plots.
switch analysisType
    case 'mean'
        
    case 'variance'
    
    case 'granulometry'
        figure,imagesc(spetotal),colorbar,xlabel('radius of cell bodies (pixels)')
        %figure,imagesc(out'),colormap(gray),colorbar, axis equal

        out2=[];  %repmat(zeros(floor(full_width/256),floor(full_height/256)), 5, 5 );

        for cellrad=1:36
            specount=0;
            out=zeros(floor(full_width/256),floor(full_height/256));
            for tileindex=1:x
                if (tilehold(tileindex,2)==num_resolutions)
                    out( tilehold(tileindex,3)+1, tilehold(tileindex,4)+1 ) = spetotal(specount+1,cellrad);
                    specount=specount+1;
                end
            end

            out2(  rem(cellrad-1,6)*floor(full_width/256)+1:(rem(cellrad-1,6)+1)*floor(full_width/256)+1, (ceil(cellrad/6)-1)*floor(full_height/256)+1: (ceil(cellrad/6))*floor(full_height/256)+1 )=out;%-(cellrad+1)*10000000;

        end

        [x2 y2]=size(out2);
        figure,imagesc(out2'), colormap(gray),colorbar, axis equal, ylim([0 y2])

end


toc

















%The functions below are for parsing XML



% ----- Subfunction PARSECHILDNODES -----
function children = parseChildNodes(theNode)
% Recurse over node children.
children = [];
if theNode.hasChildNodes
   childNodes = theNode.getChildNodes;
   numChildNodes = childNodes.getLength;
   allocCell = cell(1, numChildNodes);

   children = struct(             ...
      'Name', allocCell, 'Attributes', allocCell,    ...
      'Data', allocCell, 'Children', allocCell);

    for count = 1:numChildNodes
        theChild = childNodes.item(count-1);
        children(count) = makeStructFromNode(theChild);
    end
end

% ----- Subfunction MAKESTRUCTFROMNODE -----
function nodeStruct = makeStructFromNode(theNode)
% Create structure of node info.

nodeStruct = struct(                        ...
   'Name', char(theNode.getNodeName),       ...
   'Attributes', parseAttributes(theNode),  ...
   'Data', '',                              ...
   'Children', parseChildNodes(theNode));

if any(strcmp(methods(theNode), 'getData'))
   nodeStruct.Data = char(theNode.getData); 
else
   nodeStruct.Data = '';
end

% ----- Subfunction PARSEATTRIBUTES -----
function attributes = parseAttributes(theNode)
% Create attributes structure.

attributes = [];
if theNode.hasAttributes
   theAttributes = theNode.getAttributes;
   numAttributes = theAttributes.getLength;
   allocCell = cell(1, numAttributes);
   attributes = struct('Name', allocCell, 'Value', allocCell);

   for count = 1:numAttributes
      attrib = theAttributes.item(count-1);
      attributes(count).Name = char(attrib.getName);
      attributes(count).Value = char(attrib.getValue);
   end
end