MATLAB实现Catmull

MATLAB实现Catmull

终于调试好了Catmull-Clark细分（CC细分）的全部程序，将之前只适用于封闭四边形网格的程序进行了完善

主要一段代码来自于三维网格细分算法（Catmull-Clark subdivision & Loop subdivision）附源码，这个博主的很多篇博文都写的非常好，但是经常丢三落四的，像在这篇博文中他就用到了函数outline.m用来计算网格的边界，但是博主却没有给出outline函数，我自己重新编写了这个函数，并且能够成功执行^^，现在我贴出完整代码

function [VV, FF, S] = CCSubdivision(V, F, iter)  % Catmull_Clark subdivision  if ~exist('iter','var')  iter = 1;  end  VV = V;  FF = F;  for i = 1:iter   nv = size(VV,1);  nf = size(FF,1);   O = outline(FF);  original = 1:nv;  boundary = O(:,1)';  interior = original(~ismember(original, boundary));  no = length(original);  nb = length(boundary);  ni = length(interior);  %% Sv  Etmp = sort([FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)],2);  [E, ~, idx] = unique(Etmp, 'rows');  Aeven = sparse([E(:,1) E(:,2)], [E(:,2) E(:,1)], 1, no, no);  Aodd = sparse([FF(:,1) FF(:,2)], [FF(:,3) FF(:,4)], 1, no, no);  Aodd = Aodd + Aodd';  val_even = sum(Aeven,2);  beta = 3./(2*val_even);  val_odd = sum(Aodd,2);  gamma = 1./(4*val_odd);  alpha = 1 - beta - gamma;  Sv = sparse(no,no);  Sv(interior,:) = ...  sparse(1:ni, interior, alpha(interior), ni, no) + ...  bsxfun(@times, Aeven(interior,:), beta(interior)./val_even(interior)) + ...  bsxfun(@times, Aodd(interior,:), gamma(interior)./val_odd(interior));  Sboundary = ...  sparse([O(:,1);O(:,2)],[O(:,2);O(:,1)],1/8,no,no) + ...  sparse([O(:,1);O(:,2)],[O(:,1);O(:,2)],3/8,no,no);  Sv(boundary,:) = Sboundary(boundary,:);  %% Sf  Sf = 1/4 .* sparse(repmat((1:nf)',1 ,4), FF, 1);  i0 = no + (1:nf)';  %% Se  flaps = sparse([idx;idx], ...  [FF(:,3) FF(:,4);FF(:,4) FF(:,1);FF(:,1) FF(:,2);FF(:,2) FF(:,3)], ...  1);  onboundary = (sum(flaps,2) == 2);  flaps(onboundary,:) = 0;  ne = size(E,1);  Se = sparse( ...  [1:ne 1:ne]', ...  [E(:,1); E(:,2)], ...  [onboundary;onboundary].*1/2 + ~[onboundary;onboundary].*3/8, ...  ne, ...  no) + ...  flaps*1/16;  %% new faces & new vertices  i1 = no +   nf + (1:nf)';  i2 = no + 2*nf + (1:nf)';  i3 = no + 3*nf + (1:nf)';  i4 = no + 4*nf + (1:nf)';  FFtmp = [i0 i4 FF(:,1) i1; ...  i0 i1 FF(:,2) i2; ...  i0 i2 FF(:,3) i3; ...  i0 i3 FF(:,4) i4];  reidx = [(1:no)'; no+(1:nf)'; no+nf+idx];  FF = reidx(FFtmp);  S = [Sv; Sf; Se];  VV = S*VV;  end  end  

function out = outline( FF )
%OUTLINE Summary of this function goes here
%   Detailed explanation goes here
Etmp = sort([FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)],2);
[~, ~, idx] = unique(Etmp, 'rows');oriEtmp = [FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)];
hh=sortrows([oriEtmp,idx],3);x2=diff(sortrows(idx));
vector = all(x2==0, 2);index1=find(vector);
index2=index1+1;
index=[index1;index2];hh(index,:)=[];
out=hh(:, 1:2);end

[V,F]=obj__read('six.obj');
V=V';F=F';
iter=4;
[VV, FF] = CCSubdivision(V, F, iter);
%[VV, FF] = CCsub(V, F, iter);
obj_write('six1.obj',VV',FF');[V,F]=obj__read('torus.obj');
V=V';F=F';
iter=4;
[VV, FF] = CCSubdivision(V, F, iter);
%[VV, FF] = CCsub(V, F, iter);
obj_write('torus1.obj',VV',FF');