網格測地線算法（Geodesics in Heat）附源碼

　　測地線又稱為大地線，可以定義為空間曲面上兩點的局部最短路徑。測地線具有廣泛的應用，例如在工業上測地線最短的性質就意味着最優最省，在航海和航空中，輪船和飛機的運行路線就是測地線。[Crane et al. 2013]提出了利用熱運動方程來計算網格測地線的方法，可以想象一下，當一根燙的針尖接觸到曲面上的一點時，熱量會隨着時間的推移而擴散，測地距離因此可以和熱運動相聯系。具體算法過程如下圖所示：

　　第一步：熱運動方程用來描述熱的傳播狀態：

　　將熱運動方程離散化並整理后得到：

其中：id為單位矩陣，t為時間間隔，Δ為離散Laplacian算子，u_t為t時刻的熱狀態，u₀為初始時刻的熱狀態。

　　第二步：第一步計算得到的熱梯度方向與測地距離的梯度方向相同，由Eikonal方程知道測地距離的梯度為單位向量，於是通過歸一化熱梯度我們得到測地距離的梯度：

　　第三步：得到測地距離的梯度之后，測地線問題即變為求解以下式子：

　　根據變分法，上式最小化即求解泊松方程：

其中：Φ即為網格上頂點距離源點的測地距離。

 

function [D] = geodesics_in_heat(V, F, src)   
    % choose time step
    c = 5;
    t = c * mean(doublearea(V, F))/2;

    %% Step 1: Integrate the heat flow for some fixed time t
    L = cotmatrix(V, F);
    M = massmatrix(V, F, 'barycentric');
 
    nV = size(V, 1);
    u0 = zeros(nV, 1);
    u0(src) = 1;
    
    A = M - t*L;
    B = M*u0;

    nsrc = length(src);
    % 1.1 dirichlet condition
    hole = Cal_Boundary(F);
    if isempty(hole)
        boundary = [];
    else
        boundary = hole.boundary.edge(:,1)';
    end

    b = setdiff(boundary, src);
    nb = [b, src];
    Acons = sparse([1:length(nb)], nb, ones(1,length(nb)), length(nb), nV);
    Bcons = [zeros(length(b), 1); ones(nsrc, 1)];

    % 硬約束        
    r = setdiff([1:nV], nb);
    uD = [A(r,:);Acons]\[B(r,:);Bcons];

    % 1.2 neumann condition
    Acons = sparse([1:nsrc], src, ones(1,nsrc), nsrc, nV);
    Bcons = ones(nsrc, 1);

    % 硬約束  
    r = setdiff([1:nV], src);
    uN = [A(r,:);Acons]\[B(r,:);Bcons];
    
    % averaged boundary condition
    u = 0.5*(uN + uD);
    
    %% Step 2: Evaluate the vector field X
    G = grad(V, F); % nF*3 by nV matrix(梯度算子 - 所有三角片頂點基函數)
    grad_u = reshape(G*u, size(F,1), 3); % nF by nV matrix(所有三角片中u的梯度)
    grad_u_norm = sqrt(sum(grad_u.^2, 2));
    normalized_grad_u = bsxfun(@rdivide, grad_u, grad_u_norm+eps);
    X = -normalized_grad_u;
    
    %% Step 3: Solve the Poisson equation
    Div = div(V, F); % 散度算子
    div_X = Div*X(:); % #nV by #nF*3

    Lcons = sparse([1:nsrc], src, ones(1,nsrc), nsrc, nV);
    div_Xcons = zeros(nsrc, 1);

    % 硬約束
    r = setdiff([1:nV], src);
    D = [L(r,:);Lcons]\[div_X(r,:);div_Xcons];
    
    D = D - min(D);
end

 

效果：

本文為原創，轉載請注明出處：http://www.cnblogs.com/shushen。

 

 

參考文獻：

[1] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. 2013. Geodesics in heat: A new approach to computing distance based on heat flow. ACM Trans. Graph. 32, 5, Article 152 (October 2013), 11 pages.