Code covered by the BSD License

# Isocontour

### Dirk-Jan Kroon (view profile)

24 Feb 2011 (Updated )

Find ISO-contour geometry in a 2D image using marching-squares, and sort the contour objects

[Lines,Vertices,Objects]=isocontour(I,isovalue)
```function [Lines,Vertices,Objects]=isocontour(I,isovalue)
% This function ISOCONTOUR computes the isocontour geometry for
% a certain 2D image and isovalue. To Extract the isocontour geometry it
% uses Marching Squares and linear interpolation. Followed by sorting the
% contour geometry into separate sorted contours.
%
% This function is the 2D equivalent of Isosurface extraction
% using Marching Cubes in 3D.
%
%
%   [Lines,Vertices,Objects]=isocontour(I,isovalue);
%
% inputs,
%   I : An 2D image (grey-scale)
%   isovalue : The Iso-value of the contour
%
% outputs,
%   Lines : An array describing all the Line-pieces of the isocontour
%           geomtery, with a N x 2 index list of vertices
%   Vertices : Vertices (Corners) of the lines M x 2 list of X,Y
%           coordinates
%   Objects : A 1 x K  cell array with in every cell a list of indices
%           corresponding to one connect isocontour. If the isocontour
%           is closed then the last index value is equal to first index
%           value.
%
% Note : This function displays the image with isocontours if no output
%       is defined.
%
% Example,
%    isocontour(I,0.5);
%
% Example,
%    [Lines,Vertices,Objects]=isocontour(I,0.5);
%    figure('renderer','opengl'), imshow(I), hold on;
%    for i=1:length(Objects)
%         Points=Objects{i};
%         plot(Vertices(Points,2),Vertices(Points,1),'Color',rand(3,1));
%    end
%
% Example,
%    [Lines,Vertices]=isocontour(I,0.5);
%    figure, imshow(I), hold on;
%    V1=Vertices(Lines(:,1),:); V2=Vertices(Lines(:,2),:);
%    plot([V1(:,2) V2(:,2)]',[V1(:,1) V2(:,1)]','b');
%    % Calculate Contour Normals
%    N = V1-V2; L = sqrt(N(:,1).^2+N(:,2).^2)+eps;
%    N(:,1)=N(:,1)./L; N(:,2)=-N(:,2)./L;
%    plot([V1(:,2) V1(:,2)+N(:,1)*5]',[V1(:,1) V1(:,1)+N(:,2)*5]','r');
%
%
% This function is written by D.Kroon University of Twente (March 2011)

% Check Inputs
if(nargin==0), error('isocontour:input','no input image defined'); end
if(ndims(I)~=2), error('isocontour:input','image must be 2D'); end
if(size(I,3)>1), error('isocontour:input','image cannot be RGB'); end
if(nargin==1), isovalue=0.5; end

% Get the Line Pieces
[Vertices,Lines]=LookupDB(I,isovalue);

% Sort the Line Pieces into objects
if(nargout==3||nargout==0)
Objects=SortLines2Objects(Lines);
end

% Show image, if no output asked
if(nargout==0)
figure, imshow(I), hold on;
for i=1:length(Objects)
Points=Objects{i};
plot(Vertices(Points,2),Vertices(Points,1),'Color',[1 0 0]);
end
end

function Objects=SortLines2Objects(Lines)
% Object index list
Obj=zeros(100,2); Obj(1,1)=1; nObjects=1; reverse=false;
for i=1:size(Lines,1)-1
F=Lines(i,2);
Q=i+1;
R=find(any(Lines(Q:end,:)==F,2),1,'first');
if(~isempty(R))
R=R+i;
TF=Lines(Q,:);
if(Lines(R,1)==F),Lines(Q,:)=Lines(R,[1 2]); else Lines(Q,:)=Lines(R,[2 1]); end
if(R~=Q), Lines(R,:)=TF; end
else
F=Lines(Obj(nObjects,1),1);
R=find(any(Lines(Q:end,:)==F,2),1,'first');
if(~isempty(R))
reverse=true;
Lines(Obj(nObjects,1):i,:)=Lines(i:-1:Obj(nObjects,1),[2 1]);
R=R+i;
TF=Lines(Q,:);
if(Lines(R,1)==F),Lines(Q,:)=Lines(R,[1 2]); else Lines(Q,:)=Lines(R,[2 1]); end
if(R~=Q), Lines(R,:)=TF; end
else
if(reverse)
Lines(Obj(nObjects,1):i,:)=Lines(i:-1:Obj(nObjects,1),[2 1]);
reverse=false;
end
Obj(nObjects,2)=i;
nObjects=nObjects+1;
Obj(nObjects,1)=i+1;
end
end
end
Obj(nObjects,2)=i+1;

% Object index list, to real connect object lines
Objects=cell(1,nObjects);
for i=1:nObjects
% Determine if the line is closed
if(Lines(Obj(i,1),1)==Lines(Obj(i,2),2))
Objects{i}=[Lines(Obj(i,1):Obj(i,2),1);Lines(Obj(i,1),1)];
else
Objects{i}=Lines(Obj(i,1):Obj(i,2),1);
end
end

function [V,F]=LookupDB(Img,isovalue)

% Describe the base-polygons by edges
%  Edge number
%   1
%  ***
% 0* *2
%  ***
%   3
%
I=zeros(16,4);
I(1,:)= [4 4 4 4]; % [0 0;0 0]
I(2,:)= [1 0 4 4]; % [1 0;0 0]
I(3,:)= [2 1 4 4]; % [0 1;0 0]
I(4,:)= [2 0 4 4]; % [1 1;0 0]
I(5,:)= [0 3 4 4]; % [0 0;1 0]
I(6,:)= [1 3 4 4]; % [1 0;1 0]
I(7,:)= [2 1 0 3]; % [0 1;1 0] ambiguous
I(8,:)= [2 3 4 4]; % [1 1;1 0]
I(9,:)= [3 2 4 4]; % [0 0;0 1]
I(10,:)=[1 0 3 2]; % [1 0;0 1] ambiguous
I(11,:)=[3 1 4 4]; % [0 1;0 1]
I(12,:)=[3 0 4 4]; % [1 1;0 1]
I(13,:)=[0 2 4 4]; % [0 0;1 1]
I(14,:)=[1 2 4 4]; % [1 0;1 1]
I(15,:)=[0 1 4 4]; % [0 1;1 1]
I(16,:)=[4 4 4 4]; % [1 1;1 1]

% The base-edges by vertex positions
E=[0 0 1 0; 0 0 0 1; 0 1 1 1; 1 0 1 1; 4 4 4 4];

% base-Polygons by vertexpostions
IE=E(I(:)+1,:);
IE=[IE(1:16,:) IE(17:32,:); IE(33:48,:) IE(49:64,:)];

% Make a Binary image with pixels set to true above iso-treshold
B=Img>=isovalue;

% Get Elementary Cells
%  Cell
%  1 ** 2
%  *    *
%  *    *
%  4 ** 8
%
B0=B(1:end-1,1:end-1); B1=B(1:end-1,2:end); B2=B(2:end,1:end-1); B3=B(2:end,2:end);
V=B0+B1*2+B2*4+B3*8+1;
[x,y]=find((V>1)&(V<16));
v=V(x+(y-1)*size(V,1));

% Elementary cells to Edge coordinates defined by connected image grid-points
J=[IE(v,:);IE(v+16,:)];
r=J(:,1)==4;
J=J+[x y x y x y x y;x y x y x y x y];
J(r,:)=[];

% Vertices list defined by connected image grid-points
VP=[J(:,1:4);J(:,5:8)];

% Make a Face list
F=[(1:size(J,1))' (size(J,1)+1:2*size(J,1))'];

% Remove dubplicate vertices
[VP,a,Ind]=unique(VP,'rows'); F=Ind(F);

% Vertices described by image grid-points to real
% linear Interpolated vertices
Vind1=VP(:,1)+(VP(:,2)-1)*size(Img,1);
Vind2=VP(:,3)+(VP(:,4)-1)*size(Img,1);
V1=abs(double(Img(Vind1))-double(isovalue));
V2=abs(double(Img(Vind2))-double(isovalue));
alpha=V2./(V1+V2);
Vx=VP(:,1).*alpha+VP(:,3).*(1-alpha);
Vy=VP(:,2).*alpha+VP(:,4).*(1-alpha);
V=[Vx Vy];

```