-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathmeshSnap.m
269 lines (219 loc) · 9.36 KB
/
meshSnap.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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
function [NewClusterStruct] = meshSnap(inputCluster,tol)
% MESHSNAP
%
% Function to merge nearby points within a set of meshes located within a
% distance of TOL.
%
% Inputs:
% - inputCluster: a struct containing Objects, containing the point cloud
% and mesh information. The Objects are separate entities; the function
% will however snap the mesh vertices by taking into account the
% global/merged state of all Objects within the inputCluster
% - tol: distance threshold between a point with its neighbours to be
8000
% considered as the same point
%
% Outputs:
% - vertices2: list of nx3 coordinates of the simplified mesh
% - faces2: triangular faces of vertices2 resulting froma a Delaunay
% triangulation
%
% (c) Arnadi Murtiyoso (INSA Strasbourg - ICube-TRIO UMR 7357)
NewClusterStruct=inputCluster;
%% MERGE THE MESHES
nbRegions=numel(fieldnames(inputCluster));
objectList=strings(nbRegions,1);
% create a list containing Object names (useful later)
for j=1:nbRegions
objectList(j,1)=strcat('Object',num2str(j));
[nbVertices,~]=size(inputCluster.(objectList(j,1)).Mesh.Vertices);
% create a new field with index to the original object
v = zeros(nbVertices, 1);
v(:) = j;
inputCluster.(objectList(j,1)).Mesh.ObjectIndex = v;
end
% initialise variables
MergedVertices=inputCluster.Object1.Mesh.Vertices;
MergedFaces=inputCluster.Object1.Mesh.Faces;
MergedIndex=inputCluster.Object1.Mesh.ObjectIndex;
i=1;
% start the loop
while i<nbRegions
objectName1=objectList(i);
objectName2=objectList(i+1);
% stack the vertex lists
MergedVertices = vertcat(MergedVertices, inputCluster.(objectName2).Mesh.Vertices);
% update the face list
faceUpdate = inputCluster.(objectName2).Mesh.Faces+max(max(MergedFaces));
MergedFaces = vertcat(MergedFaces,faceUpdate);
% update the object index list
MergedIndex = vertcat(MergedIndex, inputCluster.(objectName2).Mesh.ObjectIndex);
i=i+1;
end
%% COMPUTE "SNAPPING"
%initialise the variables
ClustersCopy=inputCluster; % "dynamic" copy of the original Cluster struct
VertexList=MergedVertices; % "dynamic" list of vertices
MergedIndex2=MergedIndex; %"dynamic" list of indices
% start the loop, continue as long as the list of vertices is not empty
while ~isempty(VertexList)
% take the first seed
vertex1 = VertexList(1,:);
% change VertexList to Matlab point cloud struct, necessary for the
% next function
VertexPointCloud=pointCloud(VertexList);
% find the seed's 10 nearest neighbours (this value can be changed, but
% I doubt it would be necessary)
[indices,dists] = findNearestNeighbors(VertexPointCloud,vertex1,10);
% for now, delete the seed from the list of nearest neighbours
indices(1,:)=[];
dists(1,:)=[];
% number of neighbours inside the tolerance radius
nb = length(indices);
% initialise a dynamic list of neighbours and their pointer towards the
% original Object
NeighborList = []; % empty list of neighbours inside the radius
OriginList = []; % empty list of neighbours' origins
% start a loop to check each neighbour
for i=1:nb
% if the neighbour point's distance is less than the tolerance,
% consider it as the same point i.e. "snapped"
if dists(i,1)<tol
a=MergedIndex2;
NeighborList(i,1)=indices(i,1);
OriginList(i,1)= a(indices(i,1));
else
continue
end
end
% if the seed has neighbours...
if ~isempty(NeighborList)
% return the seed to both lists, necessary to compute the
% coordinates of the new "snapped" point
NeighborList=[NeighborList;1];
OriginList=[OriginList;MergedIndex2(1)];
% number of items in the list
szList=length(NeighborList);
% how many unique Objects?
UniqueOrigins=unique(OriginList);
szUnique=length(UniqueOrigins);
% create a recapitulative table
tableToDelete=zeros(szList,2);
for x=1:szList
originNm=strcat('Object',string(OriginList(x)));
% retrieve the vector to be deleted from reference list of
% vertices
VectToDelete = VertexList(NeighborList(x),:);
% find its index in the original Object vertices
[~, boolDel]=ismember(ClustersCopy.(originNm).Mesh.Vertices,...
VectToDelete);
indexDel=find(boolDel,1);
% update the recapitulative table
tableToDelete(x,1) = OriginList(x);
tableToDelete(x,2) = indexDel;
end
% delete the points who
8000
se indices are listed in NeighborList in
% respective original Objects
for x=1:szUnique
originNm=strcat('Object',string(UniqueOrigins(x)));
% find the indices of This Object to be deleted from the
% recapitulative table
ind1 = tableToDelete(:,1) == UniqueOrigins(x);
thisOriginDel = tableToDelete(ind1,2);
% delete the said point in the dynamic/copied cluster
ClustersCopy.(originNm).Mesh.Vertices(thisOriginDel,:)=[];
% keep only unique points and delete (possible) redundants
ClustersCopy.(originNm).Mesh.Vertices=...
unique(ClustersCopy.(originNm).Mesh.Vertices,'rows');
end
% initialise XYZ vectors of all the neighbours and the seed
x_vect=zeros(szList,1);
y_vect=zeros(szList,1);
z_vect=zeros(szList,1);
% retrieve the neighbours and seed's XYZ
for j=1:szList
x_vect(j)=VertexList(NeighborList(j),1);
y_vect(j)=VertexList(NeighborList(j),2);
z_vect(j)=VertexList(NeighborList(j),3);
end
% compute the median to generate a new point
% why median? more robust than mean/average
x_new = median(x_vect);
y_new = median(y_vect);
z_new = median(z_vect);
new = [x_new y_new z_new];
% add planarity constraint
listPlanes=zeros(szUnique,9);
listPlanesParameters=zeros(szUnique,4);
% if the neighbours belong to more than 2 planes...
if szUnique>2
for k=1:szUnique
originNm = strcat('Object',string(UniqueOrigins(k)));
listPlanesParameters(k,:)=...
inputCluster.(originNm).PlaneMatlab.Parameters;
end
% use Geom3D's three planes intersection function
% plane1 = listPlanes(1,:);
% plane2 = listPlanes(2,:);
% plane3 = listPlanes(3,:);
% new2 = intersectThreePlanes(plane1, plane2, plane3);
% use my least squares multiple planes function
[new2,~]=intersectMultPlanes(listPlanesParameters);
% if the neighbours belong to two planes...
elseif szUnique==2
for k=1:2
originNm = strcat('Object',string(UniqueOrigins(k)));
listPlanes(k,:)=inputCluster.(originNm).PlaneGeom3D;
end
plane1 = listPlanes(1,:);
plane2 = listPlanes(2,:);
% intersect the two planes to get a 3D line
line = intersectPlanes(plane1, plane2);
% project the computed median point into this 3D line
new2 = projPointOnLine3d(new, line);
else
new2=new;
end
% stact the new point to the un-merged objects
if szUnique>1
for x=1:szUnique
originNm = strcat('Object',string(UniqueOrigins(x)));
ClustersCopy.(originNm).Mesh.Vertices= ...
[ClustersCopy.(originNm).Mesh.Vertices;new2];
end
else
originNm = strcat('Object',string(UniqueOrigins(1)));
ClustersCopy.(originNm).Mesh.Vertices=...
[ClustersCopy.(originNm).Mesh.Vertices;new2];
end
% delete the snapped points from the two lists
indNeighbors=NeighborList';
VertexList(indNeighbors,:)=[];
MergedIndex2(indNeighbors,:)=[];
else
% delete the seed from the two lists if they don't have neighbours
VertexList(1,:)=[];
MergedIndex2(1,:)=[];
end
end
%% VISUALISE THE 'SNAPPED' RESULT AND STORE IN NEW STRUCT
figure('name','Detected Roof Vertices - Delaunay Mesh')
for x=1:nbRegions
originNm=strcat('Object',string(x));
% Delaunay triangulation
TR = delaunayTriangulation...
(double(ClustersCopy.(originNm).Mesh.Vertices(:,1)),...
double(ClustersCopy.(originNm).Mesh.Vertices(:,2)),...
double(ClustersCopy.(originNm).Mesh.Vertices(:,3)));
% retrieve the free boundary mesh
[F,xf] = freeBoundary(TR);
[V2, F2] = ensureManifoldMesh(xf, F);
% store the new vertices and faces in a new struct
NewClusterStruct.(originNm).Mesh.Vertices=V2;
NewClusterStruct.(originNm).Mesh.Faces=F2;
% plot the mesh
pcshow(ClustersCopy.(originNm).PtCloud)
hold on
drawMesh(V2,F2,'b');
hold on
end
end