Compute Curvatures and Examine Results (in 2D)

curvature_script.m

% Read the starting image into Matlab
bwData=MIJ.getCurrentImage > 0;
%% setup and run curvature
% calculate the labels
labelImage=bwlabel(bwData);

% Calculate the Curvatures
MIJ.run('Compute Curvatures','compute sigma=2 use order')
% Read the output image into Matlab as a Stack 
curvData=MIJ.getCurrentImage;
% (1st slice is highest eigenvalue, 2nd is second)
curv1=curvData(:,:,1);
curv2=curvData(:,:,2);

% transform to a list
[xx,yy]=meshgrid(1:size(curvData,1),1:size(curvData,2));
% find the non zero values inside objects
nonZeroPts=find(sum(abs(curvData)>0,3)>0 & bwData);
% keep the nonzero values and store the results in a list with
% x position, y position, label, principal curvature, secondary curvature
curvList=[xx(nonZeroPts) yy(nonZeroPts) labelImage(nonZeroPts) curv1(nonZeroPts) curv2(nonZeroPts)];

%% 3d histogram
% create 2D histograms of the principal and secondary curvatures
% from the list
hist3(curvList(:,4:5))
pause(1.0);

%% curvature histogram
% now create a histogram for each component
% we use absolute value (abs) to only focus on the magnitude of the
% curvature
% this command just creates the binsizes so they are all on the same grid
[binCounts,binSizes]=hist3((curvList(:,4:5)));
% the number of components
numComps=max(curvList(:,3));
% the expected maximum value based on the counts for the whole image
expMax=max(binCounts(:))/numComps;
figure(2);
subplot(round(numComps/3),3,1);
for curComp = 1:numComps
    subplot(round(numComps/3),3,curComp);
    curCounts=hist3((curvList(curvList(:,3)==curComp,4:5)),binSizes);
    imagesc(binSizes{1},binSizes{2},curCounts/max(curCounts(:)),[0,1]);
    xlabel('Primary Curvature');
    ylabel('Secondary Curvature');
    % use colorful colormap
    colormap jet;
    title([ ' Component: #' num2str(curComp)]);
end


%% show shapes vs curvature
figure(3);
subplot(2,numComps,1);
for curComp = 1:numComps
    subplot(2,numComps,curComp);
    
    imagesc(labelImage==curComp)
    title([ ' Component: #' num2str(curComp)]);
    subplot(2,numComps,numComps+curComp);
    curCounts=hist3((curvList(curvList(:,3)==curComp,4:5)),binSizes);
    imagesc(binSizes{1},binSizes{2},curCounts/max(curCounts(:)),[0,1]);
    xlabel('Primary Curvature');
    ylabel('Secondary Curvature');
    % use colorful colormap
    colormap jet;
    title([ ' Curvatures: #' num2str(curComp)]);
end


%% show shape vs primary curvature
figure(4);
subplot(2,numComps,1);
for curComp = 1:numComps
    subplot(2,numComps,curComp);
    
    imagesc(labelImage==curComp)
    title([ ' Component: #' num2str(curComp)]);
    subplot(2,numComps,numComps+curComp);
    curCounts=hist((curvList(curvList(:,3)==curComp,4)),binSizes{1});
    plot(binSizes{1},curCounts/sum(curCounts),'r-.')
    xlabel('Primary Curvature');
    ylabel('Frequency');
    % use colorful colormap
    title([ ' Curvature: #' num2str(curComp)]);
end