1
+ function [V ,S ] = alphavol(X ,R ,fig )
2
+ % ALPHAVOL Alpha shape of 2D or 3D point set.
3
+ % V = ALPHAVOL(X,R) gives the area or volume V of the basic alpha shape
4
+ % for a 2D or 3D point set. X is a coordinate matrix of size Nx2 or Nx3.
5
+ %
6
+ % R is the probe radius with default value R = Inf. In the default case
7
+ % the basic alpha shape (or alpha hull) is the convex hull.
8
+ %
9
+ % [V,S] = ALPHAVOL(X,R) outputs a structure S with fields:
10
+ % S.tri - Triangulation of the alpha shape (Mx3 or Mx4)
11
+ % S.vol - Area or volume of simplices in triangulation (Mx1)
12
+ % S.rcc - Circumradius of simplices in triangulation (Mx1)
13
+ % S.bnd - Boundary facets (Px2 or Px3)
14
+ %
15
+ % ALPHAVOL(X,R,1) plots the alpha shape.
16
+ %
17
+ % % 2D Example - C shape
18
+ % t = linspace(0.6,5.7,500)';
19
+ % X = 2*[cos(t),sin(t)] + rand(500,2);
20
+ % subplot(221), alphavol(X,inf,1);
21
+ % subplot(222), alphavol(X,1,1);
22
+ % subplot(223), alphavol(X,0.5,1);
23
+ % subplot(224), alphavol(X,0.2,1);
24
+ %
25
+ % % 3D Example - Sphere
26
+ % [x,y,z] = sphere;
27
+ % [V,S] = alphavol([x(:),y(:),z(:)]);
28
+ % trisurf(S.bnd,x,y,z,'FaceColor','blue','FaceAlpha',1)
29
+ % axis equal
30
+ %
31
+ % % 3D Example - Ring
32
+ % [x,y,z] = sphere;
33
+ % ii = abs(z) < 0.4;
34
+ % X = [x(ii),y(ii),z(ii)];
35
+ % X = [X; 0.8*X];
36
+ % subplot(211), alphavol(X,inf,1);
37
+ % subplot(212), alphavol(X,0.5,1);
38
+ %
39
+ % See also DELAUNAY, TRIREP, TRISURF
40
+
41
+ % Author: Jonas Lundgren <[email protected] > 2010
42
+
43
+ % 2010-09-27 First version of ALPHAVOL.
44
+ % 2010-10-05 DelaunayTri replaced by DELAUNAYN. 3D plots added.
45
+ % 2012-03-08 More output added. DELAUNAYN replaced by DELAUNAY.
46
+
47
+ if nargin < 2 || isempty(R ), R = inf ; end
48
+ if nargin < 3 , fig = 0 ; end
49
+
50
+ % Check coordinates
51
+ dim = size(X ,2 );
52
+ if dim < 2 || dim > 3
53
+ error(' alphavol:dimension' ,' X must have 2 or 3 columns.' )
54
+ end
55
+
56
+ % Check probe radius
57
+ if ~isscalar(R ) || ~isreal(R ) || isnan(R )
58
+ error(' alphavol:radius' ,' R must be a real number.' )
59
+ end
60
+
61
+ % Unique points
62
+ [X ,imap ] = unique(X ,' rows' );
63
+
64
+ % Delaunay triangulation
65
+ T = delaunay(X );
66
+
67
+ % Remove zero volume tetrahedra since
68
+ % these can be of arbitrary large circumradius
69
+ if dim == 3
70
+ n = size(T ,1 );
71
+ vol = volumes(T ,X );
72
+ epsvol = 1e- 12 * sum(vol )/n ;
73
+ T = T(vol > epsvol ,: );
74
+ holes = size(T ,1 ) < n ;
75
+ end
76
+
77
+ % Limit circumradius of simplices
78
+ [~ ,rcc ] = circumcenters(TriRep(T ,X ));
79
+ T = T(rcc < R ,: );
80
+ rcc = rcc(rcc < R );
81
+
82
+ % Volume/Area of alpha shape
83
+ vol = volumes(T ,X );
84
+ V = sum(vol );
85
+
86
+ % Return?
87
+ if nargout < 2 && ~fig
88
+ return
89
+ end
90
+
91
+ % Turn off TriRep warning
92
+ warning(' off' ,' MATLAB:TriRep:PtsNotInTriWarnId' )
93
+
94
+ % Alpha shape boundary
95
+ if ~isempty(T )
96
+ % Facets referenced by only one simplex
97
+ B = freeBoundary(TriRep(T ,X ));
98
+ if dim == 3 && holes
99
+ % The removal of zero volume tetrahedra causes false boundary
100
+ % faces in the interior of the volume. Take care of these.
101
+ B = trueboundary(B ,X );
102
+ end
103
+ else
104
+ B = zeros(0 ,dim );
105
+ end
106
+
107
+ % Plot alpha shape
108
+ if fig
109
+ figure(1 )
110
+ if dim == 2
111
+ % Plot boundary edges and point set
112
+ x = X(: ,1 );
113
+ y = X(: ,2 );
114
+ plot(x(B )' ,y(B )' ,' r' ,' linewidth' ,2 ), hold on
115
+ plot(x ,y ,' k.' ), hold off
116
+ str = ' Area' ;
117
+ elseif ~isempty(B )
118
+ % Plot boundary faces
119
+ trisurf(TriRep(B ,X ),' FaceColor' ,' red' ,' FaceAlpha' ,1 / 3 );
120
+ str = ' Volume' ;
121
+ else
122
+ cla
123
+ str = ' Volume' ;
124
+ end
125
+ axis equal
126
+ str = sprintf(' Radius = %g , %s = %g ' , R , str , V );
127
+ title(str ,' fontsize' ,12 )
128
+ end
129
+
130
+ % Turn on TriRep warning
131
+ warning(' on' ,' MATLAB:TriRep:PtsNotInTriWarnId' )
132
+
133
+ % Return structure
134
+ if nargout == 2
135
+ S = struct(' tri' ,imap(T ),' vol' ,vol ,' rcc' ,rcc ,' bnd' ,imap(B ));
136
+ end
137
+
138
+
139
+ % --------------------------------------------------------------------------
140
+ function vol = volumes(T ,X )
141
+ % VOLUMES Volumes/areas of tetrahedra/triangles
142
+
143
+ % Empty case
144
+ if isempty(T )
145
+ vol = zeros(0 ,1 );
146
+ return
147
+ end
148
+
149
+ % Local coordinates
150
+ A = X(T(: ,1 ),: );
151
+ B = X(T(: ,2 ),: ) - A ;
152
+ C = X(T(: ,3 ),: ) - A ;
153
+
154
+ if size(X ,2 ) == 3
155
+ % 3D Volume
156
+ D = X(T(: ,4 ),: ) - A ;
157
+ BxC = cross(B ,C ,2 );
158
+ vol = dot(BxC ,D ,2 );
159
+ vol = abs(vol )/6 ;
160
+ else
161
+ % 2D Area
162
+ vol = B(: ,1 ).*C(: ,2 ) - B(: ,2 ).*C(: ,1 );
163
+ vol = abs(vol )/2 ;
164
+ end
165
+
166
+
167
+ % --------------------------------------------------------------------------
168
+ function B = trueboundary(B ,X )
169
+ % TRUEBOUNDARY True boundary faces
170
+ % Remove false boundary caused by the removal of zero volume
171
+ % tetrahedra. The input B is the output of TriRep/freeBoundary.
172
+
173
+ % Surface triangulation
174
+ facerep = TriRep(B ,X );
175
+
176
+ % Find edges attached to two coplanar faces
177
+ E0 = edges(facerep );
178
+ E1 = featureEdges(facerep , 1e-6 );
179
+ E2 = setdiff(E0 ,E1 ,' rows' );
180
+
181
+ % Nothing found
182
+ if isempty(E2 )
183
+ return
184
+ end
185
+
186
+ % Get face pairs attached to these edges
187
+ % The edges connects faces into planar patches
188
+ facelist = edgeAttachments(facerep ,E2 );
189
+ pairs = cell2mat(facelist );
190
+
191
+ % Compute planar patches (connected regions of faces)
192
+ n = size(B ,1 );
193
+ C = sparse(pairs(: ,1 ),pairs(: ,2 ),1 ,n ,n );
194
+ C = C + C ' + speye(n );
195
+ [~ ,p ,r ] = dmperm(C );
196
+
197
+ % Count planar patches
198
+ iface = diff(r );
199
+ num = numel(iface );
200
+
201
+ % List faces and vertices in patches
202
+ facelist = cell(num ,1 );
203
+ vertlist = cell(num ,1 );
204
+ for k = 1 : num
205
+
206
+ % Neglect singel face patches, they are true boundary
207
+ if iface(k ) > 1
208
+
209
+ % List of faces in patch k
210
+ facelist{k } = p(r(k ): r(k + 1 )-1 );
211
+
212
+ % List of unique vertices in patch k
213
+ vk = B(facelist{k },: );
214
+ vk = sort(vk(: ))' ;
215
+ ik = [true ,diff(vk )>0];
216
+ vertlist{k } = vk(ik );
217
+
218
+ end
219
+ end
220
+
221
+ % Sort patches by number of vertices
222
+ ivert = cellfun(@numel ,vertlist );
223
+ [ivert ,isort ] = sort(ivert );
224
+ facelist = facelist(isort );
225
+ vertlist = vertlist(isort );
226
+
227
+ % Group patches by number of vertices
228
+ p = [0 ;find(diff(ivert ));num ] + 1 ;
229
+ ipatch = diff(p );
230
+
231
+ % Initiate true boundary list
232
+ ibound = true(n ,1 );
233
+
234
+ % Loop over groups
235
+ for k = 1 : numel(ipatch )
236
+
237
+ % Treat groups with at least two patches and four vertices
238
+ if ipatch(k ) > 1 && ivert(p(k )) > 3
239
+
240
+ % Find double patches (identical vertex rows)
241
+ V = cell2mat(vertlist(p(k ): p(k + 1 )-1 ));
242
+ [V ,isort ] = sortrows(V );
243
+ id = ~any(diff(V ),2 );
244
+ id = [id ;0] | [0 ;id ];
245
+ id(isort ) = id ;
246
+
247
+ % Deactivate faces in boundary list
248
+ for j = find(id ' )
249
+ ibound(facelist{j - 1 + p(k )}) = 0 ;
250
+ end
251
+
252
+ end
253
+ end
254
+
255
+ % Remove false faces
256
+ B = B(ibound ,: );
0 commit comments