www.gusucode.com > 基于matlab编程深度图象数据用平滑算子和差分方法求曲面的 > 基于matlab编程深度图象数据用平滑算子和差分方法求曲面的高斯曲率和平均曲率源码程序/code/main.m

    clear all;clc 
 
%实验一:半球深度图像的生成 
a=256;b=256;c=0;R=200; 
%本实验中取(a,b,c)=(256,256,0),R=10,即求新坐标为(256,256,0),半径为200; 
for i=1:512 
    for j=1:512 
        if R^2-((i-a)^2+(j-b)^2)>=0 
        f(i,j)=sqrt(R^2-((i-a)^2+(j-b)^2)); 
        else 
            f(i,j)=0; 
        end 
    end 
end 
imshow(uint8(f));title('深度图像');%显示半球深度图像 
 
 
%实验二:曲率图像的计算 
s0=[1,6,15,20,15,6,1]/64; 
d0=[1,1,1,1,1,1,1]/7; 
d1=[-3,-2,-1,0,1,2,3]/28; 
d2=[5,0,-3,-4,-3,0,5]/84; 
 
s=s0'*s0; 
Du=d0'*d1; 
Dv=d1'*d0; 
Duu=d0'*d2; 
Dvv=d2'*d0; 
Duv=d1'*d1; 
 
fu=conv2(conv2(Du,s),f); 
fv=conv2(conv2(Dv,s),f); 
fuu=conv2(conv2(Duu,s),f); 
fvv=conv2(conv2(Dvv,s),f); 
fuv=conv2(conv2(Duv,s),f); 
 
K=(fuu.*fvv-fuv.^2)./(1+fu.^2+fv.^2).^2; 
H=0.5*((1+fv.^2).*fuu+(1+fu.^2).*fvv-2*fu.*fv.*fuv)./(1+fu.^2+fv.^2).^1.5; 
figure;imshow(uint8(abs(100000*H)));title('平均曲率图像');%此处可试探性的找到合适的系数 
figure;imshow(uint8(abs(40000000*K)));title('高斯曲率图像');%显示高斯曲率图像 
 
 
%实验三：曲面的分类 
%在0~255之间进行均匀标号，以扩大不同标号间的反差,便于区分显示 
[M,N]=size(H);%上面得到的H,K应该大小相同 
for i=1:M 
    for j=1:N 
        if H(i,j)<0&K(i,j)>0 %峰 
            X(i,j)=0; 
        elseif H(i,j)<0&K(i,j)==0 %无 
            X(i,j)=30; 
        elseif H(i,j)<0&K(i,j)<0 %凹 
            X(i,j)=60; 
        elseif H(i,j)==0&K(i,j)>0 %脊 
            X(i,j)=90; 
        elseif H(i,j)==0&K(i,j)==0 %平面 
            X(i,j)=120; 
        elseif H(i,j)==0&K(i,j)<0 %谷 
            X(i,j)=150; 
        elseif H(i,j)>0&K(i,j)>0 %鞍脊 
            X(i,j)=180; 
        elseif H(i,j)>0&K(i,j)==0 %最小表面 
            X(i,j)=210; 
        elseif H(i,j)>0&K(i,j)<0 %鞍谷 
            X(i,j)=240; 
        end 
    end 
end 
figure;imshow(uint8(X));title('标号图像');%显示标号图像 
figure;imhist(uint8(X));title('标号分布图');