Line data Source code
1 : #
2 : # --- DEPRECATED!!! --- (see MP example 2023)
3 : #
4 :
5 : using LinearAlgebra
6 : using Statistics
7 :
8 : # TODO: use geometry(), defin helper for "marking" fearures & sizing fields
9 : # TODO: provide/store context in a IsometricRemeshing struct
10 :
11 0 : function facenormals!(mesh; pos::Symbol=:x, fnrm::Symbol=:n, farea::Symbol=:a)
12 0 : x = vattr(mesh, pos)
13 0 : fn = fattr(mesh, fnrm)
14 0 : area = fattr(mesh, farea)
15 :
16 0 : for f in faces(mesh)
17 0 : h = halfedge(mesh, f)
18 0 : v0, vc = vertices(mesh, h)
19 0 : v1 = destination(mesh, next(mesh, h))
20 :
21 0 : u = x[v0] - x[vc]
22 0 : v = x[v1] - x[vc]
23 :
24 0 : nrm = cross(u, v)
25 0 : ar2 = norm(nrm)
26 0 : area[f] = ar2 / 2
27 :
28 0 : if ar2 > eps(ar2)
29 0 : fn[f] = nrm / ar2
30 : else
31 0 : fn[f] = zeros(typeof(nrm))
32 : end
33 0 : end
34 : end
35 :
36 : # TODO: face areas
37 :
38 0 : function normals!(mesh;
39 : pos::Symbol=:x, vnrm::Symbol=:n,
40 : fnrm::Symbol=:n, farea::Symbol=:a)
41 0 : x = vattr(mesh, pos)
42 0 : vn = vattr(mesh, vnrm)
43 0 : fn = fattr(mesh, fnrm)
44 0 : area = fattr(mesh, farea)
45 :
46 0 : facenormals!(mesh; pos, fnrm, farea)
47 :
48 0 : for v in vertices(mesh)
49 0 : nrm = sum(fn[f] * area[f] for f in fan(mesh, v))
50 0 : vn[v] = nrm / norm(nrm)
51 0 : end
52 : end
53 :
54 0 : function tangentialsmoothing!(mesh, k::Int=1, μ=0.25;
55 : pos::Symbol=:x, vnrm::Symbol=:n,
56 : fnrm::Symbol=:n, farea::Symbol=:a)
57 0 : @info "tangential smoothing μ=$μ, k=$k"
58 :
59 0 : normals!(mesh; pos, vnrm, fnrm)
60 :
61 0 : for _ in 1:k
62 0 : _tangentialsmoothing!(mesh, μ; pos, vnrm)
63 0 : end
64 : end
65 :
66 0 : function _tangentialsmoothing!(mesh, μ;
67 : pos::Symbol=:x, vnrm::Symbol=:n,
68 : farea::Symbol=:a)
69 0 : x = vattr(mesh, pos)
70 0 : vn = vattr(mesh, vnrm)
71 :
72 : # TODO: store update vector?
73 0 : @assert allfinite(mesh)
74 :
75 0 : for v in vertices(mesh)
76 0 : u = (mean(x[vn] for vn in ring(mesh, v)) - x[v]) * μ
77 0 : n = vn[v]
78 :
79 0 : if !isfinite(sum(n))
80 0 : @show v
81 0 : @show x[v]
82 0 : @show degree(mesh, v)
83 0 : @show [fattr(mesh, :n)[f] for f in fan(mesh, v)]
84 0 : @show [fattr(mesh, :a)[f] for f in fan(mesh, v)]
85 :
86 0 : break
87 : end
88 :
89 0 : if !isfinite(sum(u))
90 0 : @show v
91 0 : @show [vn for vn in ring(mesh, v)]
92 0 : @show [isused(mesh, vn) for vn in ring(mesh, v)]
93 0 : @show [x[vn] for vn in ring(mesh, v)]
94 0 : @show x[v]
95 0 : @show u
96 0 : @show mean(x[vn] for vn in ring(mesh, v))
97 0 : @show sum(x[vn] for vn in ring(mesh, v))
98 0 : @show degree(mesh, v)
99 0 : break
100 : end
101 :
102 0 : x[v] += u - dot(u, n) * n
103 0 : end
104 : end
105 :
106 : # TODO: projection (w/ distance queries)
107 :
108 : #-----------------------------------------------------------------------------
109 :
110 : # TODO: pass arguments
111 :
112 : # TODO: check isfinite
113 :
114 0 : function allfinite(mesh)
115 0 : x = vattr(mesh, :x)
116 :
117 0 : all(isfinite(sum(x[v])) for v in vertices(mesh))
118 : end
119 :
120 : ###
121 :
122 0 : remesh!(mesh;
123 : k=32,
124 : pos::Symbol=:x, vnrm::Symbol=:n,
125 : fnrm::Symbol=:n, farea::Symbol=:a) =
126 : remesh!(mesh, mean(edgelengths(mesh, pos)); k, pos, vnrm, fnrm, farea)
127 :
128 0 : function remesh!(mesh, len;
129 : k=32,
130 : pos::Symbol=:x, vnrm::Symbol=:n,
131 : fnrm::Symbol=:n, farea::Symbol=:a)
132 : # TODO: no special cases for feature edges, boundaries, adaptive sizing
133 :
134 0 : for _ in 1:k
135 0 : splitlongedges!(mesh, 4len/3, pos)
136 0 : @assert allfinite(mesh)
137 0 : collapseshortedges!(mesh, 4len/5, pos)
138 0 : @assert allfinite(mesh)
139 0 : balancevalence!(mesh)
140 0 : tangentialsmoothing!(mesh; pos, vnrm, fnrm, farea)
141 0 : @assert allfinite(mesh)
142 : # TODO: projection
143 0 : end
144 : end
145 :
146 : """
147 : remesh!(mesh[, len=meanedgelength(mesh);
148 : k=32, pos=:x, vnrm=:n, fnrm=:n, farea=:a)
149 :
150 : Apply `k` iterations of isotropic remeshing (w/o projection onto the mesh).
151 :
152 : !!! warning
153 : This method is intended only for tests and benchmarks!
154 : """ remesh!
155 :
156 0 : function splitlongedges!(mesh, lmax, pos::Symbol=:x)
157 0 : @info "split long edges"
158 0 : n = 0
159 0 : x = vattr(mesh, pos)
160 0 : for _ in 1:10
161 0 : done = true
162 0 : for e in edges(mesh)
163 0 : v01 = vertices(mesh, e)
164 0 : (edgelength(mesh, v01, x) <= lmax) && continue
165 :
166 0 : done = false
167 0 : _, vnew = split!(mesh, e)
168 0 : x[vnew] = (x[v01[1]] + x[v01[2]]) / 2
169 0 : n += 1
170 0 : end
171 :
172 0 : done && break
173 0 : end
174 0 : @info "$n splits"
175 : end
176 :
177 0 : function collapseshortedges!(mesh, lmin, pos::Symbol=:x)
178 0 : @info "collapse short edges"
179 0 : n = 0
180 0 : x = vattr(mesh, pos)
181 0 : for _ in 1:10
182 0 : done = true
183 0 : for e in edges(mesh)
184 0 : isused(mesh, e) || continue
185 :
186 0 : (edgelength(mesh, e, x) >= lmin) && continue
187 :
188 0 : hs = halfedges(mesh, e)
189 0 : yes = (maycollapse(mesh, h) for h in hs)
190 0 : any(yes) || continue
191 :
192 0 : h =
193 : if all(yes)
194 0 : ds = tuple((degree(mesh, destination(mesh, h)) for h in hs)...)
195 0 : ds[1] >= ds[2] ? hs[1] : hs[2]
196 : else
197 0 : yes[1] ? hs[1] : hs[2]
198 : end
199 :
200 0 : done = false
201 0 : collapse!(mesh, h)
202 :
203 0 : n += 1
204 0 : end
205 :
206 0 : done && break
207 0 : end
208 0 : compact!(mesh; stable=false)
209 0 : @info "$n collapses & compact"
210 : end
211 :
212 0 : function balancevalence!(mesh)
213 0 : @info "balance valance"
214 0 : n = 0
215 :
216 0 : for _ in 1:10
217 0 : done = true
218 :
219 0 : for e in edges(mesh)
220 0 : mayflip(mesh, e) || continue
221 :
222 0 : h1, h2 = halfedges(mesh, e)
223 0 : vl = destination(mesh, h1)
224 0 : vr = destination(mesh, h2)
225 0 : vu = destination(mesh, next(mesh, h1))
226 0 : vd = destination(mesh, next(mesh, h2))
227 :
228 0 : ds = tuple((degree(mesh, v) for v in (vl, vr, vu, vd))...)
229 0 : before = sum((d-6)^2 for d in ds)
230 :
231 0 : ds = (ds[1]-1, ds[2]-1, ds[3]+1, ds[4]+1)
232 0 : after = sum((d-6)^2 for d in ds)
233 :
234 0 : (before <= after) && continue
235 :
236 0 : flip!(mesh, e)
237 0 : done = false
238 0 : n += 1
239 0 : end
240 :
241 0 : done && break
242 0 : end
243 0 : @info "$n flips"
244 0 : n
245 : end
|
