Skip to content
Snippets Groups Projects
Commit 22cb2c55 authored by Patrick M. Jensen's avatar Patrick M. Jensen
Browse files

Add MATLAB demo scripts

parent ad21a67b
Branches
No related tags found
No related merge requests found
function V = EnsureTransversal(V,T,t,tolfac)
% ENSURETRANSVERSAL Ensure mesh is transversal with given time plane
% V = EnsureTransversal(V,T,t)
% V = EnsureTransversal(V,T,t,tolfac)
%
% INPUT
% V - N x 4 array with vertex coordinates
% T - M x 4 array with tet. indices
% t - Time coordinate for time plane
% tolfac - (Optional, default: 10) Tolerance factor.
%
% OUTPUT
% V - N x 4 array with corrected vertex coordinates
%
% See also: MESHTIMESLICE
%
% Patrick M. Jensen, 2019, Technical University of Denmark
if nargin < 4, tolfac = 10; end
tol = tolfac*eps(class(V));
for ti = 1:size(T,1)
TIdx = T(ti,:);
Signs = RobustSign(V(TIdx,4) - t, tol);
nzer = nnz(Signs == 0);
% Check for problematic intersections
if nzer == 4
% Intersection is entire tet
V(TIdx,4) = V(TIdx,4) + 2*tolfac*eps(V(TIdx(1),4));
elseif nzer == 2
% Intersection is a line segment or triangle. If triangle, then the
% intersecting edge might be part of a non-manifold region, so we
% perturb the endpoints out anyway
Idx = find(Signs == 0);
offset = 2*tolfac*eps(V(TIdx(Idx),4));
V(TIdx(Idx),4) = V(TIdx(Idx),4) + offset;
elseif nzer == 1
% Intersection is a point or a triangle. If triangle, then
% neighborhood might have 'negative' curvature, so we perturb the
% vertex out anyway
Idx = find(Signs == 0);
offset = 2*tolfac*eps(V(TIdx(Idx),4));
V(TIdx(Idx),4) = V(TIdx(Idx),4) + offset;
end
end
\ No newline at end of file
function [F,V] = Make600Cell(Center,rs,rt,nsubdiv)
% MAKE600CELL Make 600 cell tetrahedron mesh
% [F,V] = Make600Cell()
% [F,V] = Make600Cell(Center)
% [F,V] = Make600Cell(Center,rs)
% [F,V] = Make600Cell(Center,rs,rt)
% [F,V] = Make600Cell(Center,rs,rt,nsubdiv)
%
% Patrick M. Jensen, 2019, Technical University of Denmark
if nargin < 1
Center = [0 0 0 0];
end
if nargin < 2
rs = 1;
end
if nargin < 3
rt = 1;
end
if nargin < 4
nsubdiv = 0;
end
if iscolumn(Center)
Center = Center';
end
V = [];
for i = 0:2^4-1
Signs = 2*bitget(i,1:4) - 1;
V = [V; 0.5*ones(1,4).*Signs]; %#ok<AGROW>
end
for i = 1:4
Pi = [0 0 0 0];
Pi(i) = 1;
V = [V; Pi; -Pi];
end
phi = 0.5*(1 + sqrt(5));
W = 0.5*[phi,1,1/phi,0];
PermIdxs = perms(1:4);
for i = 1:size(PermIdxs,1)
Idx = PermIdxs(i,:);
% Check if permutation is even
M = zeros(4);
for k = 1:4
M(k,Idx(k)) = 1;
end
if det(M) < 0
% Odd permuation so skip
continue;
end
for j = 0:2^4-1
Signs = 2*bitget(j,1:4) - 1;
if prod(Signs) > 0
V = [V; W(PermIdxs(i,:)).*Signs]; %#ok<AGROW>
end
end
end
F = convhulln(V);
% Subdivide
for si = 1:nsubdiv
for i = 1:size(F,1)
P1 = V(F(i,1),:);
P2 = V(F(i,2),:);
P3 = V(F(i,3),:);
P4 = V(F(i,4),:);
Vn1 = 0.5*(P2 + P1);
Vn2 = 0.5*(P3 + P1);
Vn3 = 0.5*(P4 + P1);
Vn4 = 0.5*(P2 + P3);
Vn5 = 0.5*(P2 + P4);
Vn6 = 0.5*(P3 + P4);
V = [V;
Vn1/norm(Vn1);
Vn2/norm(Vn2);
Vn3/norm(Vn3);
Vn4/norm(Vn4);
Vn5/norm(Vn5);
Vn6/norm(Vn6)]; %#ok<AGROW>
end
V = unique(V,'rows');
F = convhulln(V);
end
% Fix orientation so all normals point 'outwards'
for i = 1:size(F,1)
% Get tet points
P1 = V(F(i,1),:);
P2 = V(F(i,2),:);
P3 = V(F(i,3),:);
P4 = V(F(i,4),:);
% Get barycenter
PM = 0.25*(P1 + P2 + P3 + P4);
% Compute edge vectors
V1 = P2 - P1;
V2 = P3 - P1;
V3 = P4 - P1;
% Compute normal vector, see: https://math.stackexchange.com/a/2371039
W = [V1; V2; V3];
N = [-det(W(:,[2 3 4])) det(W(:,[1 3 4])) -det(W(:,[1 2 4])) det(W(:,[1 2 3]))];
if dot(PM,N) < 0
% Normal points inwards so swap tet. points
tmp = F(i,3);
F(i,3) = F(i,4);
F(i,4) = tmp;
end
end
V(:,1:3) = V(:,1:3) * rs;
V(:,4) = V(:,4) * rt;
V = V + Center;
\ No newline at end of file
function [Vs,Fs,TIs,VNs] = MeshTimeSlice(V,T,t,trimesh,VN)
% MESHTIMESLICE Find 3D cross section of 4D mesh at a given time
% [Vs,Fs,TIs,VNs] = MeshTimeSlice(V,T,t)
% [Vs,Fs,TIs,VNs] = MeshTimeSlice(V,T,t,trimesh)
% [Vs,Fs,TIs,VNs] = MeshTimeSlice(V,T,t,trimesh,VN)
%
% INPUT
% V - N x 4 array with vertex coodinates
% T - M x 4 array with tet. indices
% t - Time coordinate of time slice
% trimesh - (Optional, default: false) Ensure output is a triangle mesh
% VN - (Optional) N x 4 array with vertex normals. Only needed if
% VNs output is assigned.
%
% OUTPUT
% Vs - Ns x 3 array with vertex corrdinates for intersection
% Fs - Ms x (3 or 4) array with polygon indices for intersection
% TIs - Array with indices of intersecting tetrahedra
% VNs - Ns x 3 array with projected vertex normals.
%
% Patrick M. Jensen, 2019, Technical University of Denmark
if nargin < 4
trimesh = false;
end
compNorms = nargout > 3;
nvert = size(V,1);
ntet = size(T,1);
PEMap = zeros(nvert,nvert);
try
tol = 10*eps(class(V));
catch
tol = 0;
end
Vs = [];
Fs = [];
TIs = [];
if compNorms
VNs = [];
end
for ti = 1:ntet
% Check for intersection with time slice
P = V(T(ti,:),:);
Signs = RobustSign(P(:,4) - t, tol);
npos = nnz(Signs > 0);
nneg = nnz(Signs < 0);
nzer = nnz(Signs == 0);
if nzer == 4
% Whole tet is in time slice so add all points
i1 = AddPt(T(ti,1));
i2 = AddPt(T(ti,2));
i3 = AddPt(T(ti,3));
i4 = AddPt(T(ti,4));
TIs = [TIs ti]; %#ok<AGROW>
if trimesh
Fs = [Fs;
i1 i2 i3;
i1 i2 i4;
i1 i3 i4;
i2 i3 i4]; %#ok<AGROW>
else
Fs = [Fs;
i1 i2 i3 -1;
i1 i2 i4 -1;
i1 i3 i4 -1;
i2 i3 i4 -1]; %#ok<AGROW>
end
elseif nzer == 0 && npos == 2 && nneg == 2
% Intersection is a quad
IdxPos = find(Signs > 0);
IdxNeg = find(Signs < 0);
[P1,i1] = GetInt(T(ti,IdxPos(1)),T(ti,IdxNeg(1)));
[P2,i2] = GetInt(T(ti,IdxPos(1)),T(ti,IdxNeg(2)));
[P3,i3] = GetInt(T(ti,IdxPos(2)),T(ti,IdxNeg(1)));
[P4,i4] = GetInt(T(ti,IdxPos(2)),T(ti,IdxNeg(2)));
Cen = 0.25*(P1 + P2 + P3 + P4);
N = cross(P1 - Cen,P2 - Cen);
% Manual sort
if CWLess(P2,P1,Cen,N)
tmp = i1;
i1 = i2;
i2 = tmp;
end
if CWLess(P3,P4,Cen,N)
tmp = i3;
i3 = i4;
i4 = tmp;
end
TIs = [TIs ti]; %#ok<AGROW>
if trimesh
Fs = [Fs; i1,i2,i3; i2,i3,i4]; %#ok<AGROW>
else
Fs = [Fs; i1,i2,i4,i3]; %#ok<AGROW>
end
elseif nzer == 3
% Triangle face intersects slice
Idx = 1:4;
Idx(Signs ~= 0) = [];
i1 = AddPt(T(ti,Idx(1)));
i2 = AddPt(T(ti,Idx(2)));
i3 = AddPt(T(ti,Idx(3)));
TIs = [TIs ti]; %#ok<AGROW>
if trimesh
Fs = [Fs; i1 i2 i3]; %#ok<AGROW>
else
Fs = [Fs; i1 i2 i3 -1]; %#ok<AGROW>
end
elseif npos == 1 || nneg == 1
% Intersection is a triangle
if npos == 1
IdxA = find(Signs > 0);
IdxB = find(Signs <= 0);
else
IdxB = find(Signs >= 0);
IdxA = find(Signs < 0);
end
[~,i1] = GetInt(T(ti,IdxA(1)),T(ti,IdxB(1)));
[~,i2] = GetInt(T(ti,IdxA(1)),T(ti,IdxB(2)));
[~,i3] = GetInt(T(ti,IdxA(1)),T(ti,IdxB(3)));
TIs = [TIs ti]; %#ok<AGROW>
if trimesh
Fs = [Fs; i1 i2 i3]; %#ok<AGROW>
else
Fs = [Fs; i1 i2 i3 -1]; %#ok<AGROW>
end
elseif nzer == 2 && (npos == 2 || nneg == 2)
% Intersection is a line segment
Idx = find(Signs == 0);
i1 = AddPt(T(ti,Idx(1)));
i2 = AddPt(T(ti,Idx(2)));
TIs = [TIs ti]; %#ok<AGROW>
if trimesh
Fs = [Fs; i1 i2 -1]; %#ok<AGROW>
else
Fs = [Fs; i1 i2 -1 -1]; %#ok<AGROW>
end
elseif nzer == 1 && (npos == 3 || nneg == 3)
% Intersection is a point
TIs = [TIs ti]; %#ok<AGROW>
Idx = find(Signs == 0);
AddPt(T(ti,Idx(1)));
% No edges to add
end
% Else there is no intersection
end
% Remove duplicate faces
[~,FIdx] = unique(sort(Fs,2),'rows');
Fs = Fs(FIdx,:);
Fs(Fs == -1) = nan;
function i = AddPt(vi)
if PEMap(vi,vi) == 0
Vs = [Vs; V(vi,1:3)];
if compNorms
VNs = [VNs; VN(vi,:)];
end
i = size(Vs,1);
PEMap(vi,vi) = i;
else
i = PEMap(vi,vi);
end
end
function [P,i] = GetInt(i1,i2)
if V(i1,4) == t
if PEMap(i1,i1) == 0
P = V(i1,1:3);
Vs = [Vs; P];
if compNorms
VNs = [VNs; VN(i1,:)];
end
i = size(Vs,1);
PEMap(i1,i1) = i;
else
i = PEMap(i1,i1);
P = Vs(i,:);
end
elseif V(i2,4) == t
if PEMap(i2,i2) == 0
P = V(i2,1:3);
Vs = [Vs; P];
if compNorms
VNs = [VNs; VN(i2,:)];
end
i = size(Vs,1);
PEMap(i2,i2) = i;
else
i = PEMap(i2,i2);
P = Vs(i,:);
end
else
if PEMap(i1,i2) == 0
a = (V(i1,4) - t)/(V(i1,4) - V(i2,4));
P = (1 - a)*V(i1,1:3) + a*V(i2,1:3);
Vs = [Vs; P];
if compNorms
N = (1 - a)*VN(i1,:) + a*VN(i2,:);
VNs = [VNs; N];
end
i = size(Vs,1);
PEMap(i1,i2) = i;
PEMap(i2,i1) = i;
else
i = PEMap(i1,i2);
P = Vs(i,:);
end
end
end
end
function l = CWLess(P1,P2,Cen,N)
l = dot(N, cross(P1 - Cen, P2 - Cen)) > 0;
end
\ No newline at end of file
function S = RobustSign(X,tol)
S = zeros(size(X),'like',X);
S(X > tol) = 1;
S(X < -tol) = -1;
end
\ No newline at end of file
%% Add folder with mex files
% You need to set this to your build dir
addpath surfseg/build/x64-Release/surfseg/matlab
%% Make volume containing three overlapping balls
V1 = false(100,100,100);
V1(35,50,35) = true;
V1 = single(bwdist(V1) <= 30);
V2 = false(100,100,100);
V2(70,50,35) = true;
V2 = single(bwdist(V2) < 30);
V3 = false(100,100,100);
V3(50,50,70) = true;
V3 = single(bwdist(V3) < 30);
V = max(V1,V2);
V = max(V,V3);
clear V1 V2 V3;
%% Display volume
figure(1); clf;
patch(isosurface(V,0.5),'EdgeColor','none','FaceColor',[0.8,0.8,0.8]);
axis equal;
axis([1 size(V,2) 1 size(V,1) 1 size(V,3)]);
view(45,20);
lighting phong
camlight right
material dull
%% Segment balls
% Define (very) simple region cost volume
Cost = (1 - V) - V;
% Segmentation parameters
Cen = single([35 50 35; % Mesh centers
70 50 35;
50 50 70]);
R = single([10,10,10]); % Mesh radii
nsub = 3; % Subdiv. lvl. for subdiv. icosahedron
nsamples = 50; % Number of samples
step = 1; % Step size for samples
smoothness = 10; % Smoothness parameter
costtype = 1; % Cost type (1 = region costs)
Conn = ones(3) - eye(3); % Mesh connectivity
Conn = logical(sparse(Conn));
[Fcs,Vtx] = mex_surfcut_planesep_qpbo(Cost,Cen,R,nsub,nsamples,step,...
smoothness,Conn,costtype);
%% Display segmentation
figure(2); clf;
patch(isosurface(V,0.5),'EdgeColor','none','FaceColor',[0.8,0.8,0.8],...
'FaceAlpha',0.6);
hold on;
Colors = {'r','g','b'};
for i = 1:length(Fcs)
% Note that vertices need to be adjusted
patch('Faces',Fcs{i},'Vertices',Vtx{i}(:,[2 1 3])+1,...
'EdgeColor',Colors{i},'FaceColor','none');
end
axis equal;
axis([1 size(V,1) 1 size(V,1) 1 size(V,1)]);
view(45,20);
lighting phong
camlight right
material dull
%% Add folder with mex files
% You need to set this to your build dir
addpath surfseg/build/x64-Release/surfseg/matlab
%% Make volume containing a single ball
V = false(100,100,100);
V(50,50,50) = true;
V = single(bwdist(V) <= 30);
%% Display volume
figure(1); clf;
patch(isosurface(V,0.5),'EdgeColor','none','FaceColor',[0.8,0.8,0.8]);
axis equal;
axis([1 size(V,1) 1 size(V,1) 1 size(V,1)]);
view(45,20);
lighting phong
camlight right
material dull
%% Segment ball
% Define (very) simple region cost volume
Cost = (1 - V) - V;
% Segmentation parameters
Cen = single([50 50 50]); % Mesh center
r = single(10); % Mesh radius
nsub = 3; % Subdivision level for subdivided icosahedron
nsamples = 50; % Number of samples
step = 1; % Step size for samples
smoothness = 10; % Smoothness parameter
costtype = 1; % Cost type (1 = region costs)
[Fcs,Vtx] = mex_surfcut(Cost,Cen,r,nsub,nsamples,step,smoothness,costtype);
Fcs = Fcs{1};
Vtx = Vtx{1};
Vtx = Vtx(:,[2 1 3]) + 1; % Need to adjust vertices
%% Display segmentation
figure(2); clf;
patch(isosurface(V,0.5),'EdgeColor','none','FaceColor',[0.8,0.8,0.8],...
'FaceAlpha',0.6);
hold on;
patch('Faces',Fcs,'Vertices',Vtx,'EdgeColor','r','FaceColor','none');
axis equal;
axis([1 size(V,2) 1 size(V,1) 1 size(V,3)]);
view(45,20);
lighting phong
camlight right
material dull
%% Add folder with mex files
% You need to set this to your build dir
addpath surfseg/build/x64-Release/surfseg/matlab
%% Make volumetric time series with two growing balls
% This may take a few seconds...
V = ones(200,100,100,100,'single');
for t = 1:100
V1 = false(200,100,100);
V1(125,50,50) = true;
V1 = single(bwdist(V1) <= 5 + round(t*(40/100)));
V2 = false(200,100,100);
V2(75,50,50) = true;
V2 = single(bwdist(V2) <= 5 + round(t*(40/100)));
V(:,:,:,t) = max(V1,V2);
end
%% Display volumetric time series
for t = 1:5:size(V,4)
cla;
patch(isosurface(V(:,:,:,t),0.5),...
'EdgeColor','none','FaceColor',[0.8,0.8,0.8]);
axis equal;
axis([1 size(V,2) 1 size(V,1) 1 size(V,3)]);
title(sprintf('t = %d',t));
view(80,10);
lighting phong
camlight right
material dull
pause(0.5);
end
%% Compute (very) simple region cost volume
Cost = (1 - V) - V;
%% Segment balls using 600-cell as initialization
% Segmentation parameters
Cen = [100,50,50,95]; % Mesh center in 4D space
rs = 20; % Spatial radius
rt = 5; % Temporal radius
nsub = 2; % Subdivision level
nsamples = 100; % Number of samples
step = 1; % Step size for samples
smoothness = 40; % Smoothness parameter
costtype = 1; % Cost type (1 = region costs)
bend = false; % Experimental parameter; always set this to false
[FcsI,VtxI] = Make600Cell(Cen,rs,rt,nsub); % Init mesh
[Fcs,Vtx] = mex_surfcut_4d(Cost,int32(FcsI),single(VtxI),nsamples,step,...
smoothness,costtype,bend);
%% Display segmentation
for t = 1:5:size(V,4)
cla;
% Display time slice
patch(isosurface(V(:,:,:,t),0.5),...
'EdgeColor','none','FaceColor',[0.8,0.8,0.8],'FaceAlpha',0.6);
hold on;
% Compute and display cross section of segmentation hypersurface
VtxT = EnsureTransversal(Vtx,Fcs,t);
[Vs,Fs] = MeshTimeSlice(VtxT,Fcs,t);
Vs = Vs(:,[2 1 3]) + 1; % Adjust vertices
patch('Faces',Fs,'Vertices',Vs,'EdgeColor','r','FaceColor','none');
axis equal;
axis([1 size(V,2) 1 size(V,1) 1 size(V,3)]);
title(sprintf('t = %d',t));
view(80,10);
lighting phong
camlight right
material dull
pause(0.5);
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment