1# Contributed by Seva Alekseyev <sevaa@nih.gov> with National Institutes of Health, 2016 2# Cura is released under the terms of the LGPLv3 or higher. 3 4from math import pi, sin, cos, sqrt 5from typing import Dict 6 7import numpy 8 9from UM.Job import Job 10from UM.Logger import Logger 11from UM.Math.Matrix import Matrix 12from UM.Math.Vector import Vector 13from UM.Mesh.MeshBuilder import MeshBuilder 14from UM.Mesh.MeshReader import MeshReader 15from cura.Scene.CuraSceneNode import CuraSceneNode as SceneNode 16 17MYPY = False 18try: 19 if not MYPY: 20 import xml.etree.cElementTree as ET 21except ImportError: 22 import xml.etree.ElementTree as ET 23 24# TODO: preserve the structure of scenes that contain several objects 25# Use CADPart, for example, to distinguish between separate objects 26 27DEFAULT_SUBDIV = 16 # Default subdivision factor for spheres, cones, and cylinders 28EPSILON = 0.000001 29 30 31class Shape: 32 # Expects verts in MeshBuilder-ready format, as a n by 3 mdarray 33 # with vertices stored in rows 34 def __init__(self, verts, faces, index_base, name): 35 self.verts = verts 36 self.faces = faces 37 # Those are here for debugging purposes only 38 self.index_base = index_base 39 self.name = name 40 41 42class X3DReader(MeshReader): 43 def __init__(self) -> None: 44 super().__init__() 45 self._supported_extensions = [".x3d"] 46 self._namespaces = {} # type: Dict[str, str] 47 48 # Main entry point 49 # Reads the file, returns a SceneNode (possibly with nested ones), or None 50 def _read(self, file_name): 51 try: 52 self.defs = {} 53 self.shapes = [] 54 55 tree = ET.parse(file_name) 56 xml_root = tree.getroot() 57 58 if xml_root.tag != "X3D": 59 return None 60 61 scale = 1000 # Default X3D unit it one meter, while Cura's is one millimeters 62 if xml_root[0].tag == "head": 63 for head_node in xml_root[0]: 64 if head_node.tag == "unit" and head_node.attrib.get("category") == "length": 65 scale *= float(head_node.attrib["conversionFactor"]) 66 break 67 xml_scene = xml_root[1] 68 else: 69 xml_scene = xml_root[0] 70 71 if xml_scene.tag != "Scene": 72 return None 73 74 self.transform = Matrix() 75 self.transform.setByScaleFactor(scale) 76 self.index_base = 0 77 78 # Traverse the scene tree, populate the shapes list 79 self.processChildNodes(xml_scene) 80 81 if self.shapes: 82 builder = MeshBuilder() 83 builder.setVertices(numpy.concatenate([shape.verts for shape in self.shapes])) 84 builder.setIndices(numpy.concatenate([shape.faces for shape in self.shapes])) 85 builder.calculateNormals() 86 builder.setFileName(file_name) 87 mesh_data = builder.build() 88 89 # Manually try and get the extents of the mesh_data. This should prevent nasty NaN issues from 90 # leaving the reader. 91 mesh_data.getExtents() 92 93 node = SceneNode() 94 node.setMeshData(mesh_data) 95 node.setSelectable(True) 96 node.setName(file_name) 97 98 else: 99 return None 100 101 except Exception: 102 Logger.logException("e", "Exception in X3D reader") 103 return None 104 105 return node 106 107 # ------------------------- XML tree traversal 108 109 def processNode(self, xml_node): 110 xml_node = self.resolveDefUse(xml_node) 111 if xml_node is None: 112 return 113 114 tag = xml_node.tag 115 if tag in ("Group", "StaticGroup", "CADAssembly", "CADFace", "CADLayer", "Collision"): 116 self.processChildNodes(xml_node) 117 if tag == "CADPart": 118 self.processTransform(xml_node) # TODO: split the parts 119 elif tag == "LOD": 120 self.processNode(xml_node[0]) 121 elif tag == "Transform": 122 self.processTransform(xml_node) 123 elif tag == "Shape": 124 self.processShape(xml_node) 125 126 127 def processShape(self, xml_node): 128 # Find the geometry and the appearance inside the Shape 129 geometry = appearance = None 130 for sub_node in xml_node: 131 if sub_node.tag == "Appearance" and not appearance: 132 appearance = self.resolveDefUse(sub_node) 133 elif sub_node.tag in self.geometry_importers and not geometry: 134 geometry = self.resolveDefUse(sub_node) 135 136 # TODO: appearance is completely ignored. At least apply the material color... 137 if not geometry is None: 138 try: 139 self.verts = self.faces = [] # Safeguard 140 self.geometry_importers[geometry.tag](self, geometry) 141 m = self.transform.getData() 142 verts = m.dot(self.verts)[:3].transpose() 143 144 self.shapes.append(Shape(verts, self.faces, self.index_base, geometry.tag)) 145 self.index_base += len(verts) 146 147 except Exception: 148 Logger.logException("e", "Exception in X3D reader while reading %s", geometry.tag) 149 150 # Returns the referenced node if the node has USE, the same node otherwise. 151 # May return None is USE points at a nonexistent node 152 # In X3DOM, when both DEF and USE are in the same node, DEF is ignored. 153 # Big caveat: XML element objects may evaluate to boolean False!!! 154 # Don't ever use "if node:", use "if not node is None:" instead 155 def resolveDefUse(self, node): 156 USE = node.attrib.get("USE") 157 if USE: 158 return self.defs.get(USE, None) 159 160 DEF = node.attrib.get("DEF") 161 if DEF: 162 self.defs[DEF] = node 163 return node 164 165 def processChildNodes(self, node): 166 for c in node: 167 self.processNode(c) 168 Job.yieldThread() 169 170 # Since this is a grouping node, will recurse down the tree. 171 # According to the spec, the final transform matrix is: 172 # T * C * R * SR * S * -SR * -C 173 # Where SR corresponds to the rotation matrix to scaleOrientation 174 # C and SR are rather exotic. S, slightly less so. 175 def processTransform(self, node): 176 rot = readRotation(node, "rotation", (0, 0, 1, 0)) # (angle, axisVactor) tuple 177 trans = readVector(node, "translation", (0, 0, 0)) # Vector 178 scale = readVector(node, "scale", (1, 1, 1)) # Vector 179 center = readVector(node, "center", (0, 0, 0)) # Vector 180 scale_orient = readRotation(node, "scaleOrientation", (0, 0, 1, 0)) # (angle, axisVactor) tuple 181 182 # Store the previous transform; in Cura, the default matrix multiplication is in place 183 prev = Matrix(self.transform.getData()) # It's deep copy, I've checked 184 185 # The rest of transform manipulation will be applied in place 186 got_center = (center.x != 0 or center.y != 0 or center.z != 0) 187 188 T = self.transform 189 if trans.x != 0 or trans.y != 0 or trans.z != 0: 190 T.translate(trans) 191 if got_center: 192 T.translate(center) 193 if rot[0] != 0: 194 T.rotateByAxis(*rot) 195 if scale.x != 1 or scale.y != 1 or scale.z != 1: 196 got_scale_orient = scale_orient[0] != 0 197 if got_scale_orient: 198 T.rotateByAxis(*scale_orient) 199 # No scale by vector in place operation in UM 200 S = Matrix() 201 S.setByScaleVector(scale) 202 T.multiply(S) 203 if got_scale_orient: 204 T.rotateByAxis(-scale_orient[0], scale_orient[1]) 205 if got_center: 206 T.translate(-center) 207 208 self.processChildNodes(node) 209 self.transform = prev 210 211 # ------------------------- Geometry importers 212 # They are supposed to fill the self.verts and self.faces arrays, the caller will do the rest 213 214 # Primitives 215 216 def processGeometryBox(self, node): 217 (dx, dy, dz) = readFloatArray(node, "size", [2, 2, 2]) 218 dx /= 2 219 dy /= 2 220 dz /= 2 221 self.reserveFaceAndVertexCount(12, 8) 222 223 # xz plane at +y, ccw 224 self.addVertex(dx, dy, dz) 225 self.addVertex(-dx, dy, dz) 226 self.addVertex(-dx, dy, -dz) 227 self.addVertex(dx, dy, -dz) 228 # xz plane at -y 229 self.addVertex(dx, -dy, dz) 230 self.addVertex(-dx, -dy, dz) 231 self.addVertex(-dx, -dy, -dz) 232 self.addVertex(dx, -dy, -dz) 233 234 self.addQuad(0, 1, 2, 3) # +y 235 self.addQuad(4, 0, 3, 7) # +x 236 self.addQuad(7, 3, 2, 6) # -z 237 self.addQuad(6, 2, 1, 5) # -x 238 self.addQuad(5, 1, 0, 4) # +z 239 self.addQuad(7, 6, 5, 4) # -y 240 241 # The sphere is subdivided into nr rings and ns segments 242 def processGeometrySphere(self, node): 243 r = readFloat(node, "radius", 0.5) 244 subdiv = readIntArray(node, "subdivision", None) 245 if subdiv: 246 if len(subdiv) == 1: 247 nr = ns = subdiv[0] 248 else: 249 (nr, ns) = subdiv 250 else: 251 nr = ns = DEFAULT_SUBDIV 252 253 lau = pi / nr # Unit angle of latitude (rings) for the given tesselation 254 lou = 2 * pi / ns # Unit angle of longitude (segments) 255 256 self.reserveFaceAndVertexCount(ns*(nr*2 - 2), 2 + (nr - 1)*ns) 257 258 # +y and -y poles 259 self.addVertex(0, r, 0) 260 self.addVertex(0, -r, 0) 261 262 # The non-polar vertices go from x=0, negative z plane counterclockwise - 263 # to -x, to +z, to +x, back to -z 264 for ring in range(1, nr): 265 for seg in range(ns): 266 self.addVertex(-r*sin(lou * seg) * sin(lau * ring), 267 r*cos(lau * ring), 268 -r*cos(lou * seg) * sin(lau * ring)) 269 270 vb = 2 + (nr - 2) * ns # First vertex index for the bottom cap 271 272 # Faces go in order: top cap, sides, bottom cap. 273 # Sides go by ring then by segment. 274 275 # Caps 276 # Top cap face vertices go in order: down right up 277 # (starting from +y pole) 278 # Bottom cap goes: up left down (starting from -y pole) 279 for seg in range(ns): 280 self.addTri(0, seg + 2, (seg + 1) % ns + 2) 281 self.addTri(1, vb + (seg + 1) % ns, vb + seg) 282 283 # Sides 284 # Side face vertices go in order: down right upleft, downright up left 285 for ring in range(nr - 2): 286 tvb = 2 + ring * ns 287 # First vertex index for the top edge of the ring 288 bvb = tvb + ns 289 # First vertex index for the bottom edge of the ring 290 for seg in range(ns): 291 nseg = (seg + 1) % ns 292 self.addQuad(tvb + seg, bvb + seg, bvb + nseg, tvb + nseg) 293 294 def processGeometryCone(self, node): 295 r = readFloat(node, "bottomRadius", 1) 296 height = readFloat(node, "height", 2) 297 bottom = readBoolean(node, "bottom", True) 298 side = readBoolean(node, "side", True) 299 n = readInt(node, "subdivision", DEFAULT_SUBDIV) 300 301 d = height / 2 302 angle = 2 * pi / n 303 304 self.reserveFaceAndVertexCount((n if side else 0) + (n-2 if bottom else 0), n+1) 305 306 # Vertex 0 is the apex, vertices 1..n are the bottom 307 self.addVertex(0, d, 0) 308 for i in range(n): 309 self.addVertex(-r * sin(angle * i), -d, -r * cos(angle * i)) 310 311 # Side face vertices go: up down right 312 if side: 313 for i in range(n): 314 self.addTri(1 + (i + 1) % n, 0, 1 + i) 315 if bottom: 316 for i in range(2, n): 317 self.addTri(1, i, i+1) 318 319 def processGeometryCylinder(self, node): 320 r = readFloat(node, "radius", 1) 321 height = readFloat(node, "height", 2) 322 bottom = readBoolean(node, "bottom", True) 323 side = readBoolean(node, "side", True) 324 top = readBoolean(node, "top", True) 325 n = readInt(node, "subdivision", DEFAULT_SUBDIV) 326 327 nn = n * 2 328 angle = 2 * pi / n 329 hh = height/2 330 331 self.reserveFaceAndVertexCount((nn if side else 0) + (n - 2 if top else 0) + (n - 2 if bottom else 0), nn) 332 333 # The seam is at x=0, z=-r, vertices go ccw - 334 # to pos x, to neg z, to neg x, back to neg z 335 for i in range(n): 336 rs = -r * sin(angle * i) 337 rc = -r * cos(angle * i) 338 self.addVertex(rs, hh, rc) 339 self.addVertex(rs, -hh, rc) 340 341 if side: 342 for i in range(n): 343 ni = (i + 1) % n 344 self.addQuad(ni * 2 + 1, ni * 2, i * 2, i * 2 + 1) 345 346 for i in range(2, nn-3, 2): 347 if top: 348 self.addTri(0, i, i+2) 349 if bottom: 350 self.addTri(1, i+1, i+3) 351 352 # Semi-primitives 353 354 def processGeometryElevationGrid(self, node): 355 dx = readFloat(node, "xSpacing", 1) 356 dz = readFloat(node, "zSpacing", 1) 357 nx = readInt(node, "xDimension", 0) 358 nz = readInt(node, "zDimension", 0) 359 height = readFloatArray(node, "height", False) 360 ccw = readBoolean(node, "ccw", True) 361 362 if nx <= 0 or nz <= 0 or len(height) < nx*nz: 363 return # That's weird, the wording of the standard suggests grids with zero quads are somehow valid 364 365 self.reserveFaceAndVertexCount(2*(nx-1)*(nz-1), nx*nz) 366 367 for z in range(nz): 368 for x in range(nx): 369 self.addVertex(x * dx, height[z*nx + x], z * dz) 370 371 for z in range(1, nz): 372 for x in range(1, nx): 373 self.addTriFlip((z - 1)*nx + x - 1, z*nx + x, (z - 1)*nx + x, ccw) 374 self.addTriFlip((z - 1)*nx + x - 1, z*nx + x - 1, z*nx + x, ccw) 375 376 def processGeometryExtrusion(self, node): 377 ccw = readBoolean(node, "ccw", True) 378 begin_cap = readBoolean(node, "beginCap", True) 379 end_cap = readBoolean(node, "endCap", True) 380 cross = readFloatArray(node, "crossSection", (1, 1, 1, -1, -1, -1, -1, 1, 1, 1)) 381 cross = [(cross[i], cross[i+1]) for i in range(0, len(cross), 2)] 382 spine = readFloatArray(node, "spine", (0, 0, 0, 0, 1, 0)) 383 spine = [(spine[i], spine[i+1], spine[i+2]) for i in range(0, len(spine), 3)] 384 orient = readFloatArray(node, "orientation", None) 385 if orient: 386 # This converts X3D's axis/angle rotation to a 3x3 numpy matrix 387 def toRotationMatrix(rot): 388 (x, y, z) = rot[:3] 389 a = rot[3] 390 s = sin(a) 391 c = cos(a) 392 t = 1-c 393 return numpy.array(( 394 (x * x * t + c, x * y * t - z*s, x * z * t + y * s), 395 (x * y * t + z*s, y * y * t + c, y * z * t - x * s), 396 (x * z * t - y * s, y * z * t + x * s, z * z * t + c))) 397 398 orient = [toRotationMatrix(orient[i:i+4]) if orient[i+3] != 0 else None for i in range(0, len(orient), 4)] 399 400 scale = readFloatArray(node, "scale", None) 401 if scale: 402 scale = [numpy.array(((scale[i], 0, 0), (0, 1, 0), (0, 0, scale[i+1]))) 403 if scale[i] != 1 or scale[i+1] != 1 else None for i in range(0, len(scale), 2)] 404 405 406 # Special treatment for the closed spine and cross section. 407 # Let's save some memory by not creating identical but distinct vertices; 408 # later we'll introduce conditional logic to link the last vertex with 409 # the first one where necessary. 410 crossClosed = cross[0] == cross[-1] 411 if crossClosed: 412 cross = cross[:-1] 413 nc = len(cross) 414 cross = [numpy.array((c[0], 0, c[1])) for c in cross] 415 ncf = nc if crossClosed else nc - 1 416 # Face count along the cross; for closed cross, it's the same as the 417 # respective vertex count 418 419 spine_closed = spine[0] == spine[-1] 420 if spine_closed: 421 spine = spine[:-1] 422 ns = len(spine) 423 spine = [Vector(*s) for s in spine] 424 nsf = ns if spine_closed else ns - 1 425 426 # This will be used for fallback, where the current spine point joins 427 # two collinear spine segments. No need to recheck the case of the 428 # closed spine/last-to-first point juncture; if there's an angle there, 429 # it would kick in on the first iteration of the main loop by spine. 430 def findFirstAngleNormal(): 431 for i in range(1, ns - 1): 432 spt = spine[i] 433 z = (spine[i + 1] - spt).cross(spine[i - 1] - spt) 434 if z.length() > EPSILON: 435 return z 436 # All the spines are collinear. Fallback to the rotated source 437 # XZ plane. 438 # TODO: handle the situation where the first two spine points match 439 if len(spine) < 2: 440 return Vector(0, 0, 1) 441 v = spine[1] - spine[0] 442 orig_y = Vector(0, 1, 0) 443 orig_z = Vector(0, 0, 1) 444 if v.cross(orig_y).length() > EPSILON: 445 # Spine at angle with global y - rotate the z accordingly 446 a = v.cross(orig_y) # Axis of rotation to get to the Z 447 (x, y, z) = a.normalized().getData() 448 s = a.length()/v.length() 449 c = sqrt(1-s*s) 450 t = 1-c 451 m = numpy.array(( 452 (x * x * t + c, x * y * t + z*s, x * z * t - y * s), 453 (x * y * t - z*s, y * y * t + c, y * z * t + x * s), 454 (x * z * t + y * s, y * z * t - x * s, z * z * t + c))) 455 orig_z = Vector(*m.dot(orig_z.getData())) 456 return orig_z 457 458 self.reserveFaceAndVertexCount(2*nsf*ncf + (nc - 2 if begin_cap else 0) + (nc - 2 if end_cap else 0), ns*nc) 459 460 z = None 461 for i, spt in enumerate(spine): 462 if (i > 0 and i < ns - 1) or spine_closed: 463 snext = spine[(i + 1) % ns] 464 sprev = spine[(i - 1 + ns) % ns] 465 y = snext - sprev 466 vnext = snext - spt 467 vprev = sprev - spt 468 try_z = vnext.cross(vprev) 469 # Might be zero, then all kinds of fallback 470 if try_z.length() > EPSILON: 471 if z is not None and try_z.dot(z) < 0: 472 try_z = -try_z 473 z = try_z 474 elif not z: # No z, and no previous z. 475 # Look ahead, see if there's at least one point where 476 # spines are not collinear. 477 z = findFirstAngleNormal() 478 elif i == 0: # And non-crossed 479 snext = spine[i + 1] 480 y = snext - spt 481 z = findFirstAngleNormal() 482 else: # last point and not crossed 483 sprev = spine[i - 1] 484 y = spt - sprev 485 # If there's more than one point in the spine, z is already set. 486 # One point in the spline is an error anyway. 487 488 z = z.normalized() 489 y = y.normalized() 490 x = y.cross(z) # Already normalized 491 m = numpy.array(((x.x, y.x, z.x), (x.y, y.y, z.y), (x.z, y.z, z.z))) 492 493 # Columns are the unit vectors for the xz plane for the cross-section 494 if orient: 495 mrot = orient[i] if len(orient) > 1 else orient[0] 496 if not mrot is None: 497 m = m.dot(mrot) # Tested against X3DOM, the result matches, still not sure :( 498 499 if scale: 500 mscale = scale[i] if len(scale) > 1 else scale[0] 501 if not mscale is None: 502 m = m.dot(mscale) 503 504 # First the cross-section 2-vector is scaled, 505 # then rotated (which may make it a 3-vector), 506 # then applied to the xz plane unit vectors 507 508 sptv3 = numpy.array(spt.getData()[:3]) 509 for cpt in cross: 510 v = sptv3 + m.dot(cpt) 511 self.addVertex(*v) 512 513 if begin_cap: 514 self.addFace([x for x in range(nc - 1, -1, -1)], ccw) 515 516 # Order of edges in the face: forward along cross, forward along spine, 517 # backward along cross, backward along spine, flipped if now ccw. 518 # This order is assumed later in the texture coordinate assignment; 519 # please don't change without syncing. 520 521 for s in range(ns - 1): 522 for c in range(ncf): 523 self.addQuadFlip(s * nc + c, s * nc + (c + 1) % nc, 524 (s + 1) * nc + (c + 1) % nc, (s + 1) * nc + c, ccw) 525 526 if spine_closed: 527 # The faces between the last and the first spine points 528 b = (ns - 1) * nc 529 for c in range(ncf): 530 self.addQuadFlip(b + c, b + (c + 1) % nc, 531 (c + 1) % nc, c, ccw) 532 533 if end_cap: 534 self.addFace([(ns - 1) * nc + x for x in range(0, nc)], ccw) 535 536# Triangle meshes 537 538 # Helper for numerous nodes with a Coordinate subnode holding vertices 539 # That all triangle meshes and IndexedFaceSet 540 # num_faces can be a function, in case the face count is a function of vertex count 541 def startCoordMesh(self, node, num_faces): 542 ccw = readBoolean(node, "ccw", True) 543 self.readVertices(node) # This will allocate and fill the vertex array 544 if hasattr(num_faces, "__call__"): 545 num_faces = num_faces(self.getVertexCount()) 546 self.reserveFaceCount(num_faces) 547 548 return ccw 549 550 551 def processGeometryIndexedTriangleSet(self, node): 552 index = readIntArray(node, "index", []) 553 num_faces = len(index) // 3 554 ccw = int(self.startCoordMesh(node, num_faces)) 555 556 for i in range(0, num_faces*3, 3): 557 self.addTri(index[i + 1 - ccw], index[i + ccw], index[i+2]) 558 559 def processGeometryIndexedTriangleStripSet(self, node): 560 strips = readIndex(node, "index") 561 ccw = int(self.startCoordMesh(node, sum([len(strip) - 2 for strip in strips]))) 562 563 for strip in strips: 564 sccw = ccw # Running CCW value, reset for each strip 565 for i in range(len(strip) - 2): 566 self.addTri(strip[i + 1 - sccw], strip[i + sccw], strip[i+2]) 567 sccw = 1 - sccw 568 569 def processGeometryIndexedTriangleFanSet(self, node): 570 fans = readIndex(node, "index") 571 ccw = int(self.startCoordMesh(node, sum([len(fan) - 2 for fan in fans]))) 572 573 for fan in fans: 574 for i in range(1, len(fan) - 1): 575 self.addTri(fan[0], fan[i + 1 - ccw], fan[i + ccw]) 576 577 def processGeometryTriangleSet(self, node): 578 ccw = int(self.startCoordMesh(node, lambda num_vert: num_vert // 3)) 579 for i in range(0, self.getVertexCount(), 3): 580 self.addTri(i + 1 - ccw, i + ccw, i+2) 581 582 def processGeometryTriangleStripSet(self, node): 583 strips = readIntArray(node, "stripCount", []) 584 ccw = int(self.startCoordMesh(node, sum([n-2 for n in strips]))) 585 586 vb = 0 587 for n in strips: 588 sccw = ccw 589 for i in range(n-2): 590 self.addTri(vb + i + 1 - sccw, vb + i + sccw, vb + i + 2) 591 sccw = 1 - sccw 592 vb += n 593 594 def processGeometryTriangleFanSet(self, node): 595 fans = readIntArray(node, "fanCount", []) 596 ccw = int(self.startCoordMesh(node, sum([n-2 for n in fans]))) 597 598 vb = 0 599 for n in fans: 600 for i in range(1, n-1): 601 self.addTri(vb, vb + i + 1 - ccw, vb + i + ccw) 602 vb += n 603 604 # Quad geometries from the CAD module, might be relevant for printing 605 606 def processGeometryQuadSet(self, node): 607 ccw = self.startCoordMesh(node, lambda num_vert: 2*(num_vert // 4)) 608 for i in range(0, self.getVertexCount(), 4): 609 self.addQuadFlip(i, i+1, i+2, i+3, ccw) 610 611 def processGeometryIndexedQuadSet(self, node): 612 index = readIntArray(node, "index", []) 613 num_quads = len(index) // 4 614 ccw = self.startCoordMesh(node, num_quads*2) 615 616 for i in range(0, num_quads*4, 4): 617 self.addQuadFlip(index[i], index[i+1], index[i+2], index[i+3], ccw) 618 619 # 2D polygon geometries 620 # Won't work for now, since Cura expects every mesh to have a nontrivial convex hull 621 # The only way around that is merging meshes. 622 623 def processGeometryDisk2D(self, node): 624 innerRadius = readFloat(node, "innerRadius", 0) 625 outerRadius = readFloat(node, "outerRadius", 1) 626 n = readInt(node, "subdivision", DEFAULT_SUBDIV) 627 628 angle = 2 * pi / n 629 630 self.reserveFaceAndVertexCount(n*4 if innerRadius else n-2, n*2 if innerRadius else n) 631 632 for i in range(n): 633 s = sin(angle * i) 634 c = cos(angle * i) 635 self.addVertex(outerRadius*c, outerRadius*s, 0) 636 if innerRadius: 637 self.addVertex(innerRadius*c, innerRadius*s, 0) 638 ni = (i+1) % n 639 self.addQuad(2*i, 2*ni, 2*ni+1, 2*i+1) 640 641 if not innerRadius: 642 for i in range(2, n): 643 self.addTri(0, i-1, i) 644 645 def processGeometryRectangle2D(self, node): 646 (x, y) = readFloatArray(node, "size", (2, 2)) 647 self.reserveFaceAndVertexCount(2, 4) 648 self.addVertex(-x/2, -y/2, 0) 649 self.addVertex(x/2, -y/2, 0) 650 self.addVertex(x/2, y/2, 0) 651 self.addVertex(-x/2, y/2, 0) 652 self.addQuad(0, 1, 2, 3) 653 654 def processGeometryTriangleSet2D(self, node): 655 verts = readFloatArray(node, "vertices", ()) 656 num_faces = len(verts) // 6 657 verts = [(verts[i], verts[i+1], 0) for i in range(0, 6 * num_faces, 2)] 658 self.reserveFaceAndVertexCount(num_faces, num_faces * 3) 659 for vert in verts: 660 self.addVertex(*vert) 661 662 # The front face is on the +Z side, so CCW is a variable 663 for i in range(0, num_faces*3, 3): 664 a = Vector(*verts[i+2]) - Vector(*verts[i]) 665 b = Vector(*verts[i+1]) - Vector(*verts[i]) 666 self.addTriFlip(i, i+1, i+2, a.x*b.y > a.y*b.x) 667 668 # General purpose polygon mesh 669 670 def processGeometryIndexedFaceSet(self, node): 671 faces = readIndex(node, "coordIndex") 672 ccw = self.startCoordMesh(node, sum([len(face) - 2 for face in faces])) 673 674 for face in faces: 675 if len(face) == 3: 676 self.addTriFlip(face[0], face[1], face[2], ccw) 677 elif len(face) > 3: 678 self.addFace(face, ccw) 679 680 geometry_importers = { 681 "IndexedFaceSet": processGeometryIndexedFaceSet, 682 "IndexedTriangleSet": processGeometryIndexedTriangleSet, 683 "IndexedTriangleStripSet": processGeometryIndexedTriangleStripSet, 684 "IndexedTriangleFanSet": processGeometryIndexedTriangleFanSet, 685 "TriangleSet": processGeometryTriangleSet, 686 "TriangleStripSet": processGeometryTriangleStripSet, 687 "TriangleFanSet": processGeometryTriangleFanSet, 688 "QuadSet": processGeometryQuadSet, 689 "IndexedQuadSet": processGeometryIndexedQuadSet, 690 "TriangleSet2D": processGeometryTriangleSet2D, 691 "Rectangle2D": processGeometryRectangle2D, 692 "Disk2D": processGeometryDisk2D, 693 "ElevationGrid": processGeometryElevationGrid, 694 "Extrusion": processGeometryExtrusion, 695 "Sphere": processGeometrySphere, 696 "Box": processGeometryBox, 697 "Cylinder": processGeometryCylinder, 698 "Cone": processGeometryCone 699 } 700 701 # Parses the Coordinate.@point field, fills the verts array. 702 def readVertices(self, node): 703 for c in node: 704 if c.tag == "Coordinate": 705 c = self.resolveDefUse(c) 706 if not c is None: 707 pt = c.attrib.get("point") 708 if pt: 709 # allow the list of float values in 'point' attribute to 710 # be separated by commas or whitespace as per spec of 711 # XML encoding of X3D 712 # Ref ISO/IEC 19776-1:2015 : Section 5.1.2 713 co = [float(x) for vec in pt.split(',') for x in vec.split()] 714 num_verts = len(co) // 3 715 self.verts = numpy.empty((4, num_verts), dtype=numpy.float32) 716 self.verts[3,:] = numpy.ones((num_verts), dtype=numpy.float32) 717 # Group by three 718 for i in range(num_verts): 719 self.verts[:3,i] = co[3*i:3*i+3] 720 721 # Mesh builder helpers 722 723 def reserveFaceAndVertexCount(self, num_faces, num_verts): 724 # Unlike the Cura MeshBuilder, we use 4-vectors stored as columns for easier transform 725 self.verts = numpy.zeros((4, num_verts), dtype=numpy.float32) 726 self.verts[3,:] = numpy.ones((num_verts), dtype=numpy.float32) 727 self.num_verts = 0 728 self.reserveFaceCount(num_faces) 729 730 def reserveFaceCount(self, num_faces): 731 self.faces = numpy.zeros((num_faces, 3), dtype=numpy.int32) 732 self.num_faces = 0 733 734 def getVertexCount(self): 735 return self.verts.shape[1] 736 737 def addVertex(self, x, y, z): 738 self.verts[0, self.num_verts] = x 739 self.verts[1, self.num_verts] = y 740 self.verts[2, self.num_verts] = z 741 self.num_verts += 1 742 743 # Indices are 0-based for this shape, but they won't be zero-based in the merged mesh 744 def addTri(self, a, b, c): 745 self.faces[self.num_faces, 0] = self.index_base + a 746 self.faces[self.num_faces, 1] = self.index_base + b 747 self.faces[self.num_faces, 2] = self.index_base + c 748 self.num_faces += 1 749 750 def addTriFlip(self, a, b, c, ccw): 751 if ccw: 752 self.addTri(a, b, c) 753 else: 754 self.addTri(b, a, c) 755 756 # Needs to be convex, but not necessaily planar 757 # Assumed ccw, cut along the ac diagonal 758 def addQuad(self, a, b, c, d): 759 self.addTri(a, b, c) 760 self.addTri(c, d, a) 761 762 def addQuadFlip(self, a, b, c, d, ccw): 763 if ccw: 764 self.addTri(a, b, c) 765 self.addTri(c, d, a) 766 else: 767 self.addTri(a, c, b) 768 self.addTri(c, a, d) 769 770 771 # Arbitrary polygon triangulation. 772 # Doesn't assume convexity and doesn't check the "convex" flag in the file. 773 # Works by the "cutting of ears" algorithm: 774 # - Find an outer vertex with the smallest angle and no vertices inside its adjacent triangle 775 # - Remove the triangle at that vertex 776 # - Repeat until done 777 # Vertex coordinates are supposed to be already set 778 def addFace(self, indices, ccw): 779 # Resolve indices to coordinates for faster math 780 face = [Vector(data=self.verts[0:3, i]) for i in indices] 781 782 # Need a normal to the plane so that we can know which vertices form inner angles 783 normal = findOuterNormal(face) 784 785 if not normal: # Couldn't find an outer edge, non-planar polygon maybe? 786 return 787 788 # Find the vertex with the smallest inner angle and no points inside, cut off. Repeat until done 789 n = len(face) 790 vi = [i for i in range(n)] # We'll be using this to kick vertices from the face 791 while n > 3: 792 max_cos = EPSILON # We don't want to check anything on Pi angles 793 i_min = 0 # max cos corresponds to min angle 794 for i in range(n): 795 inext = (i + 1) % n 796 iprev = (i + n - 1) % n 797 v = face[vi[i]] 798 next = face[vi[inext]] - v 799 prev = face[vi[iprev]] - v 800 nextXprev = next.cross(prev) 801 if nextXprev.dot(normal) > EPSILON: # If it's an inner angle 802 cos = next.dot(prev) / (next.length() * prev.length()) 803 if cos > max_cos: 804 # Check if there are vertices inside the triangle 805 no_points_inside = True 806 for j in range(n): 807 if j != i and j != iprev and j != inext: 808 vx = face[vi[j]] - v 809 if pointInsideTriangle(vx, next, prev, nextXprev): 810 no_points_inside = False 811 break 812 813 if no_points_inside: 814 max_cos = cos 815 i_min = i 816 817 self.addTriFlip(indices[vi[(i_min + n - 1) % n]], indices[vi[i_min]], indices[vi[(i_min + 1) % n]], ccw) 818 vi.pop(i_min) 819 n -= 1 820 self.addTriFlip(indices[vi[0]], indices[vi[1]], indices[vi[2]], ccw) 821 822 823# ------------------------------------------------------------ 824# X3D field parsers 825# ------------------------------------------------------------ 826def readFloatArray(node, attr, default): 827 s = node.attrib.get(attr) 828 if not s: 829 return default 830 return [float(x) for x in s.split()] 831 832def readIntArray(node, attr, default): 833 s = node.attrib.get(attr) 834 if not s: 835 return default 836 return [int(x, 0) for x in s.split()] 837 838def readFloat(node, attr, default): 839 s = node.attrib.get(attr) 840 if not s: 841 return default 842 return float(s) 843 844def readInt(node, attr, default): 845 s = node.attrib.get(attr) 846 if not s: 847 return default 848 return int(s, 0) 849 850def readBoolean(node, attr, default): 851 s = node.attrib.get(attr) 852 if not s: 853 return default 854 return s.lower() == "true" 855 856def readVector(node, attr, default): 857 v = readFloatArray(node, attr, default) 858 return Vector(v[0], v[1], v[2]) 859 860def readRotation(node, attr, default): 861 v = readFloatArray(node, attr, default) 862 return (v[3], Vector(v[0], v[1], v[2])) 863 864# Returns the -1-separated runs 865def readIndex(node, attr): 866 v = readIntArray(node, attr, []) 867 chunks = [] 868 chunk = [] 869 for i in range(len(v)): 870 if v[i] == -1: 871 if chunk: 872 chunks.append(chunk) 873 chunk = [] 874 else: 875 chunk.append(v[i]) 876 if chunk: 877 chunks.append(chunk) 878 return chunks 879 880# Given a face as a sequence of vectors, returns a normal to the polygon place that forms a right triple 881# with a vector along the polygon sequence and a vector backwards 882def findOuterNormal(face): 883 n = len(face) 884 for i in range(n): 885 for j in range(i+1, n): 886 edge = face[j] - face[i] 887 if edge.length() > EPSILON: 888 edge = edge.normalized() 889 prev_rejection = Vector() 890 is_outer = True 891 for k in range(n): 892 if k != i and k != j: 893 pt = face[k] - face[i] 894 pte = pt.dot(edge) 895 rejection = pt - edge*pte 896 if rejection.dot(prev_rejection) < -EPSILON: # points on both sides of the edge - not an outer one 897 is_outer = False 898 break 899 elif rejection.length() > prev_rejection.length(): # Pick a greater rejection for numeric stability 900 prev_rejection = rejection 901 902 if is_outer: # Found an outer edge, prev_rejection is the rejection inside the face. Generate a normal. 903 return edge.cross(prev_rejection) 904 905 return False 906 907 908# Given two *collinear* vectors a and b, returns the coefficient that takes b to a. 909# No error handling. 910# For stability, taking the ration between the biggest coordinates would be better... 911def ratio(a, b): 912 if b.x > EPSILON or b.x < -EPSILON: 913 return a.x / b.x 914 elif b.y > EPSILON or b.y < -EPSILON: 915 return a.y / b.y 916 else: 917 return a.z / b.z 918 919 920def pointInsideTriangle(vx, next, prev, nextXprev): 921 vxXprev = vx.cross(prev) 922 r = ratio(vxXprev, nextXprev) 923 if r < 0: 924 return False 925 vxXnext = vx.cross(next) 926 s = -ratio(vxXnext, nextXprev) 927 return s > 0 and (s + r) < 1 928 929