Grid geodesic algorithm (geodesics in Heat) attached source code

Geodesic, also known as the geodetic line, can be defined as the local shortest path of two points on a space surface. Geodesic has a wide range of applications, for example, the shortest property of the industrial geodesic means that the best and most provincial, in navigation and aviation, ships and aircraft operating route is geodesic. [Crane et al. 2013] A method of calculating the grid geodesic using the thermal equation of motion is proposed, and it can be imagined that when a hot tip touches a point on the surface, the heat will spread over time and the geodesic distance can be associated with the thermal movement. The specific algorithm process is as follows:

The first step: the thermal equation of motion is used to describe the propagation state of heat:

The equations of thermal motion are discretized and collated to obtain:

Where: ID is the unit matrix, T is the time interval, δ is a discrete Laplacian operator, UT is the thermal state of T moment, and u0 is the thermal state at the initial moment.

The second step: The first step calculated the thermal gradient direction is the same as the gradient direction of geodesic distance, the eikonal equation knows the gradient of the geodesic distance is the unit vector, so through the normalized thermal gradient we get the ground distance gradient:

Step three: after getting the gradient of geodesic distance, the geodesic problem becomes the solution of the following equation:

On the basis of the Variational method, the Poisson equation is solved by minimizing the upper formula:

Where: Φ is the geodesic distance from the source point on the vertex of the grid.

function [D] =geodesics_in_heat (V, F, SRC)%Choose time Step c=5; T= c * Mean (Doublearea (V, F))/2; Percent Step1: Integrate the Heat flow forSomefixedTime T L=Cotmatrix (V, F); M= Massmatrix (V, F,'Barycentric'); NV= Size (V,1); U0= Zeros (NV,1); U0 (SRC)=1; A= m-t*m; B= m*u0; NSRC=length (SRC); %1.1Dirichlet condition Hole=cal_boundary (F); ifIsEmpty (hole) boundary= []; ElseBoundary= 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)]; %Hard constrained R= Setdiff ([1: NV], NB); UD=[A (R,:); Acons]\[b (R,:);    Bcons]; %1.2Neumann condition Acons= Sparse ([1: Nsrc], SRC, ones (1, NSRC), NSRC, NV); Bcons= Ones (NSRC,1); %Hard constrained R= Setdiff ([1: NV], SRC); UN=[A (R,:); Acons]\[b (R,:);        Bcons]; %averaged boundary condition U=0.5* (UN +UD); Percent Step2: Evaluate the vector field X G= Grad (V, F); % nf*3by NV Matrix (gradient operator-all triangular vertex base functions) Grad_u= Reshape (g*u, size (F,1),3); %NF by NV Matrix (gradient of u in all triangles) Grad_u_norm= sqrt (sum (grad_u.^2,2)); Normalized_grad_u= Bsxfun (@rdivide, Grad_u, grad_u_norm+EPS); X= -Normalized_grad_u; Percent Step3: Solve the Poisson equation Div= Div (V, F); %div_x of divergence operator= Div*x (:);% #nV by #nF *3lcons= Sparse ([1: Nsrc], SRC, ones (1, NSRC), NSRC, NV); Div_xcons= Zeros (NSRC,1); %Hard constrained R= Setdiff ([1: NV], SRC); D=[L (R,:);        Lcons]\[div_x (R,:);d Iv_xcons]; D= D-min (D); end

Effect:

Welcome everyone to discuss computer graphics algorithm problem, email: [Email protected]

This article is original, reprint please indicate source: Http://www.cnblogs.com/shushen.

Reference documents:

[1] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. Geodesics in heat:a New approach to computing distance based on heat flow. ACM Trans. Graph. 5, article (October), pages.

Grid geodesic algorithm (geodesics in Heat) attached source code