Coverage for /usr/local/lib/python3.10/site-packages/spam/mesh/projection.py: 76%
58 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-24 17:26 +0000
1# Library of SPAM functions for projecting morphological field onto tetrahedral
2# meshes
3# Copyright (C) 2020 SPAM Contributors
4#
5# This program is free software: you can redistribute it and/or modify it
6# under the terms of the GNU General Public License as published by the Free
7# Software Foundation, either version 3 of the License, or (at your option)
8# any later version.
9#
10# This program is distributed in the hope that it will be useful, but WITHOUT
11# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13# more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
18"""
19 The ``spam.mesh.projection`` module offers a two functions that enables the construction of three dimensional unstructured meshes (four noded tetrahedra) that embbed heterogeneous data
20 (namely, **subvolumes** and interface **orientation**) based on eihter:
21 - the distance field of an **segmented images** (see ``spaM.filters.distanceField``),
22 - ideal **geometrical objects** (see ``spam.filters.objects``).
24 >>> import spam.mesh
25 >>> spam.mesh.projectField(mesh, fields)
26 >>> spam.mesh.projectObjects(mesh, objects)
28 NOTE
29 ----
31 Objects array conventions:
33 - Sphere: ``[radius, centerX, centerY, centerZ, phase]``
34 - Cylinder: ``[radius, centerX, centerY, centerZ, directionX, directionY, directionZ, phase]``
36 Distance field and orientation convention:
38 - The distance fields inside connected components have positive values
39 - The normal vectors are pointing outside (to negative values)
40 - The subvolumes returned correspond to the outside (negative part)
42"""
44import numpy
45from spambind.mesh.meshToolkit import projmorpho
48def _refactor_mesh(mesh):
49 """
50 Helper private function for projection.
51 Remove additional nodes and change connectivity matrix accordingly.
52 """
53 points = mesh["points"]
54 cells = mesh["cells"]
56 # DEBUG: Artificially add unused node
57 # points = numpy.insert(points, 10, [-1.0,-1.0,-1.0], axis=0)
58 # for i, node_number_x4 in enumerate(cells):
59 # for j, node_number in enumerate(node_number_x4):
60 # if cells[i, j] >= 10:
61 # cells[i, j] += 1
63 # test if unused nodes
64 cells_set = set(cells.ravel()) # ordered numbering in connectivity
65 n_points_used = len(cells_set)
67 if n_points_used != points.shape[0]:
68 print("Deleting points before projection")
70 # removing unused nodes
71 points_to_delete = []
72 n_deleted = 0
73 for new, old in enumerate(cells_set):
74 # check if holes in the connectivity matrix
75 if new != old - n_deleted:
76 points_to_delete.append(new + n_deleted)
77 n_deleted += 1
78 # print(new, old, n_deleted, old - n_deleted)
79 points = numpy.delete(points, points_to_delete, axis=0)
81 print("Renumbering connectivity matrix")
82 # renumbering connectivity matrix accordingly
83 new_node_numbering = {old: new + 1 for new, old in enumerate(cells_set)}
84 for i, node_number_x4 in enumerate(cells):
85 for j, node_number in enumerate(node_number_x4):
86 if new_node_numbering[cells[i, j]] != cells[i, j]:
87 cells[i, j] = new_node_numbering[cells[i, j]]
88 del cells_set
90 # check if cell number start by 1 or 0
91 if cells.min() == 0:
92 print("Shift connectivity by 1 to start at 1")
93 cells += 1
95 return points, cells
98def projectObjects(
99 mesh,
100 objects,
101 analyticalOrientation=True,
102 cutoff=1e-6,
103 writeConnectivity=None,
104 vtkMesh=False,
105):
106 """
107 This function project a set of objects onto an unstructured mesh.
108 Each objects can be attributed a phase.
110 Parameters
111 ----------
112 mesh: dict
113 Dictionary containing points coordinates and mesh connectivity.
115 .. code-block:: python
117 mesh = {
118 # numpy array of size number of nodes x 3
119 "points": points,
120 # numpy array of size number of elements x 4
121 "cells": connectivity,
122 }
124 objects: 2D array
125 The list of objects. Each line corresponds to an object encoded as follow:
127 .. code-block:: text
129 - spheres: radius, oZ, oY, oX, phase
130 - cylinders: radius, oZ, oY, oX, nZ, nY, nX, phase
132 analyticalOrientation: bool, default=True
133 Change how the interface's normals are computed:
135 - `True`: as the direction from the centroid of the tetrahedron to the closest highest value of the object distance field (ie, the center of a sphere, the axis of a cylinder).
136 - `False`: as the normals of the facets which points lie on the object surface.
138 cutoff: float, default=1e-6
139 Volume ratio below wich elements with interfaces are ignored and considered being part of a single phase.
141 writeConnectivity: string, default=None
142 When not None, it writes a text file called `writeConnectivity` the
143 list of node and the list of elements
144 which format is:
146 .. code-block:: text
148 COORdinates ! number of nodes
149 nodeId, 0, x, y, z
150 ...
152 ELEMents ! number of elemens
153 elemId, 0, elemType, n1, n2, n3, n4, subVol(-), normalX,
154 normalY, normalZ
155 ...
157 where:
159 - ``n1, n2, n3, n4`` is the connectivity (node numbers).
160 - ``subVol(-)`` is the sub volume of the terahedron inside the inclusion.
161 - ``normalX, normalY, normalZ`` are to components of the interface normal vector.
162 - ``elemType`` is the type of element. Their meaning depends on the their phase.
163 Correspondance can be found in the function output
164 after the key word **MATE** like:
166 .. code-block:: text
168 <projmorpho::set_materials
169 . field 1
170 . . MATE,1: background
171 . . MATE,2: phase 1
172 . . MATE,3: interface phase 1 with
173 background
174 . field 2
175 . . MATE,1: background
176 . . MATE,4: phase 2
177 . . MATE,5: interface phase 2 with
178 background
179 >
181 Sub volumes and interface vector are only relevant for interfaces. Default value is 0 and [1, 0, 0].
183 vtkMesh: bool, default=False
184 Writes the VTK of interpolated fields and materials. The files take the same name as ``writeConnectivity``.
186 Returns
187 -------
188 (nElem x 5) array
189 For each element: ``elemType, subVol(-), normalX, normalY, normalZ`` (see outputFile for details).
191 NOTE
192 ----
193 Only four noded tetrahedra meshes are supported.
195 >>> import spam.mesh
196 >>> # unstructured mesh
197 >>> points, cells = spam.mesh.createCuboid([1, 1, 1], 0.1)
198 >>> mesh = {"points": points, "cells": cells}
199 >>> # one centered sphere (0.5, 0.5, 0.5) of radius 0.2 (phase 1)
200 >>> sphere = [[0.2, 0.8, 0.1, 0.1, 1]]
201 >>> # projection
202 >>> materials = spam.mesh.projectObjects(mesh, sphere, writeConnectivity="mysphere", vtkMesh=True)
204 """
206 # init projmorpho
207 pr = projmorpho(
208 name="spam" if writeConnectivity is None else writeConnectivity,
209 cutoff=cutoff,
210 )
212 # STEP 1: objects
213 objects = numpy.array(objects)
215 # if objects.shape[1] == 5:
216 # swaps = [
217 # [1, 3], # swap ox <-> oz
218 # ]
219 # for a, b in swaps:
220 # # swap axis for origin point
221 # tmp = objects[:, a].copy()
222 # objects[:, a] = objects[:, b]
223 # objects[:, b] = tmp
224 # elif objects.shape[1] == 7:
225 # # swap axis ellipsoids
226 # tmp = objects[:, 0].copy()
227 # objects[:, 0] = objects[:, 2]
228 # objects[:, 2] = tmp
229 # tmp = objects[:, 3].copy()
230 # objects[:, 3] = objects[:, 5]
231 # objects[:, 5] = tmp
232 # elif objects.shape[1] == 8: # cylinders
233 # swaps = [
234 # [1, 3], # swap ox <-> oz
235 # [4, 6], # swap nx <-> nz
236 # ]
237 # for a, b in swaps:
238 # # swap axis for origin point
239 # tmp = objects[:, a].copy()
240 # objects[:, a] = objects[:, b]
241 # objects[:, b] = tmp
243 pr.setObjects(objects)
245 # STEP 2: Mesh
246 points, cells = _refactor_mesh(mesh)
247 pr.setMesh(points, cells.ravel())
249 # STEP 3: projection
250 pr.computeFieldFromObjects()
251 pr.projection(analytical_orientation=analyticalOrientation)
253 # STEP 4: write VTK
254 if writeConnectivity:
255 pr.writeFEAP()
256 if vtkMesh:
257 pr.writeVTK()
258 pr.writeInterfacesVTK()
260 return numpy.array(pr.getMaterials())
263def projectField(mesh, fields, thresholds=[0.0], cutoff=1e-6, writeConnectivity=None, vtkMesh=False):
264 """
265 This function project a set of distance fields onto an unstructured mesh.
266 Each distance field corresponds to a phase and the interface between
267 the two phases is set by the thresholds.
269 Parameters
270 ----------
271 mesh: dict
272 Dictionary containing points coordinates and mesh connectivity.
274 .. code-block:: python
276 mesh = {
277 # numpy array of size number of nodes x 3
278 "points": points,
279 # numpy array of size number of elements x 4
280 "cells": connectivity,
281 }
283 fields: list of dicts
284 The fields should be continuous (like a distance
285 field) for a better projection. They are discretised over a regular
286 mesh (lexicographical order) and eahc one corresponds to a phase.
287 The dictionary containing the fields data is defined as follow
289 .. code-block:: python
291 fields = {
292 # coordinates of the origin of the field (3 x 1)
293 "origin": origin,
294 # lengths of fields domain of definition (3 x 1)
295 "lengths": lengths,
296 # list of fields values (list of 3D arrays)
297 "values": [phase1, phase2],
298 }
300 thresholds: list of floats, default=[0]
301 The list of thresholds.
303 cutoff: float, default=1e-6
304 Volume ratio below wich elements with interfaces are ignored and considered being part of a single phase.
306 writeConnectivity: string, default=None
307 When not None, it writes a text file called `writeConnectivity` the
308 list of node and the list of elements
309 which format is:
311 .. code-block:: text
313 COORdinates ! number of nodes
314 nodeId, 0, x, y, z
315 ...
317 ELEMents ! number of elemens
318 elemId, 0, elemType, n1, n2, n3, n4, subVol(-), normalX,
319 normalY, normalZ
320 ...
322 where:
324 - ``n1, n2, n3, n4`` is the connectivity (node numbers).
325 - ``subVol(-)`` is the sub volume of the terahedron inside the inclusion.
326 - ``normalX, normalY, normalZ`` are to components of the interface normal vector.
327 - ``elemType`` is the type of element. Their meaning depends on the their phase.
328 Correspondance can be found in the function output
329 after the key word **MATE** like:
331 .. code-block:: text
333 <projmorpho::set_materials
334 . field 1
335 . . MATE,1: background
336 . . MATE,2: phase 1
337 . . MATE,3: interface phase 1 with
338 background
339 . field 2
340 . . MATE,1: background
341 . . MATE,4: phase 2
342 . . MATE,5: interface phase 2 with
343 background
344 >
346 Sub volumes and interface vector are only relevant for interfaces. Default value is 0 and [1, 0, 0].
348 vtkMesh: bool, default=False
349 Writes the VTK of interpolated fields and materials. The files take the same name as ``writeConnectivity``.
351 Returns
352 -------
353 (nElem x 5) array
354 For each element: ``elemType, subVol(-), normalX, normalY, normalZ`` (see outputFile for details).
356 NOTE
357 ----
358 Only four noded tetrahedra meshes are supported.
360 >>> import spam.mesh
361 >>> import spam.kalisphera
362 >>> # unstructured mesh
363 >>> points, cells = spam.mesh.createCuboid([1, 1, 1], 0.1)
364 >>> mesh = {"points": points, "cells": cells}
365 >>> # create an image of a sphere and its distance field
366 >>> sphereBinary = numpy.zeros([101] * 3, dtype=float)
367 >>> spam.kalisphera.makeSphere(sphereBinary, [50.5] * 3, 20)
368 >>> sphereField = spam.mesh.distanceField(one_sphere_binary.astype(bool).astype(int))
369 >>> # format field object
370 >>> fields = {
371 >>> # coordinates of the origin of the field (3 x 1)
372 >>> "origin": [0] * 3,
373 >>> # lengths of fields domain of definition (3 x 1)
374 >>> "lengths": [1] * 3,
375 >>> # list of fields
376 >>> "values": [sphereField],
377 >>> }
378 >>> # projection
379 >>> materials = spam.mesh.projectField(mesh, fields, writeConnectivity="mysphere", vtkMesh=True)
381 """
382 # init projmorpho
383 pr = projmorpho(
384 name="spam" if writeConnectivity is None else writeConnectivity,
385 thresholds=thresholds,
386 cutoff=cutoff,
387 )
389 # STEP 1: Fields
391 # origin
392 # origin = [_ for _ in reversed(fields["origin"])]
393 origin = [_ for _ in fields["origin"]]
394 # nPoints
395 # nPoints = [n for n in reversed(fields["values"][0].shape)]
396 nPoints = [n for n in fields["values"][0].shape]
397 # field size
398 # lengths = [_ for _ in reversed(fields["lengths"])]
399 lengths = [_ for _ in fields["lengths"]]
400 # ravel fields values (ravel in F direction to swap zyx to xyz)
401 values = [v.flatten(order="F") for v in fields["values"]]
402 # values = [v.flatten(order='C') for v in fields["values"]]
404 pr.setImage(values, lengths, nPoints, origin)
406 # STEP 2: Mesh
407 points, cells = _refactor_mesh(mesh)
408 pr.setMesh(points, cells.ravel())
410 # STEP 3: projection
411 pr.computeFieldFromImages()
412 pr.projection()
414 # STEP 4: write VTK
415 if writeConnectivity:
416 pr.writeFEAP()
417 if vtkMesh:
418 pr.writeVTK()
419 pr.writeInterfacesVTK()
421 return numpy.array(pr.getMaterials())