www.gusucode.com > 三维人体重构源码程序 > 三维人体重构源码程序/3L-IBS/IBSLevelSurf.m

    function IBSLevelSurf(P,pnt,RES)
%% ====================================================================
% Author: Mohammad Rouhani, Morpheo Team, INRIA Rhone Alpes, (2013)
% Email: mohammad.rouhani@inria.fr
% Title: IBS value calculation.
% Place of publication: Grenoble, France
% Available from: URL
% http://www.iis.ee.ic.ac.uk/~rouhani/mycodes/Geometric.zip
%====================================================================
% When using this software, PLEASE ACKNOWLEDGE the effort that went 
% into development BY REFERRING THE PAPER:
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 
% Rouhani M. and Sappa A.D., Implicit B-spline ?tting using the 3L 
% algorithm, IEEE Conference on on Image Processing (ICIP'11), 2011.
%::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 
%% ====================================================================
%About this code:
%Input: P: IBS coefficients; [x,y,z]: point coordinates.
%Output: f: IBS value at the given points; Buvw: IBS monomials.
ShowOtherLevels=0;
run_mine=0;
%surface plot:---------
%Xrange=-0:.1:10; Yrange=-0:.1:10;
box=[.1 .9 .1 .9 .1 .9];
box(1:2:5)=min(pnt)-.05; box(2:2:6)=max(pnt)+.05;
%box=[ .1 .85 .4 .6 .15 .85]
%Xrange=0:step:boundx; Yrange=0:step:boundy; Zrange=0:step:boundz;
step=1/RES;
Xrange=box(1):step:box(2); Yrange=box(3):step:box(4); Zrange=box(5):step:box(6);

%Xrange=-0:step:boundx; Yrange=-0-.2:step:boundy; Zrange=-0:step:boundz;
%Xrange=-.8:step:.8; Yrange=-.6:step:1; Zrange=.45:step:1;
[X,Y,Z] = meshgrid(Xrange,Yrange,Zrange);
tic;
F=BSpline3DF(P,X(:),Y(:),Z(:));
Time_REC=toc
W=reshape(F,size(X));
%Isosurface from MatLab:--------------------
%W=postprocess(W,X,Y,Z,x,y,z,5);

%Isosurface from MatLab:--------------------
%W=postprocess(W,X,Y,Z,x,y,z,5);
stop=0; Level=0.0; vw=2; [-58 48]; [90 0];
ch=[.5 .6 .8]; %ch=[.9 .8 .1];
while ~stop    
    hh=figure; set(hh,'Position',[900 200 700 800]); 
    hold on; axis equal off; view(vw); camlight;
    p = patch(isosurface(X, Y, Z, W, Level));
    isonormals(X,Y,Z,W, p); 
    set(p, 'FaceColor', ch, 'EdgeColor', 'none'); lighting phong;
    set(gca, 'LooseInset', [0,0,0,0]);axis equal off
    if ShowOtherLevels
        Level=input('Which iso-surface to show? (exit: 0) ');
    end
    if (Level==0); stop=1; end
end

%My own function without rendering:--------------------
if run_mine
    figure
    cntr=0;
    for i=2:size(W,1)-1
        levelNo(i)=cntr;
        if cntr==0
            firstInd=i;
        end
        for j=2:size(W,2)-1
            for k=2:size(W,3)-1
                if (W(i,j,k)*W(i-1,j,k)<0)||(W(i,j,k)*W(i,j-1,k)<0)||(W(i,j,k)*W(i,j,k-1)<0)
                    if (W(i,j,k-1)*W(i,j,k+1)<0)
                        cntr=cntr+1;
                        %scatter3(X(i,j,k),Y(i,j,k),Z(i,j,k),'r.')
                        x(cntr)=X(i,j,k); y(cntr)=Y(i,j,k); z(cntr)=Z(i,j,k);
                        %myPatch(c,step,x(cntr),y(cntr),z(cntr));
                    end
                end
            end
        end
    end
end
%box on;

function F=BSpline3DF(P,x,y,z)
%%
n=numel(x); N=size(P,1); %N=size(P,2); O=size(P,3); 
stepX=1/(N-3);
%MM=zeros(n,N^3);
F=zeros(n,1);

M=size(P,1); N=size(P,2); O=size(P,3);
stepX=1/(M-3); stepY=1/(N-3); stepZ=1/(O-3);

%Every point [x,y,z] is projected to [u,v,w]-space:
i=floor(x/stepX)+1; j=floor(y/stepY)+1; k=floor(z/stepZ)+1;
u=x/stepX-i+1; v=y/stepY-j+1; w=z/stepZ-k+1;
%Blending patches together
B=[-1 3 -3 1; 3 -6 3 0; -3 0 3 0; 1 4 1 0]/6;
bu=[u.^3 u.^2 u 1+0*u]*B; bv=[v.^3 v.^2 v 1+0*v]*B; bw=[w.^3 w.^2 w 1+0*w]*B;
for ii=0:3
    for jj=0:3
        for kk=0:3 
            %Matlab standard order: layer-by-layer, column-by...           
            Col=N^2*(k+kk-1)+N*(j+jj-1)+(i+ii); 
            F=F+P(Col).*bu(:,ii+1).*bv(:,jj+1).*bw(:,kk+1);
        end
    end
end

%removing zeros from the head:
function myTriangle(normal,P1,P2,P3)
hold on;
% normal=cross(P1-P2,P1-P3);
% if norm(normal)
%     normal=abs(normal)/norm(normal);
% else
%     normal=[0 0 0];
% end
col=abs(normal*[1 0 0]')*[1 0 0]
%normal*[1 0 0]'*[1 0 0]
%col=col/norm(col);
if isnan(col)
    col=[1 0 0];
end
X=[P1(1),P2(1),P3(1)]; Y=[P1(2),P2(2),P3(2)]; Z=[P1(3),P2(3),P3(3)];
patch(X,Y,Z,abs(col));



function myPatch(c,step,xi,yi,zi)
%first of all, we should find 4 points on the corners:
[f, fx,fy,fz,fxx,fxy,fxz,fyy,fyz,fzz]=Poly(deg,c,xi,yi,zi);
%I hope fz is nonzero!
DeltaZ1=(-fx/fz-fy/fz)*step/2;
DeltaZ2=(+fx/fz-fy/fz)*step/2;
swapxz=0;
if abs(DeltaZ1)+abs(DeltaZ2)>8*step
    DeltaX1=(-fz/fx-fy/fx)*step/2;
    DeltaX2=(+fz/fx-fy/fx)*step/2;
    %swap x & z:
    swapxz=1;
    tt=xi; xi=zi; zi=tt;
    tt=fx; fx=fz; fz=tt;
    DeltaZ1=DeltaX1; DeltaZ2=DeltaX2;    
end
P=[xi,yi,zi];
normal=[fx, fy, fz]; normal=normal/norm(normal);
%quiver3(P(1),P(2),P(3),n(1),n(2),n(3),'g','LineWidth',2);
P1=[xi+step/2, yi+step/2,zi+DeltaZ1];
P2=[xi-step/2, yi+step/2,zi+DeltaZ2];
P3=[xi-step/2, yi-step/2,zi-DeltaZ1];
P4=[xi+step/2, yi-step/2,zi-DeltaZ2];
if swapxz
    tt=P(1); P(1)=P(3); P(3)=tt;
    tt=P1(1); P1(1)=P1(3); P1(3)=tt;
    tt=P2(1); P2(1)=P2(3); P2(3)=tt;
    tt=P3(1); P3(1)=P3(3); P3(3)=tt;
    tt=P4(1); P4(1)=P4(3); P4(3)=tt;
end
%the extreme case:
%draw all 4 triangles:
myTriangle(normal,P,P1,P2);
myTriangle(normal,P,P2,P3);
myTriangle(normal,P,P3,P4);
myTriangle(normal,P,P4,P1);
% scatter3(P(1),P(2),P(3),120,'filled')
% scatter3(P1(1),P1(2),P1(3),'filled')
% scatter3(P2(1),P2(2),P2(3),'filled')
% scatter3(P3(1),P3(2),P3(3),'filled')
% scatter3(P4(1),P4(2),P4(3),'filled')