forked from fangq/iso2mesh
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Content.txt
2120 lines (1979 loc) · 93.1 KB
/
Content.txt
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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
=== # Streamlined mesh generation - shortcuts ===
==== function [node,elem,face]=v2m(img,isovalues,opt,maxvol,method) ====
[node,elem,face]=v2m(img,isovalues,opt,maxvol,method)
volumetric mesh generation from binary or gray-scale volumetric images
shortcut for vol2mesh
inputs and outputs are similar to those defined in vol2mesh
==== function [no,el,regions,holes]=v2s(img,isovalues,opt,method) ====
[no,el,regions,holes]=v2s(img,isovalues,opt,method)
surface mesh generation from binary or gray-scale volumetric images
shortcut for vol2surf
inputs and outputs are similar to those defined in vol2surf; In v2s,
method can be set to 'cgalmesh' in addition to those allowed by vol2surf.
==== function [node,elem,face]=s2m(v,f,keepratio,maxvol,method,regions,holes) ====
[node,elem,face]=s2m(v,f,keepratio,maxvol,method)
[node,elem,face]=s2m(v,f,keepratio,maxvol,'tetgen',regions,holes)
volumetric mesh generation from a closed surface, shortcut for surf2mesh
inputs and outputs are similar to those defined in surf2mesh
if method='cgalpoly', s2m will call cgals2m and keepratio should be a
structure (as the 'opt' input in cgals2m)
input default values:
method: if ignored, iso2mesh uses surf2mesh ('tetgen') to do the
tetrahedral mesh generation
regions,holes: if ignored, iso2mesh assumes both are empty
==== function img=s2v(node,face,div) ====
img=s2v(node,face,div)
shortcut for surf2vol, coverting a surface to a volumetric image
input:
node: node list of the triangular surface, 3 columns for x/y/z
face: triangle node indices, each row is a triangle
div: division number along the shortest edge of the mesh (resolution)
if not given, div=50
output:
img: a volumetric binary image at position of ndgrid(xi,yi,zi)
==== function varargout=m2v(varargin) ====
vol=m2v(node,face,Nxyz)
or
vol=m2v(node,face,xi,yi,zi)
shortcut for mesh2vol, rasterizing a teterahedral mesh to a volume using graphics
input/output: please see details in the help for mesh2vol
==== function newnode=sms(node,face,iter,alpha,method) ====
newnode=sms(node,face,iter,useralpha,method)
simplified version of surface mesh smoothing
input:
node: node coordinates of a surface mesh
face: face element list of the surface mesh
iter: smoothing iteration number
alpha: scaler, smoothing parameter, v(k+1)=alpha*v(k)+(1-alpha)*mean(neighbors)
method: same as in smoothsurf, default is 'laplacianhc'
output:
newnode: output, the smoothed node coordinates
=== # Streamlined mesh generation ===
==== function [node,elem,face,regions]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues) ====
[node,elem,face,regions]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues)
convert a binary (or multi-valued) volume to tetrahedral mesh
input:
img: a volumetric binary image
ix,iy,iz: subvolume selection indices in x,y,z directions
opt: as defined in vol2surf.m
maxvol: target maximum tetrahedral elem volume
when method='cgalmesh', maxvol can specify the target
for each label (subregion index) by the following syntax
'label1=size1:label2=size2:...'
dofix: 1: perform mesh validation&repair, 0: skip repairing
method: 'cgalsurf' or omit: use CGAL surface mesher
'simplify': use binsurface and then simplify
'cgalmesh': use CGAL 3.5 3D mesher for direct mesh generation [new]
generally speaking, 'cgalmesh' is the most robust path
if you want to product meshes from binary or multi-region
volumes, however, its limitations include 1) only accept
uint8 volume, and 2) can not extract meshes from gray-scale
volumes. If ones goal is to process a gray-scale volume,
he/she should use the 'cgalsurf' option. 'simplify' approach
is not recommended unless other options has failed.
isovalues: a list of isovalues where the levelset is defined
output:
node: output, node coordinates of the tetrahedral mesh
elem: output, element list of the tetrahedral mesh, the last
column is the region ID
face: output, mesh surface element list of the tetrahedral mesh
the last column denotes the boundary ID
region: optional output. if opt.autoregion is set to 1, region
saves the interior points for each closed surface component
==== function [no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues) ====
[no,el,regions,holes]=vol2surf(img,ix,iy,iz,opt,dofix,method,isovalues)
converting a 3D volumetric image to surfaces
input:
img: a volumetric binary image; if img is empty, vol2surf will
return user defined surfaces via opt.surf if it exists
ix,iy,iz: subvolume selection indices in x,y,z directions
opt: function parameters
if method is 'cgalsurf' or 'cgalpoly':
opt=a float number>1: max radius of the Delaunay sphere(element size)
opt.radbound: same as above, max radius of the Delaunay sphere
opt.distbound: maximum deviation from the specified isosurfaces
opt(1,2,...).radbound: same as above, for each levelset
if method is 'simplify':
opt=a float number<1: compression rate for surf. simplification
opt.keepratio=a float less than 1: same as above, same for all surf.
opt(1,2,..).keepratio: setting compression rate for each levelset
opt(1,2,..).maxsurf: 1 - only use the largest disjointed surface
0 - use all surfaces for that levelset
opt(1,2,..).side: - 'upper': threshold at upper interface
'lower': threshold at lower interface
opt(1,2,..).maxnode: - the maximum number of surface node per levelset
opt(1,2,..).holes: user specified holes interior pt list
opt(1,2,..).regions: user specified regions interior pt list
opt(1,2,..).surf.{node,elem}: add additional surfaces
opt(1,2,..).{A,B}: linear transformation for each surface
opt.autoregion: if set to 1, vol2surf will try to determine
the interior points for each closed surface automatically
dofix: 1: perform mesh validation&repair, 0: skip repairing
method: - if method is 'simplify', iso2mesh will first call
binsurface to generate a voxel-based surface mesh and then
use meshresample/meshcheckrepair to create a coarser mesh;
- if method is 'cgalsurf', iso2mesh will call the surface
extraction program from CGAL to make surface mesh
- if method is not specified, 'cgalsurf' is assumed by default
isovalues: a list of isovalues where the levelset is defined
output:
no: list of nodes on the resulting suface mesh, 3 columns for x,y,z
el: list of trianglular elements on the surface, [n1,n2,n3,region_id]
regions: list of interior points for all sub-region, [x,y,z]
holes: list of interior points for all holes, [x,y,z]
==== function [node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox) ====
[node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox)
create quality volumetric mesh from isosurface patches
input parameters:
v: input, isosurface node list, dimension (nn,3)
if v has 4 columns, the last column specifies mesh density near each node
f: input, isosurface face element list, dimension (be,3)
p0: input, coordinates of one corner of the bounding box, p0=[x0 y0 z0]
p1: input, coordinates of the other corner of the bounding box, p1=[x1 y1 z1]
keepratio: input, percentage of elements being kept after the simplification
maxvol: input, maximum tetrahedra element volume
regions: list of regions, specifying by an internal point for each region
holes: list of holes, similar to regions
forcebox: 1: add bounding box, 0: automatic
outputs:
node: output, node coordinates of the tetrahedral mesh
elem: output, element list of the tetrahedral mesh
face: output, mesh surface element list of the tetrahedral mesh
the last column denotes the boundary ID
==== function img=surf2vol(node,face,xi,yi,zi) ====
img=surf2vol(node,face,xi,yi,zi)
convert a triangular surface to a shell of voxels in a 3D image
input:
node: node list of the triangular surface, 3 columns for x/y/z
face: triangle node indices, each row is a triangle
xi,yi,zi: x/y/z grid for the resulting volume
output:
img: a volumetric binary image at position of ndgrid(xi,yi,zi)
==== function [mask weight]=mesh2vol(node,elem,xi,yi,zi) ====
[mask weight]=mesh2vol(node,face,Nxyz)
[mask weight]=mesh2vol(node,face,[Nx,Ny,Nz])
[mask weight]=mesh2vol(node,face,xi,yi,zi,hf)
or
newval=mesh2vol(node_val,face,...)
fast rasterization of a 3D mesh to a volume with tetrahedron index labels
date for initial version: Feb 10,2014
input:
node: node coordinates, dimension N by 2 or N by 3 array
nodeval: a 4-column array, the first 3 columns are the node coordinates,
the last column denotes the values associated with each node
face: a triangle surface, N by 3 or N by 4 array
Nx,Ny,Nxy: output image in x/y dimensions, or both
xi,yi: linear vectors for the output pixel center positions in x/y
hf: the handle of a pre-created figure window for faster rendering
output:
mask: a 3D image, the value of each pixel is the index of the
enclosing triangle, if the pixel is outside of the mesh, NaN
weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions
for the mask
newval: when the node has 4 columns, the last column represents the
values (color) at each node, the output newval is the rasterized
mesh value map over the specified grid.
note: This function only works for matlab
example:
[no,el]=meshgrid6(0:5,0:5,0:3);
mask=mesh2vol(no,el(:,1:4),0:0.1:5,0:0.1:5,0:0.1:4);
imagesc(mask(:,:,3))
=== # iso2mesh main function backend ===
==== function [node,elem]=binsurface(img,nface) ====
[node,elem]=binsurface(img,nface)
fast isosurface extraction from 3D binary images
input:
img: a 3D binary image
nface: nface=3 or ignored - for triangular faces,
nface=4 - square faces
nface=0 - return a boundary mask image via node
output:
elem: integer array with dimensions of NE x nface, each row represents
a surface mesh face element
node: node coordinates, 3 columns for x, y and z respectively
the outputs of this subroutine can be easily plotted using
patch('Vertices',node,'faces',elem,'FaceVertexCData',node(:,3),
'FaceColor','interp');
if the surface mesh has triangular faces, one can plot it with
trisurf(elem,node(:,1),node(:,2),node(:,3))
==== function [node,elem,face]=cgalv2m(vol,opt,maxvol) ====
[node,elem,face]=cgalv2m(vol,opt,maxvol)
wrapper for CGAL 3D mesher (CGAL 3.5 or up)
convert a binary (or multi-valued) volume to tetrahedral mesh
http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html
input:
vol: a volumetric binary image
ix,iy,iz: subvolume selection indices in x,y,z directions
opt: parameters for CGAL mesher, if opt is a structure, then
opt.radbound: defines the maximum surface element size
opt.angbound: defines the miminum angle of a surface triangle
opt.distbound: defines the maximum distance between the
center of the surface bounding circle and center of the
element bounding sphere
opt.reratio: maximum radius-edge ratio
if opt is a scalar, it only specifies radbound.
maxvol: target maximum tetrahedral elem volume
output:
node: output, node coordinates of the tetrahedral mesh
elem: output, element list of the tetrahedral mesh, the last
column is the region id
face: output, mesh surface element list of the tetrahedral mesh
the last column denotes the boundary ID
note: each triangle will appear twice in the face list with each
one attaches to each side of the interface. one can remove
the redundant triangles by unique(face(:,1:3),'rows')
==== function [node,elem,face]=cgals2m(v,f,opt,maxvol,varargin) ====
[node,elem,face]=cgals2m(v,f,opt,maxvol)
wrapper for CGAL 3D mesher (CGAL 3.5 and newer)
convert a triangular surface to tetrahedral mesh
http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html
input:
v: the node coordinate list of a surface mesh (nn x 3)
f: the face element list of a surface mesh (be x 3)
opt: parameters for CGAL mesher, if opt is a structure, then
opt.radbound: defines the maximum surface element size
opt.angbound: defines the miminum angle of a surface triangle
opt.distbound: defines the maximum distance between the
center of the surface bounding circle and center of the
element bounding sphere
opt.reratio: maximum radius-edge ratio
if opt is a scalar, it only specifies radbound.
maxvol: target maximum tetrahedral elem volume
output:
node: output, node coordinates of the tetrahedral mesh
elem: output, element list of the tetrahedral mesh, the last
column is the region id
face: output, mesh surface element list of the tetrahedral mesh
the last column denotes the boundary ID
==== function [node,elem]=vol2restrictedtri(vol,thres,cent,brad,ang,radbound,distbound,maxnode) ====
[node,elem]=vol2restrictedtri(vol,thres,cent,brad,ang,radbound,distbound,maxnode)
surface mesh extraction using CGAL mesher
input:
vol: a 3D volumetric image
thres: a scalar as the threshold of of the extraction
cent: a 3d position (x,y,z) which locates inside the resulting
mesh, this is automatically computed from vol2surf
brad: maximum bounding sphere squared of the resulting mesh
ang: minimum angular constrains of the resulting tranglar elements
(in degrees)
radbound: maximum triangle delaunay circle radius
distbound: maximum delaunay sphere distances
maxnode: maximum number of surface nodes (even radbound is not reached)
output:
node: the list of 3d nodes in the resulting surface (x,y,z)
elem: the element list of the resulting mesh (3 columns of integers)
==== function img=surf2volz(node,face,xi,yi,zi) ====
img=surf2volz(node,face,xi,yi,zi)
convert a triangular surface to a shell of voxels in a 3D image
along the z-axis
input:
node: node list of the triangular surface, 3 columns for x/y/z
face: triangle node indices, each row is a triangle
xi,yi,zi: x/y/z grid for the resulting volume
output:
img: a volumetric binary image at position of ndgrid(xi,yi,zi)
==== function [mask, weight]=mesh2mask(node,face,xi,yi,hf) ====
[mask weight]=mesh2mask(node,face,Nxy)
or
[mask weight]=mesh2mask(node,face,[Nx,Ny])
or
[mask weight]=mesh2mask(node,face,xi,yi,hf)
fast rasterization of a 2D mesh to an image with triangle index labels
date for initial version: July 18,2013
input:
node: node coordinates, dimension N by 2 or N by 3 array
face: a triangle surface, N by 3 or N by 4 array
Nx,Ny,Nxy: output image in x/y dimensions, or both
xi,yi: linear vectors for the output pixel center positions in x/y
hf: (optional) the handle of a pre-created figure window, for faster
rendering
output:
mask: a 2D image, the value of each pixel is the index of the
enclosing triangle, if the pixel is outside of the mesh, NaN
weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions for
the mask
note: This function only works in MATLAB when the DISPLAY is not
disabled. The maximum size of the mask output is limited by the
screen size.
example:
[no,fc]=meshgrid6(0:5,0:5);
[mask weight]=mesh2mask(no,fc,-1:0.1:5,0:0.1:5);
imagesc(mask);
=== # iso2mesh primitive meshing functions ===
==== function [node,face,elem]=meshabox(p0,p1,opt,nodesize) ====
[node,face,elem]=meshabox(p0,p1,opt,maxvol)
create the surface and tetrahedral mesh of a box geometry
input:
p0: coordinates (x,y,z) for one end of the box diagnoal
p1: coordinates (x,y,z) for the other end of the box diagnoal
opt: maximum volume of the tetrahedral elements
nodesize: 1 or a 8x1 array, size of the element near each vertex
output:
node: node coordinates, 3 columns for x, y and z respectively
face: integer array with dimensions of NB x 3, each row represents
a surface mesh face element
elem: integer array with dimensions of NE x 4, each row represents
a tetrahedron
example:
[node,face,elem]=meshabox([2 3 2],[6 12 15],0.1,1);
plotmesh(node,elem,'x>4');
==== function [node,face,elem]=meshasphere(c0,r,tsize,maxvol) ====
[node,face,elem]=meshasphere(c0,r,tsize,maxvol)
create the surface and tetrahedral mesh of a sphere
input:
c0: center coordinates (x0,y0,z0) of the sphere
r: radius of the sphere
tsize: maximum surface triangle size on the sphere
maxvol: maximu volume of the tetrahedral elements
output:
node: node coordinates, 3 columns for x, y and z respectively
face: integer array with dimensions of NB x 3, each row represents
a surface mesh face element
elem: integer array with dimensions of NE x 4, each row represents
a tetrahedron
==== function [node,face,elem]=meshanellip(c0,rr,tsize,maxvol) ====
[node,face,elem]=meshanellip(c0,rr,opt)
create the surface and tetrahedral mesh of an ellipsoid
input:
c0: center coordinates (x0,y0,z0) of the ellipsoid
rr: radii of an ellipsoid,
if rr is a scalar, this is a sphere with radius rr
if rr is a 1x3 or 3x1 vector, it specifies the ellipsoid radii [a,b,c]
if rr is a 1x5 or 5x1 vector, it specifies [a,b,c,theta,phi]
where theta and phi are the rotation angles along z and x
axes, respectively. Rotation is applied before translation.
tsize: maximum surface triangle size on the sphere
maxvol: maximu volume of the tetrahedral elements
output:
node: node coordinates, 3 columns for x, y and z respectively
face: integer array with dimensions of NB x 3, each row represents
a surface mesh face element
elem: integer array with dimensions of NE x 4, each row represents
a tetrahedron; if ignored, only produces the surface
example:
[node,face,elem]=meshanellip([10,10,-10],[30,20,10,pi/4,pi/4],0.5,0.4);
plotmesh(node,elem,'x>10');axis equal;
==== function [node,face,elem]=meshunitsphere(tsize,maxvol) ====
[node,face,elem]=meshunitsphere(tsize,maxvol)
create the surface and/or volumetric mesh of a unit sphere
centered at [0 0 0] and radius 1
input:
tsize: maximum size of the surface triangles (from 0 to 1)
maxvol: maximum volume of the tetrahedron; if one wants to return
elem without specifying maxvol, maxvol=tsize^3
output:
node: node coordinates, 3 columns for x, y and z respectively
face: integer array with dimensions of NB x 3, each row represents
a surface mesh face element
elem: integer array with dimensions of NE x 4, each row represents
a tetrahedron. If ignored, this function only produces the surface
example:
[node,face]=meshunitsphere(0.05);
[node,face,elem]=meshunitsphere(0.05,0.01);
plotmesh(node,elem,'x>0'); axis equal;
==== function [node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv) ====
[node,face]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
or
[node,face,elem]=meshacylinder(c0,c1,r,tsize,maxvol,ndiv)
create the surface and (optionally) tetrahedral mesh of a 3D cylinder
input:
c0, c1: cylinder axis end points
r: radius of the cylinder
tsize: maximum surface triangle size on the sphere
maxvol: maximu volume of the tetrahedral elements
ndiv: approximate the cylinder surface into ndiv flat pieces, if
ignored, ndiv is set to 20
output:
node: node coordinates, 3 columns for x, y and z respectively
face: integer array with dimensions of NB x 3, each row represents
a surface mesh triangle
elem: (optional) integer array with dimensions of NE x 4, each row
represents a tetrahedron
==== function [node,elem]=meshgrid5(varargin) ====
[node,elem]=meshgrid5(v1,v2,v3,...)
mesh an ND rectangular lattice by splitting
each hypercube into 5 tetrahedra
inspired by John D'Errico
URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
input:
v1,v2,v3,... - numeric vectors defining the lattice in
each dimension.
Each vector must be of length >= 1
output:
node - factorial lattice created from (v1,v2,v3,...)
Each row of this array is one node in the lattice
elem - integer array defining simplexes as references to
rows of "node".
example:
[node,elem]=meshgrid5(0:5,0:6,0:4);
plotmesh(node,elem);
==== function [node,elem]=meshgrid6(varargin) ====
[node,elem]=meshgrid6(v1,v2,v3,...)
mesh an ND rectangular lattice by splitting
each hypercube into 6 tetrahedra
author: John D'Errico
URL: http://www.mathworks.com/matlabcentral/newsreader/view_thread/107191
input:
v1,v2,v3,... - numeric vectors defining the lattice in
each dimension.
Each vector must be of length >= 1
output:
node - factorial lattice created from (v1,v2,v3,...)
Each row of this array is one node in the lattice
elem - integer array defining simplexes as references to
rows of "node".
example:
[node,elem]=meshgrid6(0:5,0:6,0:4);
plotmesh(node,elem);
==== function [node,face,centroids]=latticegrid(varargin) ====
[node,face,centroids]=latticegrid(xrange,yrange,zrange,...)
generate a 3D lattice
input:
xrange, yrange, zrange ...: 1D vectors specifying the range of each
dimension of the lattice
output:
node: the vertices of the 3D lattice
face: the list of cell faces of the lattice, including both internal
and external facets. By default, face is in the form of a cell
array, with each row representing a face. One can use
cell2mat(face) to convert it to an array
centroids: the centroids of each lattice cell
example:
% generate a 3D lattice
[node,face,c0]=latticegrid([1 2 4],1:3,1:4);
plotmesh(node,face)
% mesh the 3D lattice based on the face info
[no,el]=surf2mesh(node,face,[],[],1,0.01,c0);
figure; plotmesh(no,el)
% mesh a 2-layer structure using a simple lattice
[node,face,c0]=latticegrid([0 10],[0 5],[0 3.5 4]);
c0(:,4)=[0.01;0.001];
[no,el]=surf2mesh(node,face,[],[],1,[],c0);
figure; plotmesh(no,el)
=== # Mesh decomposition and query ===
==== function facecell=finddisconnsurf(f) ====
facecell=finddisconnsurf(f)
subroutine to extract disconnected surfaces from a
cluster of surfaces
Date: 2008/03/06
input:
f: faces defined by node indices for all surface triangles
output:
facecell: separated disconnected surface node indices
==== function openedge=surfedge(f,varargin) ====
openedge=surfedge(f)
find the edge of an open surface or surface of a volume
input:
f: input, surface face element list, dimension (be,3)
output:
openedge: list of edges of the specified surface
==== function openface=volface(t) ====
openface=volface(t)
find the surface patches of a volume
input:
t: input, volumetric element list, dimension (ne,4)
output:
openface: list of faces of the specified volume
==== function loops=extractloops(edges) ====
loops=extractloops(edges)
extract individual loop or polyline segment from a collection of edges
input:
edges: two column matrix recording the starting/ending
points of all edge segments
output:
loops: output, a single vector separated by NaN, each segment
is a 3D polyline or loop consisted of node IDs
example:
edges=[1 2;2 3;1 4;3 4;7 3;1 9;5 6;6 7;10 9; 8 10;1 8;9 3;11 11;11 12];
loops=extractloops(edges)
==== function [conn,connnum,count]=meshconn(elem,nn) ====
[conn,connnum,count]=meshconn(elem,nn)
create node neighbor list from a mesh
input:
elem: element table of a mesh
nn : total node number of the mesh
output:
conn: output, a cell structure of length nn, conn{n}
contains a list of all neighboring node ID for node n
connnum: vector of length nn, denotes the neighbor number of each node
count: total neighbor numbers
==== function centroid=meshcentroid(v,f) ====
centroid=meshcentroid(v,f)
compute the centroids of a mesh defined by nodes and elements
(surface or tetrahedra) in R^n space
input:
v: surface node list, dimension (nn,3)
f: surface face element list, dimension (be,3)
output:
centroid: centroid positions, one row for each element
==== function nodevol=nodevolume(node,elem) ====
nodevol=nodevolume(node,elem)
calculate the Voronoi volume of each node in a simplex mesh
input:
node: node coordinates
elem: element table of a mesh
output:
nodevol: volume values for all nodes
==== function vol=elemvolume(node,elem,option) ====
vol=elemvolume(node,elem,option)
calculate the volume for a list of simplexes
input:
node: node coordinates
elem: element table of a mesh
option: if option='signed', the volume is the raw determinant,
else, the results will be the absolute values
output:
vol: volume values for all elements
==== function [conn,connnum,count]=neighborelem(elem,nn); ====
[conn,connnum,count]=neighborelem(elem,nn)
create node neighbor list from a mesh
input:
elem: element table of a mesh
nn : total node number of the mesh
output:
conn: output, a cell structure of length nn, conn{n}
contains a list of all neighboring elem ID for node n
connnum: vector of length nn, denotes the neighbor number of each node
count: total neighbor numbers
==== function facenb=faceneighbors(t,opt) ====
facenb=faceneighbors(t,opt)
to find 4 face-neighboring elements of a tetrahedron
input:
t: tetrahedron element list, 4 columns of integers
opt: if opt='surface', return boundary triangle list
(should be the same as the face output from v2m)
otherwise, return the element list for each element:
each row contains 4 numbers, representing the element
indices sharing triangular faces [1 2 3],[1 2 4],[1 3 4]
and [2 3 4] in order, where 1~4 is the node local index.
if the index is 0, indicating the face has no neighbor
(i.e. a boundary face)
output:
facenb: see opt
==== function edgenb=edgeneighbors(t,opt) ====
edgenb=edgeneighbors(t,opt)
to find neighboring triangular elements in a triangule surface
input:
t: a triangular surface element list, 3 columns of integers
opt: if opt='general', return the edge neighbors for a general
triangular surface: each edge can be shared by more than 2
triangles; if ignored, we assume all triangles are shared by no
more than 2 triangles.
output:
edgenb: if opt is not supplied, edgenb is a size(t,1) by 3 array with
each element being the triangle ID of the edge neighbor of that
triangle. For each row, the order of the neighbors is listed as those
sharing edges [1 2], [2 3] and [3 1] between the triangle nodes.
when opt='general', edgenb is a cell array with a length of size(t).
each member of the cell array is a list of edge neighbors (the order
is not defined).
==== function [f maxsize]=maxsurf(facecell,node) ====
[f maxsize]=maxsurf(facecell,node)
return the surface with the maximum element number or
total area from a cell arry of surfaces
input:
facecell: a cell array, each element is a face array
node: optional, node list, if given, the output is the
surface with the largest surface area.
output:
f: the surface data (node indices) for the surface with the
most elements (or largest area when node is given)
maxsize: if node is not given, maxisize is row number of f;
otherwise, maxsize is the total area of f
==== function mask=flatsegment(node,edge) ====
mask=flatsegment(node,edge)
decompose edge loops into flat segments alone arbitrary planes of the bounding box
this code is fragile: it can not handle curves with many co-linear
nodes near the corner point
input:
node: x,y,z coordinates of each node of the mesh
edge: input, a single vector separated by NaN, each segment
is a close-polygon consisted by node IDs
output:
mask: output, a cell, each element is a close-polygon
on x/y/z plane
==== function newedge=orderloopedge(edge) ====
[newedge]=orderloopedge(edge)
order the node list of a simple loop based on connection sequence
input:
edge: a loop consisted by a sequence of edges, each row
is an edge with two integers: start/end node index
output:
newedge: reordered edge node list
==== function [X,V,E,F]=mesheuler(face) ====
[X,V,E,F]=mesheuler(face)
Euler's charastistics of a mesh
input:
face: a closed surface mesh
output:
X: Euler's charastistics
V: number of vertices
E: number of edges
F: number of faces
==== function seg=bbxflatsegment(node,loop) ====
seg=bbxflatsegment(node,loop)
decompose edge loops into flat segments along the x/y/z
planes of the bounding box
input:
node: x,y,z coordinates of each node of the mesh
loop: input, a single vector separated by NaN, each segment
is a close-polygon consisted by node IDs
output:
seg: output, a single vector separated by NaN, each segment
is a close-polygon on x/y/z plane
==== function plane=surfplane(node,face) ====
plane=surfplane(node,face)
plane equation coefficients for each face in a surface
input:
node: a list of node coordinates (nn x 3)
face: a surface mesh triangle list (ne x 3)
output:
plane: a (ne x 4) array, in each row, it has [a b c d]
to denote the plane equation as "a*x+b*y+c*z+d=0"
==== function [pt,p0,v0,t,idx]=surfinterior(node,face) ====
[pt,p0,v0,t,idx]=surfinterior(node,face)
identify a point that is enclosed by the (closed) surface
input:
node: a list of node coordinates (nn x 3)
face: a surface mesh triangle list (ne x 3)
output:
pt: the interior point coordinates [x y z]
p0: ray origin used to determine the interior point
v0: the vector used to determine the interior point
t : ray-tracing intersection distances (with signs) from p0. the
intersection coordinates can be expressed as p0+t(i)*v0
idx: index to the face elements that intersect with the ray, order
match that of t
==== function elist=surfpart(f,loopedge) ====
elist=surfpart(f,loopedge)
partition a triangular surface using a closed loop defined by existing edges
input:
f: input, surface face element list, dimension (be,3)
loopedge: a 2-column array, specifying a closed loop in CCW order
output:
elist: list of triangles that is enclosed by the loop
==== function seeds=surfseeds(node,face) ====
seeds=surfseeds(node,face)
calculate a set of interior points with each enclosed by a closed
component of a surface
input:
node: a list of node coordinates (nn x 3)
face: a surface mesh triangle list (ne x 3)
output:
seeds: the interior points coordinates for each closed-surface
component
==== function quality=meshquality(node,elem,maxnode) ====
quality=meshquality(node,elem)
compute the Joe-Liu mesh quality measure of an N-D mesh (N<=3)
input:
node: node coordinates of the mesh (nn x 3)
elem: element table of an N-D mesh (ne x (N+1))
output:
quality: a vector of the same length as size(elem,1), with
each element being the Joe-Liu mesh quality metric (0-1) of
the corresponding element. A value close to 1 represents
higher mesh quality (1 means equilateral tetrahedron);
a value close to 0 means nearly degenerated element.
reference:
A. Liu, B. Joe, Relationship between tetrahedron shape measures,
BIT 34 (2) (1994) 268-287.
==== function edges=meshedge(elem,varargin) ====
edges=meshedge(elem,opt)
return all edges in a surface or volumetric mesh
input:
elem: element table of a mesh (support N-d space element)
opt: optional input, giving the additional options. If opt
is a struct, it can have the following field:
opt.nodeorder: if 1, assuming the elem node indices is in CCW
orientation; 0 use nchoosek() output to order edges
you can replace opt by a series of ('param', value) pairs.
output:
edge: edge list; each row is an edge, specified by the starting and
ending node indices, the total edge number is
size(elem,1) x nchoosek(size(elem,2),2). All edges are ordered
by looping through each element first.
==== function snorm=surfacenorm(node,face,varargin) ====
snorm=surfacenorm(node,face)
or
snorm=surfacenorm(node,face,'Normalize',0)
compute the normal vectors for a triangular surface
input:
node: a list of node coordinates (nn x 3)
face: a surface mesh triangle list (ne x 3)
opt: a list of optional parameters, currently surfacenorm supports:
'Normalize': [1|0] if set to 1, the normal vectors will be
unitary (default)
output:
snorm: output surface normal vector at each face
==== function [edges,idx,edgemap]=uniqedges(elem) ====
[edges,idx,edgemap]=uniqedges(elem)
return the unique edge list from a surface or tetrahedral mesh
input:
elem: a list of elements, each row is a list of nodes for an element.
elem can have 2, 3 or 4 columns
output:
edge: unique edges in the mesh, denoted by a pair of node indices
idx: index of the output in the raw edge list (returned by meshedge)
edgemap: index of the raw edges in the output list (for triangular mesh)
==== function [elist,nextfront]=advancefront(edges,loop,elen) ====
[elist,nextfront]=advancefront(edges,loop,elen)
advance an edge-front on an oriented surface to the next separated by
one-element width
input:
edges: edge list of an oriented surface mesh, must be in CCW order
loop: a 2-column array, specifying a closed loop in CCW order
elen: node number inside each element, if ignored, elen is set to 3
output:
elist: list of triangles that is enclosed between the two
edge-fronts
nextfront: a new edge loop list representing the next edge-front
=== # Mesh processing and reparing ===
==== function [node,elem]=meshcheckrepair(node,elem,opt,varargin) ====
[node,elem]=meshcheckrepair(node,elem,opt)
check and repair a surface mesh
input/output:
node: input/output, surface node list, dimension (nn,3)
elem: input/output, surface face element list, dimension (be,3)
opt: options, including
'dupnode': remove duplicated nodes
'dupelem' or 'duplicated': remove duplicated elements
'dup': both above
'isolated': remove isolated nodes
'open': abort when open surface is found
'deep': call external jmeshlib to remove non-manifold vertices
'meshfix': repair a closed surface by the meshfix utility (new)
it can remove self-intersecting elements and fill holes
'intersect': test a surface for self-intersecting elements
==== function newelem=meshreorient(node,elem) ====
newelem=meshreorient(node,elem)
reorder nodes in a surface or tetrahedral mesh to ensure all
elements are oriented consistently
input:
node: list of nodes
elem: list of elements (each row are indices of nodes of each element)
output:
newelem: the element list with consistent ordering
==== function elem=removedupelem(elem) ====
elem=removedupelem(elem)
remove doubly duplicated (folded) elements
input:
elem: list of elements (node indices)
output:
elem: element list after removing the duplicated elements
==== function [newnode,newelem]=removedupnodes(node,elem) ====
[newnode,newelem]=removedupnodes(node,elem)
removing the duplicated nodes from a mesh
input:
elem: integer array with dimensions of NE x 4, each row contains
the indices of all the nodes for each tetrahedron
node: node coordinates, 3 columns for x, y and z respectively
output:
newnode: nodes without duplicates
newelem: elements with only the unique nodes
==== function [no,el]=removeisolatednode(node,elem) ====
[no,el]=removeisolatednode(node,elem)
remove isolated nodes: nodes that are not included in any element
input:
node: list of node coordinates
elem: list of elements of the mesh
output:
no: node coordinates after removing the isolated nodes
el: element list of the resulting mesh
==== function fnew=removeisolatedsurf(v,f,maxdiameter) ====
fnew=removeisolatedsurf(v,f,maxdiameter)
remove disjointed surface fragment filtered by using mesh diameter
input:
v: list of nodes of the input surface
f: list of triangles of the input surface
maxdiameter: maximum bounding box size for surface removal
ouput:
fnew: new face list after removing the components smaller than
maxdiameter
==== function f=surfaceclean(f,v) ====
f=surfaceclean(f,v)
remove surface patches that are located inside
the bounding box faces
input:
v: surface node list, dimension (nn,3)
f: surface face element list, dimension (be,3)
output:
f: faces free of those on the bounding box
==== function eid=getintersecttri(tmppath) ====
eid=getintersecttri(tmppath)
get the IDs of self-intersecting elements from tetgen
call this when tetgen complains about self-intersection
input:
tmppath: working dir, use mwpath('') in most cases
output:
eid: an array of all intersecting surface elements,
one can read the corresponding node/elem by
[no,el]=readoff(mwpath('post_vmesh.off'));
==== function elem=delendelem(elem,mask) ====
elem=delendelem(elem,mask)
delete elements whose nodes are all edge nodes
input/output:
elem: input/output, surface/volumetric element list
mask: of length of node number, =0 for internal nodes, =1 for edge nodes
==== function [newnode,newface]=surfreorient(node,face) ====
[newnode,newface]=surfreorient(node,elem)
reorder nodes in a single closed surface to ensure the norms of all
triangles are pointing outward
input:
node: list of nodes
face: list of surface triangles (each row are indices of nodes of each triangle)
output:
newnode: the output node list, in most cases it equals node
newface: the face list with consistent ordering
=== # Mesh resampling and optimization ===
==== function [node,elem]=meshresample(v,f,keepratio) ====
[node,elem]=meshresample(v,f,keepratio)
resample mesh using CGAL mesh simplification utility
input:
v: list of nodes
f: list of surface elements (each row for each triangle)
keepratio: decimation rate, a number less than 1, as the percentage
of the elements after the sampling
output:
node: the node coordinates of the sampled surface mesh
elem: the element list of the sampled surface mesh
==== function [newno,newfc]=remeshsurf(node,face,opt) ====
[newno,newfc]=remeshsurf(node,face,opt)
remesh a triangular surface and the output is guaranteed to be
free of self-intersecting element. This function is similar to
meshresample, but it can both downsample or upsample a mesh
input:
node: list of nodes on the input suface mesh, 3 columns for x,y,z
face: list of trianglular elements on the surface, [n1,n2,n3,region_id]
opt: function parameters
opt.gridsize: resolution for the voxelization of the mesh
opt.closesize: if there are openings, set the closing diameter
opt.elemsize: the size of the element of the output surface
if opt is a scalar, it defines the elemsize and gridsize=opt/4
output:
newno: list of nodes on the resulting suface mesh, 3 columns for x,y,z
newfc: list of trianglular elements on the surface, [n1,n2,n3,region_id]
==== function p=smoothsurf(node,mask,conn,iter,useralpha,usermethod,userbeta) ====
p=smoothsurf(node,mask,conn,iter,useralpha,usermethod,userbeta)
smoothing a surface mesh
input:
node: node coordinates of a surface mesh
mask: flag whether a node is movable: 0 movable, 1 non-movable
if mask=[], it assumes all nodes are movable
conn: input, a cell structure of length size(node), conn{n}
contains a list of all neighboring node ID for node n,
this can be computed from meshconn function
iter: smoothing iteration number
useralpha: scaler, smoothing parameter, v(k+1)=(1-alpha)*v(k)+alpha*mean(neighbors)
usermethod: smoothing method, including 'laplacian','laplacianhc' and 'lowpass'
userbeta: scaler, smoothing parameter, for 'laplacianhc'
output:
p: output, the smoothed node coordinates
recommendations
Based on [Bade2006], 'Lowpass' method outperforms 'Laplacian-HC' in volume
preserving and both are significantly better than the standard Laplacian method
[Bade2006] R. Bade, H. Haase, B. Preim, "Comparison of Fundamental Mesh
Smoothing Algorithms for Medical Surface Models,"
Simulation and Visualization, pp. 289-304, 2006.
==== function [no,el,fc,nodemap]=sortmesh(origin,node,elem,ecol,face,fcol) ====
[no,el,fc]=sortmesh(origin,node,elem,face)
sort nodes and elements in a mesh so that the indexed
nodes and elements are closer to each order
(this may reduce cache-miss in a calculation)
input:
origin: sorting all nodes and elements with the distance and
angles wrt this location, if origin=[], it will be
node(1,:)
node: list of nodes
elem: list of elements (each row are indices of nodes of each element)
ecol: list of columns in elem to participate sorting
face: list of surface triangles (this can be omitted)
fcol: list of columns in face to participate sorting
output:
no: node coordinates in the sorted order
el: the element list in the sorted order
fc: the surface triangle list in the sorted order (can be ignored)
nodemap: the new node mapping order, no=node(nodemap,:)
==== function [newnode,newelem]=mergemesh(node,elem,varargin) ====
[newnode,newelem]=mergemesh(node,elem,varargin)
concatenate two or more tetrahedral meshes or triangular surfaces
input:
node: node coordinates, dimension (nn,3)
elem: tetrahedral element or triangle surface (nn,3) to (nn,5)
output:
newnode: the node coordinates after merging, dimension (nn,3)
newelem: tetrahedral element or surfaces after merging (nn,4) or (nhn,5)
note: you can call meshcheckrepair for the output newnode and
newelem to remove the duplicated nodes or elements. mergemesh does
detect self-intersecting elements when merging; to remove self-intersecting
elements, you need to use mergesurf().
example:
[node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
[node2,face2,elem2]=meshasphere([5 5 13.1],3,0.3,3);
[newnode,newelem]=mergemesh(node1,elem1,node2,elem2);
plotmesh(newnode,newelem);
figure;
[newnode,newface]=mergemesh(node1,face1,node2,face2);
plotmesh(newnode,newface,'x>5');
==== function [newnode,newelem,newface]=meshrefine(node,elem,varargin) ====
[newnode,newelem,newface]=meshrefine(node,elem,face,opt)
refine a tetrahedral mesh by adding new nodes or constraints
input parameters:
node: existing tetrahedral mesh node list
elem: existing tetrahedral element list
face: (optional) existing tetrahedral mesh surface triangle list
opt: options for mesh refinement:
if opt is a Nx3 array, opt is treated as a list of new nodes to
be inserted into the mesh. the new nodes must be located on the
surface or inside the original mesh. external nodes are
discarded, unless the opt.extcmdopt is specified.