-
Notifications
You must be signed in to change notification settings - Fork 1
/
example5.m
96 lines (77 loc) · 4.82 KB
/
example5.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
% Simple application: drawing approximate piece-wise linear isodistance curves (isolines).
% In theory, isolines on a triangular mesh are piecewise second-order curves.
% Here, we draw piecewise linear approximation of these curves. The distance function
% of the mesh is sampled at the subdivision vertices within each triangle and
% represented as a linear function within each subdivision triangle.
% This particular method works well only for 'exact' geodesic method,
% because for the approximate methods the distance field is discontinuous
% on the edges of the mesh.
% Danil Kirsanov, 09/2007
subdivision_level = 5; %number of additional vertices per edge; each mesh triangle will have roughly speaking subdivision_level^2 number of subdivision triangles
number_of_isolines = 10;
global geodesic_library;
geodesic_library = 'geodesic_matlab_api'; %"release" is faster and "debug" does additional checks
rand('state', 0); %comment this statement if you want to produce random mesh every time
[vertices,faces] = create_flat_triangular_mesh(0.2, 0);
N = length(vertices);
% N = 300; %number of points in a mesh
% [vertices,faces] = create_hedgehog_mesh(N, 0.1); %create "noisy sphere" mesh; "vertices" contains 3D vertex coordinates; "faces" contains vertex id's for every triangle
hold off;
trisurf(faces,vertices(:,1),vertices(:,2),vertices(:,3), 'FaceColor', 'w', 'EdgeColor', 'k', 'FaceAlpha', 0.99); %plot the mesh
hold on;
mesh = geodesic_new_mesh(vertices,faces); %initilize new mesh
algorithm = geodesic_new_algorithm(mesh, 'exact'); %initialize new geodesic algorithm
source_points = {};
for vertex_id = [26,N-5, floor(N/2)-3, floor(N/2)+3]; %create source points
source_points{end+1} = geodesic_create_surface_point('vertex',vertex_id, vertices(vertex_id,:));
end;
geodesic_propagate(algorithm, source_points);
[sub_weights,sub_tri] = create_subdivision_pattern(subdivision_level); %coordinates of the subdivision triangles
[source_id, distances] = geodesic_distance_and_source(algorithm); %find distances for all vertices of the mesh
max_distance = max(distances);
isolines_step = max_distance/number_of_isolines;
isolines_distances = isolines_step/2 : isolines_step : max_distance;
[known_distances, tmp] = find(sub_weights > 0.999); %we already know the distances to the actual vertices of the mesh
unknown_distances = setdiff([1:size(sub_weights,1)], known_distances);
for face_id = 1:size(faces,1)
vertex_ids = faces(face_id,:);
sub_vertices = sub_weights*vertices(vertex_ids,:); %find coordinates of subdivision vertices
sub_distances = zeros(size(sub_vertices,1),1);
sub_distances(known_distances) = distances(vertex_ids); %no need to reevaluate distance finction for the vertices of the face
%for sub_index = 1:size(sub_vertices,1) %find distances at all subdivision vertices
for sub_index = unknown_distances %find distances at all internal subdivision vertices
q = geodesic_create_surface_point('face',face_id,sub_vertices(sub_index,:));
[source_id, d] = geodesic_distance_and_source(algorithm, q);
sub_distances(sub_index) = d;
end;
selected_distances = find(isolines_distances >= min(sub_distances) & ...
isolines_distances <= max(sub_distances));
selected_distances = isolines_distances(selected_distances);
for tri_index = 1:size(sub_tri,1) %intersection of the equidistance curve with every subdivision triangle is approximated as a straight line
for isodistance = selected_distances
p = [];
for i = 1:3
index1 = sub_tri(tri_index, i);
index2 = sub_tri(tri_index, mod(i,3) + 1);
d1 = sub_distances(index1);
d2 = sub_distances(index2);
if d1 == isodistance %special case when isodistance paths exactly through the vertex
p(end+1,:) = sub_vertices(index1,:);
end
if isodistance > min(d1,d2) & isodistance < max(d1,d2) %draw equidistant curve
a = (isodistance - d2)/(d1 - d2); % I believe that division by zero will never happen here
p(end+1,:) = a*sub_vertices(index1,:) + (1-a)*sub_vertices(index2,:);
end;
end
if size(p,1) == 2
plot3(p(:,1), p(:,2),p(:,3),'b', 'LineWidth', 2);
end
end
end;
end;
for i=1:length(source_points);
h = plot3(source_points{i}.x,source_points{i}.y,source_points{i}.z,'ro','MarkerSize',5, 'MarkerFaceColor', 'r');
end;
campos([0, 0, 5])
legend(h, 'sources');
geodesic_delete;