計算數學入門系列 - 距離函數（三）

 

一般來說，找出曲線的距離函數的表示方式並不簡單。舉個例來說我們如果想表示一個橢圓形， 最簡單的隱性表示方式就是找 \(f\left(x,y\right)=\frac{x^2}{a^2}+\frac{y^2}{b^2}-1\)。 但是計算一下，就可以知道他並不是一個距離函數。 甚至乎，\(f\left(x,y\right)=\sqrt{\frac{x^2}{a^2}+\frac{y^2}{b^2}}-1\)畫起來有點像，可是仍然不是我們想找的距離函數。

Figure 6. "ellipse1.jpg". An implicit representation of an ellipse.

a = 0.25; b = 0.7;

f=sqrt(x.^2/a^2+y.^2/b^2)-1;

contour(x,y,g,30); axis equal; axis([-1.1 1.1 -1.1 1.1])

hold on

contour(x,y,g,[0 0],'k--')

colormap('jet')

print -djpeg ellipse1.jpg

在Figure 6 可以看到，我們要找的橢圓形附近，y方向的等高線相對比x方向距離比較遠。要找出橢圓形的距離函數，我們就必須 要從每一個meshpoint找出在橢圓上距離最短的距離。

Figure 7. "ellipse2.jpg". The signed distance function representation of an ellipse.

[m,n]=size(x);

f = @(theta,x0,a,b) (x0(1)-a*cos(theta)).^2+(x0(2)-b*sin(theta)).^2;

x0=zeros(1,2);

for i=1:m

for j=1:n

   x0(1)=x(i,j);

   x0(2)=y(i,j);

   theta0 = atan2(x0(2),x0(1));

   theta_s = fminsearch(@(theta) f(theta,x0,a,b),theta0);

   sdist(i,j) = sign(sdist(i,j))*sqrt(f(theta_s,x0,a,b));

end

end 

contour(x,y,sdist,30); axis equal; axis([-1.1 1.1 -1.1 1.1])

hold on

contour(x,y,sdist,[0 0],'k--')

colormap('jet')

print -djpeg ellipse2.jpg

這段程式裏面，函數f 想代表某一點x0 跟在橢圓形上面的一點 (x(theta),y(theta)) = (a cos(theta), b sin(theta)) 距離的平方。要找出x0 這個位置的距離函數值，我們就需要改變theta 的數值去優化函數f。 而這個優化問題，我們就直接使用MATLAB內置的程式fminsearch 去解決。一些比橢圓形更複雜的曲線，我們還是可以使用差不多的方法把距離函數找出來。

Figure 8. "star5.jpg". The signed distance function representation of a 5-folded star-shape.

x=[-1:0.005:1];
[x,y]=meshgrid(x,x);
[m,n]=size(x);
a=0.5; b=0.1;
f = @(theta,x0,a,b) (x0(1)-a*(1+b*cos(5*theta))*cos(theta)).^2+...
   (x0(2)-a*(1+b*cos(5*theta))*sin(theta)).^2;
sdist=zeros(m,n);
x0=zeros(1,2);
for i=1:m
for j=1:n
   x0(1)=x(i,j);
   x0(2)=y(i,j);
   theta0 = atan2(x0(2),x0(1));
   theta_s = fminsearch(@(theta) f(theta,x0,a,b),theta0);
   sdist(i,j) = sqrt(x0(1)^2+x0(2)^2)-a*(1+b*cos(5*theta0));
   sdist(i,j) = sign(sdist(i,j))*sqrt(f(theta_s,x0,a,b));
end
end
contour(x,y,sdist,30); axis equal; axis([-1 1 -1 1])
hold on
contour(x,y,sdist,[0 0],'k--')
colormap('jet')
print -djpeg star5.jpg