From 4dd8c1e4e1eaa2c118d35da45b7a7019c5d8ac5d Mon Sep 17 00:00:00 2001 From: Spiros Tsalikis Date: Sat, 23 Nov 2024 20:41:48 -0500 Subject: [PATCH 1/5] vtkSurfaceNets3D: Improve variable names and use std containers --- Filters/Core/vtkSurfaceNets3D.cxx | 1070 ++++++++++++++--------------- 1 file changed, 524 insertions(+), 546 deletions(-) diff --git a/Filters/Core/vtkSurfaceNets3D.cxx b/Filters/Core/vtkSurfaceNets3D.cxx index 52f46eecc151..49620c3ac68a 100644 --- a/Filters/Core/vtkSurfaceNets3D.cxx +++ b/Filters/Core/vtkSurfaceNets3D.cxx @@ -35,13 +35,13 @@ vtkStandardNewMacro(vtkSurfaceNets3D); // vtkConstrainedSmoothingFilter performs the smoothing. // // A templated surface nets extraction algorithm implementation follows. It -// uses a edge-by-edge parallel algorithm (aka flying edges approach) for +// uses an edge-by-edge parallel algorithm (aka flying edges approach) for // performance. There are four passes to surface extraction algorithm: 1) // classify x-edges. 2) classify y-z-edges, 3) perform a prefix sum to -// determine where to write / allocate output data, and 4) a output +// determine where to write / allocate output data, and 4) an output // generation pass (i.e., generate points, polygons, and optional scalar // data). An optional fifth step smooths this output mesh to improve mesh -// quailty. +// quality. // // Some terminology: Eight voxel points (which in VTK is point-associated data) are // combined to create regular hexahedron (which in VTK are voxel @@ -85,9 +85,9 @@ vtkStandardNewMacro(vtkSurfaceNets3D); // intersected edge, then the two voxel faces using that edge will also be // "intersected" by a smoothing stencil connection. Therefore, an edge case // number can be converted into a face case number. Note because the edge -// case number exceeds what is representable by a 8-bit unsigned char, the +// case number exceeds what is representable by an 8-bit unsigned char, the // edge case number is represented by a 16-bit unsigned short (but is -// converted to a 8-bit unsigned char, face case number as needed). (One +// converted to an 8-bit unsigned char, face case number as needed). (One // interesting note about edge cases: although as implemented here the 12 // voxel edges can be set independently, this is geometrically impossible // since if one edge is intersected, than at least two others must be as @@ -97,7 +97,7 @@ vtkStandardNewMacro(vtkSurfaceNets3D); namespace { // anonymous -// The SurfaceNets struct below implements the core of the surface nets +// The SurfaceNets struct below implements the core of the surface nets' // algorithm. It uses a flying edges approach to parallel process data // edge-by-edge, which provides edge-based parallel tasking, reduces the // number of voxel lookups and eliminates costly coincident point merging. @@ -108,7 +108,7 @@ namespace // the information gathered is whether the x-edge, y-edge, and/or z-edge // requires "intersection", whether a point needs to be inserted into the // center of the voxel, and whether the voxel origin point/triad origin is -// inside of any labeled region, or outside. In Pass#3, a prefix sum is used +// inside any labeled region, or outside. In Pass#3, a prefix sum is used // to characterize the output, and allocate the appropriate output // arrays. Finally, in Pass#4, the final output points, surface, and // smoothing stencils are produced. Following the surface extraction, an @@ -118,7 +118,7 @@ namespace // since smoothing generally causes the quads to become non-planar. // // A key concept of this implementation is EdgeMetaData (often referred to as -// eMD[5]). The edge meta data maintains information about each volume x-edge +// eMD[5]). The edge metadata maintains information about each volume x-edge // (i.e., row) which is necessary for threading the implementation. The // information maintained is: eMD[0]- the number points produced along the // x-row; eMD[1]- the number of quad primitives produced from this row; @@ -133,68 +133,77 @@ namespace // padding) which keeps track of information needed to thread the algorithm // across the volume x-edges. // -// Another way to look at this: the edge meta data characterizes each row of -// voxel triads. eMD keeps track of the the number of points, quads, and stencil +// Another way to look at this: the edge metadata characterizes each row of +// voxel triads. eMD keeps track of the number of points, quads, and stencil // edges generated by a row of voxel triads. It also maintains a clipped region // [xMin_i,xMax_i) or [xL,xR) along the edge of voxel triads in which primitives -// may be generated (i.e., tracks computational trimming). The edge meta data +// may be generated (i.e., tracks computational trimming). The edge metadata // provides bookkeeping necessary to support threaded computing. // Some structs and typedefs to clarify code. using TriadType = unsigned char; using EdgeCaseType = unsigned short; using FaceCaseType = unsigned char; -using FaceCaseType = unsigned char; +using TrimmedEdgesCaseType = unsigned char; + +struct EdgeMetaDataType +{ + vtkIdType NumPoints = 0; // number of points produced along this row + vtkIdType NumQuads = 0; // number of quad primitives produced along this row + vtkIdType NumStencilEdges = 0; // number of stencil edges + vtkIdType XMin = 0; // minimum index of first intersection along this row + vtkIdType XMax = -1; // maximum index of intersection along this row +}; + +// const values to access the correct dimension of the data +enum Dim : std::uint8_t +{ + X = 0, + Y = 1, + Z = 2 +}; template struct SurfaceNets { // The triad classification carries information on five different bits. - // The first bit Bit 1 indicates whether the origin of the triad is inside - // or outside of *any* labeled region. Bit 2 indicates whether the x-edge - // needs intersection (i.e., a surface net passes through it); Bit 3 - // whether the y-edge needs intersection; and Bit 4 whether the z-edge - // needs intersection. (Triad edges require intersection when the two end + // Bit 1 indicates whether the origin of the triad is inside or outside + // *any* labeled region. Bit 2 indicates whether the x-edge needs + // intersection (i.e., a surface net passes through it); Bit 3 whether + // the y-edge needs intersection; and Bit 4 whether the z-edge needs + // intersection. (Triad edges require intersection when the two end // point values are not equal to one another, and at least one of the end // point values is "Inside" a labeled region.) Finally, the fifth bit is // used to indicate whether a point will be generated in the voxel cube/cell // associated with a triad. This fifth bit (ProducePoint) is used to // simplify and speed up code. - enum TriadClassification + enum TriadClassification : std::uint8_t { - Outside = 0, // triad origin point is outside of any labeled region - Inside = 1, // triad origin inside of some labeled region + Outside = 0, // triad origin point is outside any labeled region + Inside = 1, // triad origin inside some labeled region XIntersection = 2, // triad x-axis requires intersection YIntersection = 4, // triad y-axis requires intersection ZIntersection = 8, // triad z-axis requires intersection ProducePoint = 16 // the voxel associated with this triad will produce a point }; - enum SurfaceFaceClassification - { - SolidFace = 0, - SurfaceFace = 1, - JunctionFace = 2, - ComplexFace = 3 - }; - // Given a pointer to a voxel's triad, first determine the seven triad cases // (from the points defining a voxel cell: (x,y,z); ([x+1],y,z); (x,[y+1],z); // ([x+1],[y+1],z); (x,y,[z+1]); ([x+1],y,[z+1]); (x,[y+1],[z+1]), and then // compute the edge case number for this voxel cell. Note that a resulting // value of zero means that the voxel cell is not intersected (i.e., no edge is - // intersected). This method assumes that the tPtr is not on the boundary + // intersected). This method assumes that the triadPtr is not on the boundary // of the padded volume. - EdgeCaseType GetEdgeCase(TriadType* tPtr) + EdgeCaseType GetEdgeCase(TriadType* triadPtr) { TriadType triads[7]; - triads[0] = tPtr[0]; - triads[1] = tPtr[1]; - triads[2] = tPtr[this->TriadDims[0]]; - triads[3] = tPtr[1 + this->TriadDims[0]]; - triads[4] = tPtr[this->TriadSliceOffset]; - triads[5] = tPtr[1 + this->TriadSliceOffset]; - triads[6] = tPtr[this->TriadDims[0] + this->TriadSliceOffset]; + triads[0] = triadPtr[0]; + triads[1] = triadPtr[1]; + triads[2] = triadPtr[this->TriadDims[X]]; + triads[3] = triadPtr[this->TriadDims[X] + 1]; + triads[4] = triadPtr[this->TriadSliceOffset]; + triads[5] = triadPtr[this->TriadSliceOffset + 1]; + triads[6] = triadPtr[this->TriadSliceOffset + this->TriadDims[X]]; // Process the selected twelve edges from the seven triads to produce an // edge case number. The triad numbering is the same as a vtkVoxel point @@ -202,84 +211,83 @@ struct SurfaceNets // numbering: first the four voxel x-edges, then the four y-edges, then the // four voxel z-edges. // x-edges - EdgeCaseType eCase = (triads[0] & 2) >> 1; - eCase |= (triads[2] & 2); - eCase |= (triads[4] & 2) << 1; - eCase |= (triads[6] & 2) << 2; + EdgeCaseType edgeCase = (triads[0] & TriadClassification::XIntersection) >> 1; + edgeCase |= (triads[2] & TriadClassification::XIntersection); + edgeCase |= (triads[4] & TriadClassification::XIntersection) << 1; + edgeCase |= (triads[6] & TriadClassification::XIntersection) << 2; // y-edges - eCase |= (triads[0] & 4) << 2; - eCase |= (triads[1] & 4) << 3; - eCase |= (triads[4] & 4) << 4; - eCase |= (triads[5] & 4) << 5; + edgeCase |= (triads[0] & TriadClassification::YIntersection) << 2; + edgeCase |= (triads[1] & TriadClassification::YIntersection) << 3; + edgeCase |= (triads[4] & TriadClassification::YIntersection) << 4; + edgeCase |= (triads[5] & TriadClassification::YIntersection) << 5; // z-edges - eCase |= (triads[0] & 8) << 5; - eCase |= (triads[1] & 8) << 6; - eCase |= (triads[2] & 8) << 7; - eCase |= (triads[3] & 8) << 8; + edgeCase |= (triads[0] & TriadClassification::ZIntersection) << 5; + edgeCase |= (triads[1] & TriadClassification::ZIntersection) << 6; + edgeCase |= (triads[2] & TriadClassification::ZIntersection) << 7; + edgeCase |= (triads[3] & TriadClassification::ZIntersection) << 8; - return eCase; + return edgeCase; } // GetEdgeCase // Given a voxel cell edge case, convert it to a voxel face case. While // this could be done through a table, the size of the table is large // enough that a procedural approach simplifies the code. Basically, each // intersected voxel cell edge will activate two voxel faces. - FaceCaseType GetFaceCase(EdgeCaseType edgeCase) + static FaceCaseType GetFaceCase(EdgeCaseType edgeCase) { - FaceCaseType fCase = 0; - + FaceCaseType faceCase = 0; // Process each of the voxel's twelve edges. If edge is set, then set // the two faces using the edge. - if (edgeCase & 1) // edge 0 + if (edgeCase & 1) // edge 0, faces 2 & 4 { - fCase |= 20; // faces 2 & 4 + faceCase |= 20; } - if (edgeCase & 2) // edge 1 + if (edgeCase & 2) // edge 1, faces 3 & 4 { - fCase |= 24; // faces 3 & 4 + faceCase |= 24; } - if (edgeCase & 4) // edge 2 + if (edgeCase & 4) // edge 2, faces 2 & 5 { - fCase |= 36; // faces 2 & 5 + faceCase |= 36; } - if (edgeCase & 8) // edge 3 + if (edgeCase & 8) // edge 3, faces 3 & 5 { - fCase |= 40; // faces 3 & 5 + faceCase |= 40; } - if (edgeCase & 16) // edge 4 + if (edgeCase & 16) // edge 4, faces 0 & 4 { - fCase |= 17; // faces 0 & 4 + faceCase |= 17; } - if (edgeCase & 32) // edge 5 + if (edgeCase & 32) // edge 5, faces 1 & 4 { - fCase |= 18; // faces 1 & 4 + faceCase |= 18; } - if (edgeCase & 64) // edge 6 + if (edgeCase & 64) // edge 6, faces 0 & 5 { - fCase |= 33; // faces 0 & 5 + faceCase |= 33; } - if (edgeCase & 128) // edge 7 + if (edgeCase & 128) // edge 7, faces 1 & 5 { - fCase |= 34; // faces 1 & 5 + faceCase |= 34; } - if (edgeCase & 256) // edge 8 + if (edgeCase & 256) // edge 8, faces 0 & 2 { - fCase |= 5; // faces 0 & 2 + faceCase |= 5; } - if (edgeCase & 512) // edge 9 + if (edgeCase & 512) // edge 9, faces 1 & 2 { - fCase |= 6; // faces 1 & 2 + faceCase |= 6; } - if (edgeCase & 1024) // edge 10 + if (edgeCase & 1024) // edge 10, faces 0 & 3 { - fCase |= 9; // faces 0 & 3 + faceCase |= 9; } - if (edgeCase & 2048) // edge 11 + if (edgeCase & 2048) // edge 11, faces 1 & 3 { - fCase |= 10; // faces 1 & 3 + faceCase |= 10; } - return fCase; + return faceCase; } // GetFaceCase // Perform analysis of voxel edge case: count the number of edge intersections on each @@ -290,7 +298,7 @@ struct SurfaceNets // Process each of the voxel's twelve edges. If edge is set, then increment // the two face counts using the edge. - if (edgeCase & 1) // edge 0, affects faces 2 & 4 + if (edgeCase & 1) // edge 0, faces 2 & 4 { faceCounts[2]++; faceCounts[4]++; @@ -351,22 +359,36 @@ struct SurfaceNets faceCounts[3]++; } - return (*std::max_element(faceCounts, faceCounts + 6)); + return *std::max_element(faceCounts, faceCounts + 6); } // CountFaceIntersections // Obtain information indicating whether quad polygons are to be generated // from the triad specified. A triad may produce up to three quad polygons - // corresponding to the lower left corner of a voxel. One is a x-y quad; a + // corresponding to the lower left corner of a voxel. One is an x-y quad; an // x-z quad, and a y-z quad. - static bool GenerateXYQuad(TriadType triad) { return ((triad & SurfaceNets::ZIntersection) > 0); } - static bool GenerateXZQuad(TriadType triad) { return ((triad & SurfaceNets::YIntersection) > 0); } - static bool GenerateYZQuad(TriadType triad) { return ((triad & SurfaceNets::XIntersection) > 0); } - static bool ProducesQuad(TriadType triad) { return ((triad & 14) > 0); } - static unsigned char GetNumberOfQuads(TriadType triad) + static bool GenerateXYQuad(TriadType triad) + { + return (triad & TriadClassification::ZIntersection) > 0; + } + static bool GenerateXZQuad(TriadType triad) { - unsigned char numQuads = (SurfaceNets::GenerateXYQuad(triad) ? 1 : 0); - numQuads += (SurfaceNets::GenerateXZQuad(triad) ? 1 : 0); - numQuads += (SurfaceNets::GenerateYZQuad(triad) ? 1 : 0); + return (triad & TriadClassification::YIntersection) > 0; + } + static bool GenerateYZQuad(TriadType triad) + { + return (triad & TriadClassification::XIntersection) > 0; + } + static bool ProducesQuad(TriadType triad) + { + static constexpr TriadType triadMask = TriadClassification::XIntersection | + TriadClassification::YIntersection | TriadClassification::ZIntersection; + return (triad & triadMask) > 0; + } + static uint8_t GetNumberOfQuads(TriadType triad) + { + uint8_t numQuads = SurfaceNets::GenerateXYQuad(triad); + numQuads += SurfaceNets::GenerateXZQuad(triad); + numQuads += SurfaceNets::GenerateYZQuad(triad); return numQuads; } @@ -377,32 +399,36 @@ struct SurfaceNets // which of the six edge are to be generated. The table of stencil cases was // generated programmatically (see GenerateStencils). static const unsigned char StencilFaceCases[64][7]; - static unsigned char GetNumberOfStencilFaceEdges(FaceCaseType caseNum) + + static unsigned char GetNumberOfStencilFaceEdges(FaceCaseType faceCase) { - return SurfaceNets::StencilFaceCases[caseNum][0]; + return SurfaceNets::StencilFaceCases[faceCase][0]; } - static const unsigned char* GetStencilFaceEdges(FaceCaseType caseNum) + static const unsigned char* GetStencilFaceEdges(FaceCaseType faceCase) { - return SurfaceNets::StencilFaceCases[caseNum]; + return SurfaceNets::StencilFaceCases[faceCase]; } void GenerateFaceStencils(unsigned char stencils[][7]); // This smoothing stencil table is indexed by the voxel *edge* case. It // indexes into the face-case-based smoothing stenciles. This table is // constructed programmatically. - unsigned char GetNumberOfStencilEdges(EdgeCaseType caseNum) + unsigned char GetNumberOfStencilEdges(EdgeCaseType edgeCase) { - return SurfaceNets::StencilFaceCases[this->StencilTable[caseNum]][0]; + return SurfaceNets::StencilFaceCases[this->StencilTable[edgeCase]][0]; } - const unsigned char* GetStencilEdges(EdgeCaseType caseNum) + const unsigned char* GetStencilEdges(EdgeCaseType edgeCase) { - return SurfaceNets::StencilFaceCases[this->StencilTable[caseNum]]; + return SurfaceNets::StencilFaceCases[this->StencilTable[edgeCase]]; } void GenerateEdgeStencils(int optLevel = 0); // Return whether a triad, and its associated voxel cell, requires the // generation of a point. - bool ProducesPoint(TriadType triad) { return (triad & SurfaceNets::ProducePoint) > 0; } + static bool ProducesPoint(TriadType triad) + { + return (triad & TriadClassification::ProducePoint) > 0; + } // Input and output data. T* Scalars; // input image scalars @@ -419,24 +445,15 @@ struct SurfaceNets // Internal variables used by the various algorithm methods. Interfaces VTK // image data in an efficient form more convenient to the algorithm. vtkIdType Dims[3]; - int Min0; - int Max0; - int Inc0; - int Min1; - int Max1; - int Inc1; - int Min2; - int Max2; - int Inc2; - - // Algorithm-derived data for bookkeeping data locations - // when parallel computing. - TriadType* Triads; + int Min[3]; + int Max[3]; + int Inc[3]; + + // Algorithm-derived data for bookkeeping data locations when parallel computing. + std::vector Triads; vtkIdType TriadDims[3]; vtkIdType TriadSliceOffset; - static constexpr int EdgeMetaDataSize = 5; - vtkIdType NumberOfEdges; - vtkIdType* EdgeMetaData; + std::vector EdgeMetaData; // The stencil table used to obtain smoothing stencils from the voxel *edge // case*. This table indexes into the StencilFaceCases[64][7] using the voxel @@ -456,130 +473,111 @@ struct SurfaceNets , LabelValues(nullptr) , BackgroundLabel(0) , Dims{ 0, 0, 0 } - , Min0(0) - , Max0(0) - , Inc0(0) - , Min1(0) - , Max1(0) - , Inc1(0) - , Min2(0) - , Max2(0) - , Inc2(0) - , Triads(nullptr) + , Min{ 0, 0, 0 } + , Max{ 0, 0, 0 } + , Inc{ 0, 0, 0 } , TriadDims{ 0, 0, 0 } , TriadSliceOffset(0) - , NumberOfEdges(0) - , EdgeMetaData(nullptr) { this->GenerateEdgeStencils(1); } - // Place holder for now in case fancy bit fiddling is needed later. - void SetTriadClassification(TriadType* tPtr, TriadType tCase) { *tPtr = tCase; } - // Classify the triad y-edges. Use the triad cases at both ends of the y-edge // first; if necessary, access the voxel values. The indices i and row are // expressed in the triad coordinates. - unsigned char ClassifyYEdge( - T* inPtr, vtkIdType i, TriadType case0, vtkIdType row, TriadType case1) + TriadClassification ClassifyYEdge( + T* inPtr, vtkIdType i, TriadType triad0, vtkIdType row, TriadType triad1) { // If on padded boundary, edge is never intersected. - if (row >= (this->TriadDims[1] - 2)) + if (row >= this->Dims[Y]) { - return 0; + return TriadClassification::Outside; } - unsigned char inout0 = (case0 & 0x1); - unsigned char inout1 = (case1 & 0x1); + const TriadType inout0 = triad0 & TriadClassification::Inside; + const TriadType inout1 = triad1 & TriadClassification::Inside; if (inout0 == inout1) { - if (inout0 == Outside) - { // both triad origins are outside - return 0; + if (inout0 == TriadClassification::Outside) // both triad origins are outside + { + return TriadClassification::Outside; } - else - { // both triad origins are inside, need to check regions - vtkIdType idx = i - 1; // shift into volume (i.e, no padding) - T s0 = inPtr[idx]; - T s1 = inPtr[idx + this->Inc1]; - return (s0 == s1 ? 0 : YIntersection); + else // both triad origins are inside, need to check regions + { + const vtkIdType idx = i - 1; // shift into volume (i.e, no padding) + const T& s0 = inPtr[idx]; + const T& s1 = inPtr[idx + this->Inc[Y]]; + return s0 == s1 ? TriadClassification::Outside : TriadClassification::YIntersection; } } - else - { // one triad origin point is inside, one outside, so y-edge-intersection - return YIntersection; + else // one triad origin point is inside, one outside, so y-edge-intersection + { + return TriadClassification::YIntersection; } } // ClassifyYEdge // Classify the triad z-edges. Use the triad cases at both ends of the z-edge // first; if necessary, access the voxel values. The indices i and slice are // expressed in the triad coordinates. - unsigned char ClassifyZEdge( - T* inPtr, vtkIdType i, TriadType case0, vtkIdType slice, TriadType case1) + TriadClassification ClassifyZEdge( + T* inPtr, vtkIdType i, TriadType triad0, vtkIdType slice, TriadType triad1) { // If on padded boundary, edge is never intersected. - if (slice >= (this->TriadDims[2] - 2)) + if (slice >= this->Dims[Z]) { - return 0; + return TriadClassification::Outside; } - unsigned char inout0 = (case0 & 0x1); - unsigned char inout1 = (case1 & 0x1); + const TriadType inout0 = triad0 & TriadClassification::Inside; + const TriadType inout1 = triad1 & TriadClassification::Inside; if (inout0 == inout1) { - if (inout0 == Outside) - { // both triad origins are outside - return 0; + if (inout0 == TriadClassification::Outside) // both triad origins are outside + { + return TriadClassification::Outside; } - else - { // both triad origins are inside, need to check regions - vtkIdType idx = i - 1; // shift into volume (i.e., no padding) - T s0 = inPtr[idx]; - T s1 = inPtr[idx + this->Inc2]; - return (s0 == s1 ? 0 : ZIntersection); + else // both triad origins are inside, need to check regions + { + const vtkIdType idx = i - 1; // shift into volume (i.e., no padding) + const T& s0 = inPtr[idx]; + const T& s1 = inPtr[idx + this->Inc[Z]]; + return s0 == s1 ? TriadClassification::Outside : TriadClassification::ZIntersection; } } - else - { // one triad origin point is inside, one outside, so z-edge-intersection - return ZIntersection; + else // one triad origin point is inside, one outside, so z-edge-intersection + { + return TriadClassification::ZIntersection; } } // ClassifyZEdge // Composite trimming information to determine which portion of the // volume x-edge (row,slice) to process. In particular, gather the 2x2 // trim edge metadata that forms a row of voxel cells. - vtkIdType* Get2x2EdgeTrim(vtkIdType row, vtkIdType slice, vtkIdType& xL, vtkIdType& xR) + EdgeMetaDataType* Get2x2EdgeTrim(vtkIdType row, vtkIdType slice, vtkIdType& xMin, vtkIdType& xMax) { - // Grab the meta data for the 2x2 bundle of rows. - vtkIdType* EMD = this->EdgeMetaData; - vtkIdType* dims = this->TriadDims; - vtkIdType size = this->EdgeMetaDataSize; - - // Gather the four edge metadata rows that form a column of voxel cells. - vtkIdType* ePtrs[4] = { nullptr, nullptr, nullptr, nullptr }; - ePtrs[0] = EMD + (slice * dims[1] + row) * size; // current edge row - ePtrs[1] = ePtrs[0] + size; // to the right of the current edge - ePtrs[2] = ePtrs[0] + dims[1] * size; // above the current edge - ePtrs[3] = ePtrs[2] + size; // above and to the right of the current edge + // Gather the metadata for the four (2x2) edge rows that form a column of voxel cells. + std::array eMDPtrs; + eMDPtrs[0] = this->EdgeMetaData.data() + (slice * this->TriadDims[Y] + row); // current edge row + eMDPtrs[1] = eMDPtrs[0] + 1; // to the right + eMDPtrs[2] = eMDPtrs[0] + this->TriadDims[Y]; // above + eMDPtrs[3] = eMDPtrs[2] + 1; // above and to the right // Determine the trim over the 2x2 bundle of metadata. - xL = this->TriadDims[0]; - xR = 0; - for (auto i = 0; i < 4; ++i) + xMin = this->TriadDims[X]; + xMax = 0; + for (auto& eMDPtr : eMDPtrs) { - vtkIdType* eMD = ePtrs[i]; - xL = (eMD[3] < xL ? eMD[3] : xL); - xR = (eMD[4] > xR ? eMD[4] : xR); + xMin = std::min(xMin, eMDPtr->XMin); + xMax = std::max(xMax, eMDPtr->XMax); } - return ePtrs[0]; + return eMDPtrs[0]; } // Get2x2EdgeTrim - using TrimmedEdgesCase = unsigned char; // Composite the trimming information to determine which portion of the // volume x-edge (row,slice) to process. Since processing occurs across 3x3 // bundles of edges, we need to composite the metadata from these nine // edges to determine trimming. Also get the 3x3 triads and 3x3 bundle of - // edge meta data. This function always return not nullptr ePtrs, and tPtrs for + // edge metadata. This function always return not nullptr ePtrs, and tPtrs for // the 4,5,7,8 indices. The index 0 is not nullptr if row != 0 && slice != 0. // The index 1 is not nullptr if slice != 0. The index 2 is not nullptr if // row != 0 && slice != 0. So there are basically 4 cases. @@ -587,76 +585,70 @@ struct SurfaceNets // if return value is 1: row != 0 && slice == 0 // if return value is 2: row == 0 && slice != 0 // if return value is 3: row != 0 && slice != 0 - TrimmedEdgesCase Get3x3EdgeTrim(vtkIdType row, vtkIdType slice, vtkIdType& xL, vtkIdType& xR, - vtkIdType* ePtrs[9], TriadType* tPtrs[9]) + TrimmedEdgesCaseType Get3x3EdgeTrim(vtkIdType row, vtkIdType slice, vtkIdType& xMin, + vtkIdType& xMax, std::array& eMDPtrs, + std::array& triadPtrs) { - // Grab the meta data for the 3x3 bundle of rows. Watch out for + // Grab the metadata for the 3x3 bundle of rows. Watch out for // bundles near the (-x,-y,-z) boundaries. (The (+x,+y,+z) boundaries // are always okay due to the nature of the padding, and iteration // over rows and slices). - vtkIdType* EMD = this->EdgeMetaData; - vtkIdType* dims = this->TriadDims; - vtkIdType size = this->EdgeMetaDataSize; - TriadType* triads = this->Triads; - vtkIdType sliceOffset = this->TriadSliceOffset; - TrimmedEdgesCase trimmedEdgesCase = (row != 0) + ((slice != 0) << 1); - - // Initialize the triads and edge meta data. This simplifies - // the code. - std::fill_n(ePtrs, 9, nullptr); - std::fill_n(tPtrs, 9, nullptr); + const vtkIdType& sliceOffset = this->TriadSliceOffset; + TrimmedEdgesCaseType trimmedEdgesCase = (row != 0) + ((slice != 0) << 1); + + // Initialize the triads and edge metadata. This simplifies the code. + std::fill_n(eMDPtrs.begin(), 9, nullptr); + std::fill_n(triadPtrs.begin(), 9, nullptr); // These portions of the bundle are always valid, with no boundary issues. - ePtrs[4] = EMD + (slice * dims[1] + row) * size; // current edge row - tPtrs[4] = triads + row * dims[0] + slice * sliceOffset; + eMDPtrs[4] = this->EdgeMetaData.data() + (slice * this->TriadDims[Y] + row); // current edge row + triadPtrs[4] = this->Triads.data() + row * this->TriadDims[X] + slice * sliceOffset; - ePtrs[5] = ePtrs[4] + size; // to the right of the current edge - tPtrs[5] = tPtrs[4] + dims[0]; + eMDPtrs[5] = eMDPtrs[4] + 1; // to the right of the current edge + triadPtrs[5] = triadPtrs[4] + this->TriadDims[X]; - ePtrs[7] = ePtrs[4] + dims[1] * size; // above the current edge - tPtrs[7] = tPtrs[4] + sliceOffset; + eMDPtrs[7] = eMDPtrs[4] + this->TriadDims[Y]; // above the current edge + triadPtrs[7] = triadPtrs[4] + sliceOffset; - ePtrs[8] = ePtrs[7] + size; // above and to the right of the current edge - tPtrs[8] = tPtrs[7] + dims[0]; + eMDPtrs[8] = eMDPtrs[7] + 1; // above and to the right of the current edge + triadPtrs[8] = triadPtrs[7] + this->TriadDims[X]; // May be near the -x,-y,-z boundaries. // If at origin of y-z plane. if (row != 0 && slice != 0) { - ePtrs[0] = ePtrs[4] - size - dims[1] * size; - tPtrs[0] = tPtrs[4] - dims[0] - sliceOffset; + eMDPtrs[0] = eMDPtrs[4] - 1 - this->TriadDims[Y]; + triadPtrs[0] = triadPtrs[4] - this->TriadDims[X] - sliceOffset; } if (slice != 0) // if not on -z boundary { - ePtrs[1] = ePtrs[4] - dims[1] * size; - tPtrs[1] = tPtrs[4] - sliceOffset; + eMDPtrs[1] = eMDPtrs[4] - this->TriadDims[Y]; + triadPtrs[1] = triadPtrs[4] - sliceOffset; - ePtrs[2] = ePtrs[4] + size - dims[1] * size; - tPtrs[2] = tPtrs[4] + dims[0] - sliceOffset; + eMDPtrs[2] = eMDPtrs[4] + 1 - this->TriadDims[Y]; + triadPtrs[2] = triadPtrs[4] + this->TriadDims[X] - sliceOffset; } if (row != 0) // if not on -y boundary { - ePtrs[3] = ePtrs[4] - size; - tPtrs[3] = tPtrs[4] - dims[0]; + eMDPtrs[3] = eMDPtrs[4] - 1; + triadPtrs[3] = triadPtrs[4] - this->TriadDims[X]; - ePtrs[6] = ePtrs[4] - size + dims[1] * size; - tPtrs[6] = tPtrs[4] - dims[0] + sliceOffset; + eMDPtrs[6] = eMDPtrs[4] - 1 + this->TriadDims[Y]; + triadPtrs[6] = triadPtrs[4] - this->TriadDims[X] + sliceOffset; } // Determine the trim over 3x3 bundle of metadata. This relies // on the earlier initialization of eMD. - vtkIdType* eMD; - xL = this->TriadDims[0]; - xR = 0; - for (auto i = 0; i < 9; ++i) + xMin = this->TriadDims[X]; + xMax = 0; + for (const auto& eMDPtr : eMDPtrs) { - eMD = ePtrs[i]; - if (eMD != nullptr) + if (eMDPtr != nullptr) { - xL = (eMD[3] < xL ? eMD[3] : xL); - xR = (eMD[4] > xR ? eMD[4] : xR); + xMin = std::min(xMin, eMDPtr->XMin); + xMax = std::max(xMax, eMDPtr->XMax); } } return trimmedEdgesCase; @@ -667,17 +659,18 @@ struct SurfaceNets // for each row of voxel cells. This avoids having to use a locator to merge // coincident points. The x-row iterator works across 3x3 bundles of // volume x-edges, with the current edge being processed in the center of - // the bundle. The edge bundle meta data is passed in to initialize the + // the bundle. The edge bundle metadata is passed in to initialize the // point ids. Note that while a 3x3 bundle is advanced, only six of the nine // values actually are used. (Empirical performance experiments show that // trying to optimize to six values actually slows execution slightly, // probably a compiler optimization hit.) - void InitRowIterator(vtkIdType* ePtrs[9], vtkIdType pIds[9]) + void InitRowIterator( + const std::array& eMDPtrs, std::array& pointIds) { - for (auto idx = 0; idx < 9; ++idx) + for (int idx = 0; idx < 9; ++idx) { - vtkIdType* eMD = ePtrs[idx]; - pIds[idx] = (eMD != nullptr ? eMD[0] : -1); + const auto eMDPtr = eMDPtrs[idx]; + pointIds[idx] = (eMDPtr != nullptr ? eMDPtr->NumPoints : -1); } } @@ -686,75 +679,75 @@ struct SurfaceNets // surrounding it, have points generated inside of them. Note that the // point ids refer to the nine edges in the 3x3 edge bundle centered around // the current edge being processed. - void AdvanceRowIterator( - vtkIdType i, TriadType* tPtrs[9], vtkIdType pIds[9], TrimmedEdgesCase trimmedEdgesCase) + void AdvanceRowIterator(vtkIdType i, std::array& triadPtrs, + std::array& pointIds, const TrimmedEdgesCaseType& trimmedEdgesCase) { - if (this->ProducesPoint(tPtrs[4][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[4][i])) { - pIds[4]++; + pointIds[4]++; } - if (this->ProducesPoint(tPtrs[5][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[5][i])) { - pIds[5]++; + pointIds[5]++; } - if (this->ProducesPoint(tPtrs[7][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[7][i])) { - pIds[7]++; + pointIds[7]++; } - if (this->ProducesPoint(tPtrs[8][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[8][i])) { - pIds[8]++; + pointIds[8]++; } switch (trimmedEdgesCase) { case 1: // if not on -y boundary { // do the checks only for 3 and 6 ids - if (this->ProducesPoint(tPtrs[3][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[3][i])) { - pIds[3]++; + pointIds[3]++; } - if (this->ProducesPoint(tPtrs[6][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[6][i])) { - pIds[6]++; + pointIds[6]++; } break; } case 2: // if not on -z boundary { // do the checks only for 1 and 2 ids - if (this->ProducesPoint(tPtrs[1][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[1][i])) { - pIds[1]++; + pointIds[1]++; } - if (this->ProducesPoint(tPtrs[2][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[2][i])) { - pIds[2]++; + pointIds[2]++; } break; } case 3: // if not on -y boundary and not on -z boundary { - // do the checks for 0, 1, 2, 3,6 ids - if (this->ProducesPoint(tPtrs[0][i])) + // do the checks for 0, 1, 2, 3, 6 ids + if (SurfaceNets::ProducesPoint(triadPtrs[0][i])) { - pIds[0]++; + pointIds[0]++; } - if (this->ProducesPoint(tPtrs[1][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[1][i])) { - pIds[1]++; + pointIds[1]++; } - if (this->ProducesPoint(tPtrs[2][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[2][i])) { - pIds[2]++; + pointIds[2]++; } - if (this->ProducesPoint(tPtrs[3][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[3][i])) { - pIds[3]++; + pointIds[3]++; } - if (this->ProducesPoint(tPtrs[6][i])) + if (SurfaceNets::ProducesPoint(triadPtrs[6][i])) { - pIds[6]++; + pointIds[6]++; } break; } @@ -765,17 +758,17 @@ struct SurfaceNets } // Given an i,j,k triad index, create a new point in the center of the - // triad. It is possible for some points to be generated outside of the + // triad. It is possible for some points to be generated outside the // actual image (i.e., in the padded boundary triads). The point is // generated in image space, later it will be transformed into world space // via vtkImageTransform. (Recall that the volume is padded out in the // x-y-z directions.) void GeneratePoint(vtkIdType ptId, vtkIdType i, vtkIdType j, vtkIdType k) { - float* x = this->NewPts + 3 * ptId; - x[0] = this->Min0 + static_cast(i) - 0.5; - x[1] = this->Min1 + static_cast(j) - 0.5; - x[2] = this->Min2 + static_cast(k) - 0.5; + float* p = this->NewPts + 3 * ptId; + p[X] = this->Min[X] + static_cast(i) - 0.5; + p[Y] = this->Min[Y] + static_cast(j) - 0.5; + p[Z] = this->Min[Z] + static_cast(k) - 0.5; } // Produce the output polygons (quads) for this triad. Note that at most @@ -787,46 +780,46 @@ struct SurfaceNets { template void operator()(OffsetsT* vtkNotUsed(offsets), ConnectivityT* conn, vtkIdType i, vtkIdType row, - vtkIdType slice, struct SurfaceNets* snet, TriadType triad, vtkIdType* pIds, + vtkIdType slice, SurfaceNets* snet, TriadType triad, std::array& pointIds, vtkIdType& quadId) { auto connRange = GetRange(conn); auto connIter = connRange.begin() + (quadId * 4); // Prepare to write scalar data. s0 is the triad origin. - T backgroundLabel = snet->BackgroundLabel; - T s0 = snet->GetVoxelForTriad(i, row, slice); + const T backgroundLabel = snet->BackgroundLabel; + const T s0 = snet->GetVoxelForTriad(i, row, slice); if (SurfaceNets::GenerateXYQuad(triad)) { - *connIter++ = pIds[4]; - *connIter++ = pIds[4] - 1; - *connIter++ = pIds[3] - 1; - *connIter++ = pIds[3]; // normal to the z triad edge + *connIter++ = pointIds[4]; + *connIter++ = pointIds[4] - 1; + *connIter++ = pointIds[3] - 1; + *connIter++ = pointIds[3]; // normal to the z triad edge - T s1 = snet->GetVoxelForTriad(i, row, slice + 1); + const T s1 = snet->GetVoxelForTriad(i, row, slice + 1); snet->WriteScalarTuple(s0, s1, backgroundLabel, quadId++); } if (SurfaceNets::GenerateXZQuad(triad)) { - *connIter++ = pIds[4]; - *connIter++ = pIds[1]; - *connIter++ = pIds[1] - 1; - *connIter++ = pIds[4] - 1; // normal to the y edge + *connIter++ = pointIds[4]; + *connIter++ = pointIds[1]; + *connIter++ = pointIds[1] - 1; + *connIter++ = pointIds[4] - 1; // normal to the y edge - T s1 = snet->GetVoxelForTriad(i, row + 1, slice); + const T s1 = snet->GetVoxelForTriad(i, row + 1, slice); snet->WriteScalarTuple(s0, s1, backgroundLabel, quadId++); } if (SurfaceNets::GenerateYZQuad(triad)) { - *connIter++ = pIds[4]; - *connIter++ = pIds[3]; - *connIter++ = pIds[0]; - *connIter++ = pIds[1]; // quad is normal to the x edge + *connIter++ = pointIds[4]; + *connIter++ = pointIds[3]; + *connIter++ = pointIds[0]; + *connIter++ = pointIds[1]; // quad is normal to the x edge - T s1 = snet->GetVoxelForTriad(i + 1, row, slice); + const T s1 = snet->GetVoxelForTriad(i + 1, row, slice); snet->WriteScalarTuple(s0, s1, backgroundLabel, quadId++); } @@ -837,68 +830,68 @@ struct SurfaceNets struct GenerateStencilImpl : public vtkCellArray::DispatchUtilities { template - void operator()(OffsetsT* offsets, ConnectivityT* conn, SurfaceNets* snet, EdgeCaseType eCase, - vtkIdType pIds[9], vtkIdType& sOffset) + void operator()(OffsetsT* offsets, ConnectivityT* conn, SurfaceNets* snet, + EdgeCaseType edgeCase, std::array& pointIds, vtkIdType& sOffset) { // The point on which the stencil operates - vtkIdType pId = pIds[4]; + const vtkIdType& pointId = pointIds[4]; auto offsetRange = GetRange(offsets); - auto offsetIter = offsetRange.begin() + pId; + auto offsetIter = offsetRange.begin() + pointId; auto connRange = GetRange(conn); auto connIter = connRange.begin() + sOffset; // Create the stencil. Note that for stencils with just one connection // (e.g., on the boundary of the image), the stencil point is "locked" // in place to prevent any motion to avoid shrinkage etc. - vtkIdType numEdges = snet->GetNumberOfStencilEdges(eCase); - *offsetIter++ = sOffset; + const vtkIdType numEdges = snet->GetNumberOfStencilEdges(edgeCase); + *offsetIter = sOffset; sOffset += numEdges; if (numEdges == 1) { - *connIter = pId; + *connIter = pointId; return; } // Create up to six stencil edges connecting the voxel edge face // neighbors. - const unsigned char* sEdges = snet->GetStencilEdges(eCase); + const unsigned char* stencilEdges = snet->GetStencilEdges(edgeCase); // -x face - if (sEdges[1]) + if (stencilEdges[1]) { - *connIter++ = pIds[4] - 1; + *connIter++ = pointIds[4] - 1; } // +x face - if (sEdges[2]) + if (stencilEdges[2]) { - *connIter++ = pIds[4] + 1; + *connIter++ = pointIds[4] + 1; } // -y face - if (sEdges[3]) + if (stencilEdges[3]) { - *connIter++ = pIds[3]; + *connIter++ = pointIds[3]; } // +y face - if (sEdges[4]) + if (stencilEdges[4]) { - *connIter++ = pIds[5]; + *connIter++ = pointIds[5]; } // -z face - if (sEdges[5]) + if (stencilEdges[5]) { - *connIter++ = pIds[1]; + *connIter++ = pointIds[1]; } // +z face - if (sEdges[6]) + if (stencilEdges[6]) { - *connIter++ = pIds[7]; + *connIter = pointIds[7]; } } // operator() }; // GenerateStencilImpl @@ -923,7 +916,8 @@ struct SurfaceNets // with boundary triads. T GetVoxelForTriad(vtkIdType i, vtkIdType row, vtkIdType slice) { - return this->Scalars[(slice - 1) * this->Inc2 + (row - 1) * this->Inc1 + (i - 1) * this->Inc0]; + return this + ->Scalars[(slice - 1) * this->Inc[Z] + (row - 1) * this->Inc[Y] + (i - 1) * this->Inc[X]]; } // Helper function writes the scalar 2-tuple. @@ -937,8 +931,8 @@ struct SurfaceNets std::swap(s0, s1); } - *scalars++ = s0; // write 2-tuple - *scalars++ = s1; + scalars[0] = s0; // write 2-tuple + scalars[1] = s1; } // WriteScalarTuple // The following are methods supporting the four passes of the @@ -1055,9 +1049,9 @@ void SurfaceNets::GenerateFaceStencils(unsigned char stencils[][7]) for (FaceCaseType faceCase = 0; faceCase < 64; ++faceCase) { stencils[faceCase][0] = 0; - for (auto faceNum = 1; faceNum <= 6; ++faceNum) + for (int faceNum = 1; faceNum <= 6; ++faceNum) { - stencils[faceCase][faceNum] = ((faceCase & (1 << (faceNum - 1))) > 0 ? 1 : 0); + stencils[faceCase][faceNum] = (faceCase & (1 << (faceNum - 1))) > 0; stencils[faceCase][0] += stencils[faceCase][faceNum]; // count the faces/stencil edges } } // for all face cases @@ -1071,7 +1065,7 @@ void SurfaceNets::GenerateEdgeStencils(int optLevel) { // Create the basic stencils without optimization. Basically, convert from // the 2^12 edge cases to the 2^6 stencil face cases. - constexpr int numEdgeCases = 4096; + static constexpr int numEdgeCases = 4096; for (auto edgeCase = 0; edgeCase < numEdgeCases; ++edgeCase) { this->StencilTable[edgeCase] = this->GetFaceCase(edgeCase); @@ -1091,7 +1085,7 @@ void SurfaceNets::GenerateEdgeStencils(int optLevel) // more intersections, then a different smoothing stencil should be associated // with the edge case (i.e., the stencil enhanced surface edge smoothing). unsigned char faceCount[6]; // the number of edge intersections on each of the voxels' six faces. - for (auto edgeCase = 0; edgeCase < numEdgeCases; ++edgeCase) + for (int edgeCase = 0; edgeCase < numEdgeCases; ++edgeCase) { unsigned char maxInts = this->CountFaceIntersections(edgeCase, faceCount); if (maxInts <= 2) @@ -1111,7 +1105,7 @@ void SurfaceNets::GenerateEdgeStencils(int optLevel) this->StencilTable[edgeCase] = faceStencilCase; // reference new stencil } -} // GenerareEdgeStencils +} // GenerateEdgeStencils // Implementations of the four passes of the surface nets boundary extraction // process. @@ -1133,24 +1127,24 @@ void SurfaceNets::ClassifyXEdges( { T s0, s1 = (*inPtr); // s1 is the first voxel value in the current row bool isLV0, isLV1 = lMap->IsLabelValue(s1); - vtkIdType numTriads = this->TriadDims[0]; - TriadType* rowTriadPtr = this->Triads + row * this->TriadDims[0] + slice * this->TriadSliceOffset; - vtkIdType* eMD = this->EdgeMetaData + (slice * this->TriadDims[1] + row) * this->EdgeMetaDataSize; - vtkIdType minInt = numTriads, maxInt = 0; + const vtkIdType& numTriads = this->TriadDims[X]; + TriadType* rowTriadPtr = + this->Triads.data() + row * this->TriadDims[X] + slice * this->TriadSliceOffset; + EdgeMetaDataType& eMD = this->EdgeMetaData[slice * this->TriadDims[Y] + row]; + vtkIdType xMin = numTriads, xMax = 0; + const vtkIdType numTriadsMinus1 = numTriads - 1, numTriadsMinus2 = numTriads - 2; // Run along the entire x-edge classifying the triad x axes. Be careful // with the padded triads: the 0th and (n-1) triads will not produce // intersections because they are in a padded voxel. Note that the ith // triad corresponds to the (i-1) image voxel. - for (vtkIdType i = 0; i < (numTriads - 1); ++i) + for (vtkIdType i = 0; i < numTriadsMinus1; ++i) { - // This handles the left-hand edge of the slice as well as setting up for - // the next triad. - TriadType* tPtr = rowTriadPtr + i; + // This handles the left-hand edge of the slice as well as setting up for the next triad. s0 = s1; isLV0 = isLV1; - if (i == (numTriads - 2)) + if (i == numTriadsMinus2) { // Edge of slice, voxel value s1 does not exist due to padding s1 = s0; @@ -1159,33 +1153,33 @@ void SurfaceNets::ClassifyXEdges( else { // Processing triads which are associated with voxels. - s1 = inPtr[i * this->Inc0]; - isLV1 = (s0 == s1 ? isLV0 : lMap->IsLabelValue(s1)); + s1 = inPtr[i * this->Inc[X]]; + isLV1 = s0 == s1 ? isLV0 : lMap->IsLabelValue(s1); } // Is the current triad origin vertex a label value? - TriadType tCase = (isLV0 ? SurfaceNets::Inside : SurfaceNets::Outside); + TriadType triad = isLV0 ? TriadClassification::Inside : TriadClassification::Outside; // Is the current x-edge split (i.e., different labels on each end). // Also update edge trim. if ((isLV0 || isLV1) && s0 != s1) { - tCase |= SurfaceNets::XIntersection; - minInt = (i < minInt ? i : minInt); - maxInt = i + 1; + triad |= TriadClassification::XIntersection; + xMin = std::min(xMin, i); + xMax = i + 1; } // If non-initialized (zero) state, update classification. - if (tCase > SurfaceNets::Outside) + if (triad > TriadClassification::Outside) { - this->SetTriadClassification(tPtr, tCase); + rowTriadPtr[i] = triad; } // if contour interacts with this triad } // for all triad-x-edges along this image x-edge - // The beginning and ending of intersections [xL,xR) along the edge is used + // The beginning and ending of intersections [xMin, xMax) along the edge is used // for computational trimming. - eMD[3] = minInt; - eMD[4] = (maxInt < numTriads ? maxInt : numTriads); + eMD.XMin = xMin; + eMD.XMax = std::min(xMax, numTriads); } // ClassifyXEdges //------------------------------------------------------------------------------ @@ -1196,45 +1190,48 @@ template void SurfaceNets::ClassifyYZEdges(T* inPtr, vtkIdType row, vtkIdType slice) { // Classify y- and z- triad edges. - // Triad cases(tPtr,tCase): this row, the next row (y-classification), and the + // Triad cases(triadPtr,tCase): this row, the next row (y-classification), and the // next slice (z-classification). - vtkIdType numTriads = this->TriadDims[0]; - TriadType* tPtr = this->Triads + row * numTriads + slice * this->TriadSliceOffset; - TriadType* tPtrY = tPtr + this->TriadDims[0]; - TriadType* tPtrZ = tPtr + this->TriadSliceOffset; + const vtkIdType& numTriads = this->TriadDims[X]; + TriadType* triadPtr = this->Triads.data() + row * numTriads + slice * this->TriadSliceOffset; + TriadType* triadPtrY = triadPtr + this->TriadDims[X]; + TriadType* triadPtrZ = triadPtr + this->TriadSliceOffset; // Edge metadata: this edge eMD, in the y-direction, and the z-direction. - vtkIdType* eMD = this->EdgeMetaData + (row + slice * this->TriadDims[1]) * this->EdgeMetaDataSize; - vtkIdType* eMDY = eMD + this->EdgeMetaDataSize; - vtkIdType* eMDZ = eMD + this->TriadDims[1] * this->EdgeMetaDataSize; + EdgeMetaDataType& eMD = this->EdgeMetaData[row + slice * this->TriadDims[Y]]; + EdgeMetaDataType& eMDY = (&eMD)[1]; + EdgeMetaDataType& eMDZ = (&eMD)[this->TriadDims[Y]]; + const vtkIdType numTriadsMinus1 = numTriads - 1, numTriadsMinus2 = numTriads - 2; // By default, all non-padded voxels on this volume-x-row will be // processed. However, based on the edge trim from the first pass, or the - // particulars of the data surrounding this edge, the edge trim (xL,xR) may - // be modified. - vtkIdType xL = 1, xR = (numTriads - 1); + // particulars of the data surrounding this edge, the edge trim (xMin,xMax) may be modified. + vtkIdType xMin = 1, xMax = numTriadsMinus1; // A quick check to determine whether this row of voxels needs processing // (this is a relatively common situation). If no x-edge intersections // exist (eMD[3]==numTriads) in this row or the rows to the right and // above, and the x-, y-, and z-rows are the same value, then the row can - // be skipped. Note that tPtr[1] is the first triad with an associated + // be skipped. Note that triadPtr[1] is the first triad with an associated // voxel value (due to padding). - bool xInts = !((eMD[3] >= numTriads) && (eMDY[3] >= numTriads) && (eMDZ[3] >= numTriads)); + const bool xInts = !(eMD.XMin >= numTriads && eMDY.XMin >= numTriads && eMDZ.XMin >= numTriads); if (!xInts) { - if (tPtr[1] == SurfaceNets::Outside) + if (triadPtr[1] == TriadClassification::Outside) { - if (tPtrY[1] == SurfaceNets::Outside && tPtrZ[1] == SurfaceNets::Outside) + if (triadPtrY[1] == TriadClassification::Outside && + triadPtrZ[1] == TriadClassification::Outside) { return; // This is a fairly common situation } } else // no volume-x-edge intersections, voxel values inside the same labeled region { - unsigned char yEdgeClass = this->ClassifyYEdge(inPtr, 1, tPtr[1], row, tPtrY[1]); - unsigned char zEdgeClass = this->ClassifyZEdge(inPtr, 1, tPtr[1], slice, tPtrZ[1]); - if (!yEdgeClass && !zEdgeClass) + TriadClassification yEdgeClass = + this->ClassifyYEdge(inPtr, 1, triadPtr[1], row, triadPtrY[1]); + TriadClassification zEdgeClass = + this->ClassifyZEdge(inPtr, 1, triadPtr[1], slice, triadPtrZ[1]); + if (yEdgeClass == TriadClassification::Outside && zEdgeClass == TriadClassification::Outside) { return; // no volume-x-edge ints, and voxel values are in the same material } @@ -1247,42 +1244,45 @@ void SurfaceNets::ClassifyYZEdges(T* inPtr, vtkIdType row, vtkIdType slice) // same material. If not, leave edge trim to default values. Otherwise, // reset the edge trim to the trim values determined in Pass1 / ClassifyXEdges // because this is the only place where voxel value changes can occur. - unsigned char yEdgeClass = this->ClassifyYEdge(inPtr, 1, tPtr[1], row, tPtrY[1]); - unsigned char zEdgeClass = this->ClassifyZEdge(inPtr, 1, tPtr[1], slice, tPtrZ[1]); - if (!yEdgeClass && !zEdgeClass) - { - xL = eMD[3]; - xL = (eMDY[3] < xL ? eMDY[3] : xL); - xL = (eMDZ[3] < xL ? eMDZ[3] : xL); - xL = (xL < 1 ? 1 : xL); - } - vtkIdType lastValue = numTriads - 2; - yEdgeClass = this->ClassifyYEdge(inPtr, lastValue, tPtr[lastValue], row, tPtrY[lastValue]); - zEdgeClass = this->ClassifyZEdge(inPtr, lastValue, tPtr[lastValue], slice, tPtrZ[lastValue]); - if (!yEdgeClass && !zEdgeClass) - { - xR = eMD[4]; - xR = (eMDY[4] > xR ? eMDY[4] : xR); - xR = (eMDZ[4] > xR ? eMDZ[4] : xR); - xR = (xR >= numTriads ? (numTriads - 1) : xR); + TriadClassification yEdgeClass = this->ClassifyYEdge(inPtr, 1, triadPtr[1], row, triadPtrY[1]); + TriadClassification zEdgeClass = + this->ClassifyZEdge(inPtr, 1, triadPtr[1], slice, triadPtrZ[1]); + if (yEdgeClass == TriadClassification::Outside && zEdgeClass == TriadClassification::Outside) + { + xMin = eMD.XMin; + xMin = std::min(xMin, eMDY.XMin); + xMin = std::min(xMin, eMDZ.XMin); + xMin = std::max(xMin, static_cast(1)); + } + const vtkIdType& lastValue = numTriadsMinus2; + yEdgeClass = + this->ClassifyYEdge(inPtr, lastValue, triadPtr[lastValue], row, triadPtrY[lastValue]); + zEdgeClass = + this->ClassifyZEdge(inPtr, lastValue, triadPtr[lastValue], slice, triadPtrZ[lastValue]); + if (yEdgeClass == TriadClassification::Outside && zEdgeClass == TriadClassification::Outside) + { + xMax = eMD.XMax; + xMax = std::max(xMax, eMDY.XMax); + xMax = std::max(xMax, eMDZ.XMax); + xMax = std::min(xMax, numTriadsMinus1); } } // Intersections along one of the volume-x-edges // Classify all the triad y- and z-edges, excluding the padded triads. - for (vtkIdType i = xL; i < xR; ++i) + for (vtkIdType i = xMin; i < xMax; ++i) { - TriadType tCase = tPtr[i]; - tCase |= this->ClassifyYEdge(inPtr, i, tPtr[i], row, tPtrY[i]); - tCase |= this->ClassifyZEdge(inPtr, i, tPtr[i], slice, tPtrZ[i]); - if (tPtr[i] != tCase) + TriadType tCase = triadPtr[i]; + tCase |= this->ClassifyYEdge(inPtr, i, triadPtr[i], row, triadPtrY[i]); + tCase |= this->ClassifyZEdge(inPtr, i, triadPtr[i], slice, triadPtrZ[i]); + if (triadPtr[i] != tCase) { - this->SetTriadClassification(tPtr + i, tCase); + triadPtr[i] = tCase; } } // for all voxels in this volume x-row // Update the edge trim - eMD[3] = xL; - eMD[4] = xR; + eMD.XMin = xMin; + eMD.XMax = xMax; } // ClassifyYZEdges // Process the voxels in a row, combining triads to determine the voxel @@ -1291,27 +1291,27 @@ void SurfaceNets::ClassifyYZEdges(T* inPtr, vtkIdType row, vtkIdType slice) // the code, a bit is set in the triad corresponding to the voxel // (ProducePoint). Because the triads from four rows are combined to produce // a voxel case, setting this bit could produce a race condition. Thus, the -// processing of voxels is 4-way interleaved/checkerboarded to avoid race +// processing of voxels is 4-way interleaved in a checkerboard way to avoid race // conditions i.e., 0<=whichEdge<4 with a group defined as the bundle of four // neighboring edges with origin (x,y,z) in the +y,+z directions. template void SurfaceNets::ProduceVoxelCases(vtkIdType group, int whichEdge, vtkIdType numRowPairs) { - vtkIdType numTriads = this->TriadDims[0]; - vtkIdType row = 2 * (group % numRowPairs) + (whichEdge % 2); - vtkIdType slice = 2 * (group / numRowPairs) + (whichEdge / 2); + const vtkIdType& numTriads = this->TriadDims[X]; + const vtkIdType row = 2 * (group % numRowPairs) + (whichEdge % 2); + const vtkIdType slice = 2 * (group / numRowPairs) + (whichEdge / 2); // Make sure we don't process bogus padded triads. - if (row >= (this->TriadDims[1] - 1) || slice >= (this->TriadDims[2] - 1)) + if (row >= (this->TriadDims[Y] - 1) || slice >= (this->TriadDims[Z] - 1)) { return; // don't process +y,+z padded boundaries } - // Grab the triad data for this row; and the meta data for this row, and the rows + // Grab the triad data for this row; and the metadata for this row, and the rows // that are needed to form a column of voxel cells. - vtkIdType xL, xR; // computational edge trim - vtkIdType* eMD = this->Get2x2EdgeTrim(row, slice, xL, xR); - TriadType* tPtr = this->Triads + row * numTriads + slice * this->TriadSliceOffset; + vtkIdType xMin, xMax; // computational edge trim + EdgeMetaDataType& eMD = *this->Get2x2EdgeTrim(row, slice, xMin, xMax); + TriadType* triadPtr = this->Triads.data() + row * numTriads + slice * this->TriadSliceOffset; // Loop across voxels in this row. We need to determine the case number of // each voxel from the seven triads that contribute to the complete edge @@ -1319,27 +1319,25 @@ void SurfaceNets::ProduceVoxelCases(vtkIdType group, int whichEdge, vtkIdType // because the smoothing stencils may include +/-x points before and after // the current voxel, the left edge trim is started one before the current // location along the voxel-x-edge. - xL = (xL <= 0 ? 0 : xL - 1); - for (vtkIdType i = xL; i < xR; ++i) + xMin = std::max(xMin - 1, static_cast(0)); + for (vtkIdType i = xMin; i < xMax; ++i) { - EdgeCaseType eCase = this->GetEdgeCase(tPtr + i); - if (eCase > 0) // then a point must be generated in this voxel + const EdgeCaseType edgeCase = this->GetEdgeCase(triadPtr + i); + if (edgeCase > 0) // then a point must be generated in this voxel { // Set the bit indicating the triad's voxel cell will generate a point - TriadType triad = tPtr[i]; - triad |= SurfaceNets::ProducePoint; - this->SetTriadClassification(tPtr + i, triad); + triadPtr[i] |= TriadClassification::ProducePoint; // Update metadata for this volume edge - eMD[0]++; // number of points generated - eMD[1] += this->GetNumberOfQuads(triad); // number of quads - eMD[2] += this->GetNumberOfStencilEdges(eCase); // stencil edges - } // if produces a point - } // for all triads on this row + eMD.NumPoints++; // number of points generated + eMD.NumQuads += SurfaceNets::GetNumberOfQuads(triadPtr[i]); // number of quads + eMD.NumStencilEdges += this->GetNumberOfStencilEdges(edgeCase); // stencil edges + } // if produces a point + } // for all triads on this row // Update the edge trim - eMD[3] = xL; - eMD[4] = xR; + eMD.XMin = xMin; + eMD.XMax = xMax; } // ProduceVoxelCases //------------------------------------------------------------------------------ @@ -1357,15 +1355,15 @@ void SurfaceNets::ConfigureOutput( // case, sum up the number of points, quads, and stencils generated for // each row. Note that to avoid race conditions, row processing is interleaved // (i.e., groups of four rows: +/-y +/-z volume edges). - vtkIdType numRows = this->TriadDims[1]; - vtkIdType numRowPairs = numRows / 2; - vtkIdType numSlices = this->TriadDims[2]; - vtkIdType numSlicePairs = numSlices / 2; - vtkIdType numGroups = numRowPairs * numSlicePairs; + const vtkIdType numRows = this->TriadDims[Y]; + const vtkIdType numRowPairs = numRows / 2; + const vtkIdType numSlices = this->TriadDims[Z]; + const vtkIdType numSlicePairs = numSlices / 2; + const vtkIdType numGroups = numRowPairs * numSlicePairs; // Process the four edges that compose a group in order. The four edges form // a 2x2 bundle, in the order (j,k),(j+1,k),(j,k+1),(j+1,k+1). - for (auto edgeNum = 0; edgeNum < 4; ++edgeNum) + for (int edgeNum = 0; edgeNum < 4; ++edgeNum) { // Edge groups consist of four neighboring volume x-edges (+/-y,+/-z). Process // in interleaving fashion (i.e., checkerboard) to avoid races (ProduceVoxelCases() @@ -1380,58 +1378,54 @@ void SurfaceNets::ConfigureOutput( }); } - // Begin prefix sum to determine the point, quad, and stencil number - // offsets for each row. - vtkIdType* eMD; - vtkIdType row, slice, numPts, numQuads, numSEdges; - - // Accumulate the total number of points, quads, and stencil edges - // across all the image x-rows. - vtkIdType numOutPts = 0; - vtkIdType numOutQuads = 0; - vtkIdType numOutSEdges = 0; + // Begin prefix sum to determine the point, quad, and stencil number offsets for each row. + EdgeMetaDataType tempEMD; + // Accumulate the total number of points, quads, and stencil edges across all the image x-rows. + EdgeMetaDataType outputEMD; // Prefix sum to build offsets into the output points, quads, and // stencils. We process all edge metadata. - for (slice = 0; slice < numSlices; ++slice) + for (vtkIdType slice = 0; slice < numSlices; ++slice) { - for (row = 0; row < numRows; ++row) + const vtkIdType sliceOffset = slice * this->TriadDims[Y]; + for (vtkIdType row = 0; row < numRows; ++row) { - eMD = this->EdgeMetaData + (slice * this->TriadDims[1] + row) * this->EdgeMetaDataSize; - numPts = eMD[0]; - numQuads = eMD[1]; - numSEdges = eMD[2]; + EdgeMetaDataType& eMD = this->EdgeMetaData[sliceOffset + row]; + tempEMD.NumPoints = eMD.NumPoints; + tempEMD.NumQuads = eMD.NumQuads; + tempEMD.NumStencilEdges = eMD.NumStencilEdges; - eMD[0] = numOutPts; - eMD[1] = numOutQuads; - eMD[2] = numOutSEdges; + eMD.NumPoints = outputEMD.NumPoints; + eMD.NumQuads = outputEMD.NumQuads; + eMD.NumStencilEdges = outputEMD.NumStencilEdges; - numOutPts += numPts; - numOutQuads += numQuads; - numOutSEdges += numSEdges; + outputEMD.NumPoints += tempEMD.NumPoints; + outputEMD.NumQuads += tempEMD.NumQuads; + outputEMD.NumStencilEdges += tempEMD.NumStencilEdges; } // for rows in this slice } // for slices // Output can now be allocated. - if (numOutPts > 0) + if (outputEMD.NumPoints > 0) { // Points, which are floats - newPts->SetNumberOfPoints(numOutPts); + newPts->SetNumberOfPoints(outputEMD.NumPoints); vtkFloatArray* fPts = static_cast(newPts->GetData()); this->NewPts = fPts->GetPointer(0); // Boundaries, a set of quads contained in vtkCellArray newQuads->UseFixedSizeDefaultStorage(4); - newQuads->ResizeExact(numOutQuads, 4 * numOutQuads); + newQuads->ResizeExact(outputEMD.NumQuads, 4 * outputEMD.NumQuads); this->NewQuads = newQuads; // Scalars, which are of type T and 2-components - newScalars->SetNumberOfTuples(numOutQuads); + newScalars->SetNumberOfTuples(outputEMD.NumQuads); this->NewScalars = newScalars->GetPointer(0); // Smoothing stencils, which are represented by a vtkCellArray - stencils->ResizeExact(numOutPts, numOutSEdges); - stencils->Dispatch(FinalizeStencilsOffsetsImpl{}, numOutPts, numOutSEdges); + stencils->ResizeExact(outputEMD.NumPoints, outputEMD.NumStencilEdges); + stencils->Dispatch( + FinalizeStencilsOffsetsImpl{}, outputEMD.NumPoints, outputEMD.NumStencilEdges); this->NewStencils = stencils; } } // ConfigureOutput @@ -1446,22 +1440,22 @@ void SurfaceNets::ConfigureOutput( // to generate points, quads, and stencils, the point ids are determined by // advancing the starting point ids from the current triad row, as well as the // rows immediately surrounding the current row (i.e., all those rows to -// to which stencil edges connect to, as well as generated quads). This forms +// which stencil edges connect to, as well as generated quads). This forms // a 3x3 bundle of volume edges / voxel rows centered on the current x-row. // The edge trim is used to reduce the amount of computation to perform. template void SurfaceNets::GenerateOutput(vtkIdType row, vtkIdType slice) { - // This volume edge's meta data, and the neighboring edge. - vtkIdType* eMD = this->EdgeMetaData + (slice * this->TriadDims[1] + row) * this->EdgeMetaDataSize; - vtkIdType* eMDNei = eMD + this->EdgeMetaDataSize; + // This volume edge's metadata, and the neighboring edge. + const EdgeMetaDataType& eMD = this->EdgeMetaData[slice * this->TriadDims[Y] + row]; + const EdgeMetaDataType& eMDNei = (&eMD)[1]; // Return if there is nothing to do (i.e., no points, quads or stencils to // generate). After the prefix sum in Pass3, eMD[0] is the starting number // of the points generated along this volume x-edge. So if (eMDNei[0]-eMD[0])<=0 // nothing is produced along the edge (i.e., no points generated) so it can be // skipped. - if (eMDNei[0] <= eMD[0]) + if (eMDNei.NumPoints <= eMD.NumPoints) { return; } @@ -1469,54 +1463,55 @@ void SurfaceNets::GenerateOutput(vtkIdType row, vtkIdType slice) // Given a volume x-edge to process (defined by [row,slice]), determine the // trim edges and the 3x3 row triad cases centered around the current // x-edge. - vtkIdType xL, xR; // computational trim edges - TriadType* tPtrs[9]; // pointers to the 3x3 bundle of row triad cases - vtkIdType* eMDPtrs[9]; // pointers to the 3x3 bundle of edge meta data - TrimmedEdgesCase trimmedEdgesCase = this->Get3x3EdgeTrim(row, slice, xL, xR, eMDPtrs, tPtrs); - TriadType* tPtr = tPtrs[4]; // triad pointers for current row + vtkIdType xMin, xMax; // computational trim edges + std::array triadPtrs; // pointers to the 3x3 bundle of row triad cases + std::array eMDPtrs; // pointers to the 3x3 bundle of edge metadata + TrimmedEdgesCaseType trimmedEdgesCase = + this->Get3x3EdgeTrim(row, slice, xMin, xMax, eMDPtrs, triadPtrs); + TriadType* triadPtr = triadPtrs[4]; // triad pointers for current row // Initialize the point numbering process using a row iterator. This uses // the information gathered from the prefix sum (Pass3) and contained in - // the edge meta data to obtain point numbers/ids, and the number/size of - // quads and stencils. The pIds[9] are the current starting point ids for + // the edge metadata to obtain point numbers/ids, and the number/size of + // quads and stencils. The nine pointIds are the current starting point ids for // rows surrounding the current edge (in total, a 3x3 stencil, which - // includes in the center of the stencil, the current edge). The point ids - // are initialized with the edge meta data, and advanced as a function of - // the nine triads[9] along the nine edges. - vtkIdType pIds[9]; - this->InitRowIterator(eMDPtrs, pIds); - vtkIdType quadId = eMD[1]; // starting quad id for this row - vtkIdType sOffset = eMD[2]; // starting stencil offset for this row + // includes in the center of the stencil, the current edge). The pointIds + // are initialized with the edge metadata, and advanced as a function of + // the nine triadPtrs along the nine edges. + std::array pointIds; + this->InitRowIterator(eMDPtrs, pointIds); + vtkIdType quadId = eMD.NumQuads; // starting quad id for this row + vtkIdType sOffset = eMD.NumStencilEdges; // starting stencil offset for this row // Now traverse all the voxels in this row, generating points, quads, // stencils, and optional scalar data. Points are only generated from the // current row; quad segments from the triad x-y-z edges; and stencils // connecting a voxel's point to six possible face neighbors. - for (auto i = xL; i < xR; ++i) + for (vtkIdType i = xMin; i < xMax; ++i) { // See if any points or quads are to be generated in this voxel. - TriadType triad = tPtr[i]; - if (this->ProducesPoint(triad)) + const TriadType& triad = triadPtr[i]; + if (SurfaceNets::ProducesPoint(triad)) { // Output a point in the center of the voxel. - this->GeneratePoint(pIds[4], i, row, slice); + this->GeneratePoint(pointIds[4], i, row, slice); // Produce quads if necessary. - if (this->ProducesQuad(triad)) + if (SurfaceNets::ProducesQuad(triad)) { - this->NewQuads->Dispatch(GenerateQuadsImpl{}, i, row, slice, this, triad, pIds, quadId); + this->NewQuads->Dispatch(GenerateQuadsImpl{}, i, row, slice, this, triad, pointIds, quadId); } // If a point is generated, then smoothing stencils are required (i.e., // stencils indicate how the generated point is connected to other // points). Up to six connections corresponding to six face neighbors // may be generated. - EdgeCaseType eCase = this->GetEdgeCase(tPtr + i); - this->NewStencils->Dispatch(GenerateStencilImpl{}, this, eCase, pIds, sOffset); - } // if need to generate a point + const EdgeCaseType edgeCase = this->GetEdgeCase(triadPtr + i); + this->NewStencils->Dispatch(GenerateStencilImpl{}, this, edgeCase, pointIds, sOffset); + } // if you need to generate a point // Need to increment the point ids. - this->AdvanceRowIterator(i, tPtrs, pIds, trimmedEdgesCase); + this->AdvanceRowIterator(i, triadPtrs, pointIds, trimmedEdgesCase); } // for all triads on this row } // GenerateOutput @@ -1527,7 +1522,7 @@ void SurfaceNets::GenerateOutput(vtkIdType row, vtkIdType slice) struct NetsWorker { // PASS 1: Process all triads on the given x-rows to classify triad - // x-axis. Interface to vtkSMPTools::For(). Note that triad row i + // x-axis. Interface to vtkSMPTools::For(). Note that triad row "i" // corresponds to image row (i-1) (due to padding). Also note that looking // up labels can be expensive, so a vtkLabelMapLookup is used to accelerate // the lookup process. Note that the parallel for (vtkSMPTools::For()) is @@ -1550,7 +1545,7 @@ struct NetsWorker void operator()(vtkIdType slice, vtkIdType endSlice) { vtkLabelMapLookup* lMap = this->LMap.Local(); - T *rowPtr, *slicePtr = this->Algo->Scalars + (slice - 1) * this->Algo->Inc2; + T *rowPtr, *slicePtr = this->Algo->Scalars + (slice - 1) * this->Algo->Inc[Z]; // Process slice-by-slice. Note that the bottom and top slices are not // processed (they have been 0-initialized). @@ -1558,30 +1553,30 @@ struct NetsWorker { // Process only triads associated with voxels (i.e., no padded voxels). rowPtr = slicePtr; - for (vtkIdType row = 1; row < (this->Algo->TriadDims[1] - 1); ++row) + for (vtkIdType row = 1, rowEnd = this->Algo->TriadDims[Y] - 1; row < rowEnd; ++row) { this->Algo->ClassifyXEdges(rowPtr, row, slice, lMap); - rowPtr += this->Algo->Inc1; + rowPtr += this->Algo->Inc[Y]; } // for all rows in this slice - slicePtr += this->Algo->Inc2; + slicePtr += this->Algo->Inc[Z]; } // for all slices in this batch } void Reduce() { - // Delete all of the label map lookups + // Delete all the label map lookups for (auto lmItr = this->LMap.begin(); lmItr != this->LMap.end(); ++lmItr) { delete *lmItr; } // over all threads - // Note that unlike vtkSurfaceNets2D, the edge meta data has been + // Note that unlike vtkSurfaceNets2D, the edge metadata has been // initialized to a "do not process" state so nothing else needs be // done. } }; // Pass1 dispatch // PASS 2: Process all voxels on the given x-rows to classify triad y-z-axes, - // and classify voxels. Interface to vtkSMPTools::For(). Note that triad row i + // and classify voxels. Interface to vtkSMPTools::For(). Note that triad row "i" // corresponds to image row (i-1) (i.e., the triads pad out the volume). template struct Pass2 @@ -1590,7 +1585,7 @@ struct NetsWorker Pass2(SurfaceNets* algo) { this->Algo = algo; } void operator()(vtkIdType slice, vtkIdType endSlice) { - T *rowPtr, *slicePtr = this->Algo->Scalars + (slice - 1) * this->Algo->Inc2; + T *rowPtr, *slicePtr = this->Algo->Scalars + (slice - 1) * this->Algo->Inc[Z]; // Process slice-by-slice. Note that the bottom and top slices are not // processed (they have been 0-initialized). @@ -1598,12 +1593,12 @@ struct NetsWorker { // Process only triads associated with voxels (i.e., no padded voxels). rowPtr = slicePtr; - for (vtkIdType row = 1; row < (this->Algo->TriadDims[1] - 1); ++row) + for (vtkIdType row = 1, rowEnd = this->Algo->TriadDims[Y] - 1; row < rowEnd; ++row) { this->Algo->ClassifyYZEdges(rowPtr, row, slice); - rowPtr += this->Algo->Inc1; + rowPtr += this->Algo->Inc[Y]; } // for all rows in this slice - slicePtr += this->Algo->Inc2; + slicePtr += this->Algo->Inc[Z]; } // for all slices in this batch } }; // Pass2 dispatch @@ -1628,24 +1623,23 @@ struct NetsWorker { // Note that there is no need to process the last (padded) slice. SurfaceNets* algo = this->Algo; - vtkIdType sliceOffset = algo->TriadDims[1] * algo->EdgeMetaDataSize; - vtkIdType* eMD0 = algo->EdgeMetaData + (slice * sliceOffset); - vtkIdType* eMD1 = eMD0 + sliceOffset; + EdgeMetaDataType* eMD0Ptr = algo->EdgeMetaData.data() + (slice * algo->TriadDims[Y]); + EdgeMetaDataType* eMD1Ptr = eMD0Ptr + algo->TriadDims[Y]; for (; slice < endSlice; ++slice) { // Make sure that some data is actually generated on this slice. Skip entire // slice if possible. - if (eMD1[0] > eMD0[0]) // Are points generated? + if (eMD1Ptr->NumPoints > eMD0Ptr->NumPoints) // Are points generated? { // Note that there is no need to process the last (padded) triads on this row. - for (vtkIdType row = 0; row < (this->Algo->TriadDims[1] - 1); ++row) + for (vtkIdType row = 0, rowEnd = this->Algo->TriadDims[Y] - 1; row < rowEnd; ++row) { this->Algo->GenerateOutput(row, slice); } // for all rows in this slice } // if points are generated - eMD0 = eMD1; - eMD1 = eMD0 + sliceOffset; + eMD0Ptr = eMD1Ptr; + eMD1Ptr = eMD0Ptr + algo->TriadDims[Y]; } // for all slices in this batch } }; // Pass4 dispatch @@ -1660,12 +1654,12 @@ struct NetsWorker using ValueType = vtk::GetAPIType; auto newScalars = ST::FastDownCast(newScalarsDataArray); - // The update extent may be different than the extent of the image. + // The update extent may be different from the extent of the image. // The only problem with using the update extent is that one or two // sources enlarge the update extent. This behavior is slated to be // eliminated. - vtkIdType incs[3]; - input->GetIncrements(incs); + vtkIdType increments[3]; + input->GetIncrements(increments); int* ext = input->GetExtent(); // Capture information for subsequent processing. Make sure that we are @@ -1677,15 +1671,15 @@ struct NetsWorker vtkErrorWithObjectMacro(self, "Expecting 3D data (volume)."); } - algo.Min0 = updateExt[0]; - algo.Max0 = updateExt[1]; - algo.Inc0 = incs[0]; - algo.Min1 = updateExt[2]; - algo.Max1 = updateExt[3]; - algo.Inc1 = incs[1]; - algo.Min2 = updateExt[4]; - algo.Max2 = updateExt[5]; - algo.Inc2 = incs[2]; + algo.Min[X] = updateExt[0]; + algo.Max[X] = updateExt[1]; + algo.Inc[X] = increments[X]; + algo.Min[Y] = updateExt[2]; + algo.Max[Y] = updateExt[3]; + algo.Inc[Y] = increments[Y]; + algo.Min[Z] = updateExt[4]; + algo.Max[Z] = updateExt[5]; + algo.Inc[Z] = increments[Z]; // Now allocate the working arrays. The Triads array tracks case# for // each voxel triad (and the corresponding voxel). Note that each input @@ -1695,46 +1689,48 @@ struct NetsWorker // but be aware that the triads on the boundaries of the volume are // treated specially. Note that the allocation of the triads initializes // them to zero; we depend on this as the initial triad classification. - algo.Dims[0] = algo.Max0 - algo.Min0 + 1; - algo.Dims[1] = algo.Max1 - algo.Min1 + 1; - algo.Dims[2] = algo.Max2 - algo.Min2 + 1; - algo.TriadDims[0] = algo.Dims[0] + 2; // padded in the +/-x direction - algo.TriadDims[1] = algo.Dims[1] + 2; // padded in the +/-y direction - algo.TriadDims[2] = algo.Dims[2] + 2; // padded in the +/-z direction - algo.TriadSliceOffset = algo.TriadDims[0] * algo.TriadDims[1]; - algo.Triads = new TriadType[algo.TriadSliceOffset * algo.TriadDims[2]](); - - // Also allocate the characterization (metadata) array for all of the x + algo.Dims[X] = algo.Max[X] - algo.Min[X] + 1; + algo.Dims[Y] = algo.Max[Y] - algo.Min[Y] + 1; + algo.Dims[Z] = algo.Max[Z] - algo.Min[Z] + 1; + algo.TriadDims[X] = algo.Dims[X] + 2; // padded in the +/-x direction + algo.TriadDims[Y] = algo.Dims[Y] + 2; // padded in the +/-y direction + algo.TriadDims[Z] = algo.Dims[Z] + 2; // padded in the +/-z direction + algo.TriadSliceOffset = algo.TriadDims[X] * algo.TriadDims[Y]; + algo.Triads.resize(static_cast(algo.TriadSliceOffset * algo.TriadDims[Z]), 0); + + // Also allocate the characterization (metadata) array for all the x // volume edges, including the padded out triads. So the x-edge metadata is // defined on the y-z plane. This edge metadata array (often referred to // as eMD[5]) tracks 0) the number points added along each x-row; as well // as 1) the number of quad primitives; 2) the number of stencil edges; // and the 3) xMin_i and 4) xMax_i (minimum index of first intersection, // maximum index of intersection for row i, so-called trim edges used for - // computational trimming). Note that the edge meta data eMD[0-2] is zero + // computational trimming). Note that the edge metadata eMD[0-2] is zero // initialized, while eMD[3,4] is initialized to a "do not process" state // which will likely by updated in pass1, pass2, or pass3. - algo.NumberOfEdges = algo.TriadDims[1] * algo.TriadDims[2]; // y-z plane of edges - algo.EdgeMetaData = new vtkIdType[algo.NumberOfEdges * algo.EdgeMetaDataSize](); - vtkSMPTools::For(0, algo.NumberOfEdges, + + // y-z plane of edges + algo.EdgeMetaData.resize(static_cast(algo.TriadDims[Y] * algo.TriadDims[Z])); + vtkSMPTools::For(0, static_cast(algo.EdgeMetaData.size()), [&](vtkIdType begin, vtkIdType end) { - for (vtkIdType eNum = begin; eNum < end; ++eNum) + for (vtkIdType edgeId = begin; edgeId < end; ++edgeId) { - vtkIdType* eMD = algo.EdgeMetaData + eNum * algo.EdgeMetaDataSize; - eMD[3] = algo.TriadDims[0]; - eMD[4] = 0; + EdgeMetaDataType& eMD = algo.EdgeMetaData[edgeId]; + eMD.XMin = algo.TriadDims[0]; + eMD.XMax = 0; } }); // Compute the starting offset location for scalar data. We may be operating // on a part of the volume. ValueType* scalars = static_cast(scalarsArray->GetPointer(0)); - algo.Scalars = scalars + incs[0] * (updateExt[0] - ext[0]) + incs[1] * (updateExt[2] - ext[2]) + - incs[2] * (updateExt[4] - ext[4]) + self->GetArrayComponent(); + algo.Scalars = scalars + increments[X] * (updateExt[0] - ext[0]) + + increments[Y] * (updateExt[2] - ext[2]) + increments[Z] * (updateExt[4] - ext[4]) + + self->GetArrayComponent(); // This algorithm executes just once no matter how many contour/label - // values, requiring a fast lookup as to whether a data/voxel value is a + // values, requiring a fast lookup whether a data/voxel value is a // contour value, or should be considered part of the background. In // Pass1, instances of vtkLabelMapLookup are created (per thread) // which performs the fast label lookup. @@ -1750,12 +1746,12 @@ struct NetsWorker // a time. Empirically this performs a little better than processing each // edge separately. Pass1 pass1(&algo); - vtkSMPTools::For(1, algo.TriadDims[2] - 1, pass1); + vtkSMPTools::For(1, algo.TriadDims[Z] - 1, pass1); // Classify the triad y-z-edges; finalize the triad classification. // Note that the last padded z-slice of triads is not processed. Pass2 pass2(&algo); - vtkSMPTools::For(1, algo.TriadDims[2] - 1, pass2); + vtkSMPTools::For(1, algo.TriadDims[Z] - 1, pass2); // Prefix sum to determine the size and character of the output, and // then allocate it. @@ -1765,11 +1761,10 @@ struct NetsWorker // data slice-by-slice. Note that the last (padded) slice is not // processed. Pass4 pass4(&algo); - vtkSMPTools::For(0, algo.TriadDims[2] - 1, pass4); + vtkSMPTools::For(0, algo.TriadDims[Z] - 1, pass4); - // Clean up and return - delete[] algo.Triads; - delete[] algo.EdgeMetaData; + algo.Triads.clear(); + algo.EdgeMetaData.clear(); } }; // NetsWorker @@ -1866,7 +1861,7 @@ struct ScalarsWorker // Functor to drive the threaded conversion of a quad output mesh to // a different type (i.e., triangles). -struct TransformMesh +struct TransformMeshToTris { vtkFloatArray* Points; vtkCellArray* QuadMesh; @@ -1878,14 +1873,16 @@ struct TransformMesh // Each thread has a cell array iterator to avoid constant allocation. vtkSMPThreadLocalObject TLIdList; - TransformMesh(vtkFloatArray* pts, vtkCellArray* qMesh, int triStrategy) + TransformMeshToTris(vtkFloatArray* pts, vtkCellArray* qMesh, int triStrategy) : Points(pts) , QuadMesh(qMesh) - , OutputMesh(nullptr) + , OutputMesh(vtkSmartPointer::New()) , TriStrategy(triStrategy) - , NumOutputCells(0) - , OutputConnSize(0) + , NumOutputCells(2 * qMesh->GetNumberOfCells()) + , OutputConnSize(6 * qMesh->GetNumberOfCells()) { + this->OutputMesh->UseFixedSizeDefaultStorage(3); + this->OutputMesh->ResizeExact(this->NumOutputCells, this->OutputConnSize); } void operator()(vtkIdType cellId, vtkIdType endCellId) @@ -1896,7 +1893,7 @@ struct TransformMesh double x0[3], x1[3], x2[3], x3[3]; vtkIdType conn[4]; int triStrategy = this->TriStrategy; - bool d02 = true; // diagonal 02 + bool d02; // diagonal 02 for (; cellId < endCellId; ++cellId) { @@ -1909,15 +1906,13 @@ struct TransformMesh if (triStrategy == vtkSurfaceNets3D::TRIANGULATION_MIN_EDGE) { - d02 = (vtkMath::Distance2BetweenPoints(x0, x2) < vtkMath::Distance2BetweenPoints(x1, x3) - ? true - : false); + d02 = vtkMath::Distance2BetweenPoints(x0, x2) < vtkMath::Distance2BetweenPoints(x1, x3); } else if (triStrategy == vtkSurfaceNets3D::TRIANGULATION_MIN_AREA) { double a02 = vtkTriangle::TriangleArea(x0, x2, x1) + vtkTriangle::TriangleArea(x0, x2, x3); double a13 = vtkTriangle::TriangleArea(x1, x3, x0) + vtkTriangle::TriangleArea(x1, x3, x2); - d02 = (a02 < a13 ? true : false); + d02 = a02 < a13; } else // if ( triStrategy == vtkSurfaceNets3D::TRIANGULATION_GREEDY ) { @@ -1944,23 +1939,6 @@ struct TransformMesh this->OutputMesh->Dispatch(ConvertToTrisImpl{}, cellId, conn); } // over this batch of cells } -}; // TransformMesh - -// Transform quad mesh to triangles. Also transform the cell data associated with the -// quads. -struct TransformMeshToTris : public TransformMesh -{ - // This constructor sets up the transformation process. - TransformMeshToTris(vtkFloatArray* pts, vtkCellArray* qMesh, int triStrategy) - : TransformMesh(pts, qMesh, triStrategy) - { - this->NumOutputCells = 2 * qMesh->GetNumberOfCells(); - this->OutputConnSize = 6 * qMesh->GetNumberOfCells(); - - this->OutputMesh = vtkSmartPointer::New(); - this->OutputMesh->UseFixedSizeDefaultStorage(3); - this->OutputMesh->ResizeExact(this->NumOutputCells, this->OutputConnSize); - } }; // TransformMeshToTris // This function is used to triangulate the output quads produced by the -- GitLab From b07d9f8a7a5c21d90475eafc31c99a44aa228d11 Mon Sep 17 00:00:00 2001 From: Spiros Tsalikis Date: Sat, 23 Nov 2024 23:57:02 -0500 Subject: [PATCH 2/5] vtkSurfaceNets3D: Convert AdvanceRowIterator to be branch-less Considering pointIds is a local variable, branch-less code is faster. --- Filters/Core/vtkSurfaceNets3D.cxx | 65 +++++++------------------------ 1 file changed, 13 insertions(+), 52 deletions(-) diff --git a/Filters/Core/vtkSurfaceNets3D.cxx b/Filters/Core/vtkSurfaceNets3D.cxx index 49620c3ac68a..798bc04b60e1 100644 --- a/Filters/Core/vtkSurfaceNets3D.cxx +++ b/Filters/Core/vtkSurfaceNets3D.cxx @@ -682,73 +682,34 @@ struct SurfaceNets void AdvanceRowIterator(vtkIdType i, std::array& triadPtrs, std::array& pointIds, const TrimmedEdgesCaseType& trimmedEdgesCase) { - if (SurfaceNets::ProducesPoint(triadPtrs[4][i])) - { - pointIds[4]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[5][i])) - { - pointIds[5]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[7][i])) - { - pointIds[7]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[8][i])) - { - pointIds[8]++; - } + pointIds[4] += SurfaceNets::ProducesPoint(triadPtrs[4][i]); + pointIds[5] += SurfaceNets::ProducesPoint(triadPtrs[5][i]); + pointIds[7] += SurfaceNets::ProducesPoint(triadPtrs[7][i]); + pointIds[8] += SurfaceNets::ProducesPoint(triadPtrs[8][i]); switch (trimmedEdgesCase) { case 1: // if not on -y boundary { // do the checks only for 3 and 6 ids - if (SurfaceNets::ProducesPoint(triadPtrs[3][i])) - { - pointIds[3]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[6][i])) - { - pointIds[6]++; - } + pointIds[3] += SurfaceNets::ProducesPoint(triadPtrs[3][i]); + pointIds[6] += SurfaceNets::ProducesPoint(triadPtrs[6][i]); break; } case 2: // if not on -z boundary { // do the checks only for 1 and 2 ids - if (SurfaceNets::ProducesPoint(triadPtrs[1][i])) - { - pointIds[1]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[2][i])) - { - pointIds[2]++; - } + pointIds[1] += SurfaceNets::ProducesPoint(triadPtrs[1][i]); + pointIds[2] += SurfaceNets::ProducesPoint(triadPtrs[2][i]); break; } case 3: // if not on -y boundary and not on -z boundary { // do the checks for 0, 1, 2, 3, 6 ids - if (SurfaceNets::ProducesPoint(triadPtrs[0][i])) - { - pointIds[0]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[1][i])) - { - pointIds[1]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[2][i])) - { - pointIds[2]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[3][i])) - { - pointIds[3]++; - } - if (SurfaceNets::ProducesPoint(triadPtrs[6][i])) - { - pointIds[6]++; - } + pointIds[0] += SurfaceNets::ProducesPoint(triadPtrs[0][i]); + pointIds[1] += SurfaceNets::ProducesPoint(triadPtrs[1][i]); + pointIds[2] += SurfaceNets::ProducesPoint(triadPtrs[2][i]); + pointIds[3] += SurfaceNets::ProducesPoint(triadPtrs[3][i]); + pointIds[6] += SurfaceNets::ProducesPoint(triadPtrs[6][i]); break; } case 0: -- GitLab From dc29a98f0066dca227d390a99c11e4e55b8fc56d Mon Sep 17 00:00:00 2001 From: Spiros Tsalikis Date: Thu, 16 Jan 2025 14:18:51 -0500 Subject: [PATCH 3/5] vtkSurfaceNets3D: Fix orientation consistency and scalar data handling - Ensured quads are generated with consistent orientation based on voxel label comparison. - Added logic to swap scalar values and vertex indices (c1 and c3) when the voxel labels indicate a reversed order. - Updated `GenerateXYQuad`, `GenerateXZQuad`, and `GenerateYZQuad` blocks to account for proper orientation checks and label prioritization. - Improved clarity in quad connectivity assignments to ensure consistent winding order. - Updated `WriteScalarTuple` to simply write the tuple - Generalized changes to maintain compatibility across all quad generation planes (XY, XZ, YZ). These changes improve the reliability of quad generation by enforcing consistent labeling, winding order, and surface normals, reducing errors in downstream processing. --- Filters/Core/vtkSurfaceNets3D.cxx | 85 +++++++++++++++++++++---------- 1 file changed, 57 insertions(+), 28 deletions(-) diff --git a/Filters/Core/vtkSurfaceNets3D.cxx b/Filters/Core/vtkSurfaceNets3D.cxx index 798bc04b60e1..8c6b5c800207 100644 --- a/Filters/Core/vtkSurfaceNets3D.cxx +++ b/Filters/Core/vtkSurfaceNets3D.cxx @@ -749,39 +749,75 @@ struct SurfaceNets // Prepare to write scalar data. s0 is the triad origin. const T backgroundLabel = snet->BackgroundLabel; - const T s0 = snet->GetVoxelForTriad(i, row, slice); + const T s0Origin = snet->GetVoxelForTriad(i, row, slice); if (SurfaceNets::GenerateXYQuad(triad)) { - *connIter++ = pointIds[4]; - *connIter++ = pointIds[4] - 1; - *connIter++ = pointIds[3] - 1; - *connIter++ = pointIds[3]; // normal to the z triad edge + const vtkIdType& c0 = pointIds[4]; + vtkIdType c1 = pointIds[4] - 1; + const vtkIdType c2 = pointIds[3] - 1; + vtkIdType c3 = pointIds[3]; + + T s0 = s0Origin; + T s1 = snet->GetVoxelForTriad(i, row, slice + 1); + if (s0 == backgroundLabel || (s1 != backgroundLabel && s0 > s1)) + { + std::swap(s0, s1); + std::swap(c1, c3); + } + + *connIter++ = c0; + *connIter++ = c1; + *connIter++ = c2; + *connIter++ = c3; - const T s1 = snet->GetVoxelForTriad(i, row, slice + 1); - snet->WriteScalarTuple(s0, s1, backgroundLabel, quadId++); + snet->WriteScalarTuple(s0, s1, quadId++); } if (SurfaceNets::GenerateXZQuad(triad)) { - *connIter++ = pointIds[4]; - *connIter++ = pointIds[1]; - *connIter++ = pointIds[1] - 1; - *connIter++ = pointIds[4] - 1; // normal to the y edge + const vtkIdType& c0 = pointIds[4]; + vtkIdType c1 = pointIds[1]; + const vtkIdType c2 = pointIds[1] - 1; + vtkIdType c3 = pointIds[4] - 1; + + T s0 = s0Origin; + T s1 = snet->GetVoxelForTriad(i, row + 1, slice); + if (s0 == backgroundLabel || (s1 != backgroundLabel && s0 > s1)) + { + std::swap(s0, s1); + std::swap(c1, c3); + } + + *connIter++ = c0; + *connIter++ = c1; + *connIter++ = c2; + *connIter++ = c3; - const T s1 = snet->GetVoxelForTriad(i, row + 1, slice); - snet->WriteScalarTuple(s0, s1, backgroundLabel, quadId++); + snet->WriteScalarTuple(s0, s1, quadId++); } if (SurfaceNets::GenerateYZQuad(triad)) { - *connIter++ = pointIds[4]; - *connIter++ = pointIds[3]; - *connIter++ = pointIds[0]; - *connIter++ = pointIds[1]; // quad is normal to the x edge + const vtkIdType& c0 = pointIds[4]; + vtkIdType c1 = pointIds[3]; + const vtkIdType c2 = pointIds[0]; + vtkIdType c3 = pointIds[1]; + + T s0 = s0Origin; + T s1 = snet->GetVoxelForTriad(i + 1, row, slice); + if (s0 == backgroundLabel || (s1 != backgroundLabel && s0 > s1)) + { + std::swap(s0, s1); + std::swap(c1, c3); + } + + *connIter++ = c0; + *connIter++ = c1; + *connIter++ = c2; + *connIter = c3; - const T s1 = snet->GetVoxelForTriad(i + 1, row, slice); - snet->WriteScalarTuple(s0, s1, backgroundLabel, quadId++); + snet->WriteScalarTuple(s0, s1, quadId++); } } // operator() @@ -882,17 +918,10 @@ struct SurfaceNets } // Helper function writes the scalar 2-tuple. - void WriteScalarTuple(T s0, T s1, T backgroundLabel, vtkIdType quadId) + void WriteScalarTuple(T s0, T s1, vtkIdType quadId) { T* scalars = this->NewScalars + 2 * quadId; - - if (s0 == backgroundLabel || (s1 != backgroundLabel && s0 > s1)) - { - // Background label is placed last; s0 Date: Thu, 31 Jul 2025 15:06:43 -0400 Subject: [PATCH 4/5] Add TestSurfaceNets3DNormalsConsistency --- Filters/Core/Testing/Cxx/CMakeLists.txt | 1 + .../TestSurfaceNets3DNormalsConsistency.cxx | 156 ++++++++++++++++++ 2 files changed, 157 insertions(+) create mode 100644 Filters/Core/Testing/Cxx/TestSurfaceNets3DNormalsConsistency.cxx diff --git a/Filters/Core/Testing/Cxx/CMakeLists.txt b/Filters/Core/Testing/Cxx/CMakeLists.txt index acaffa7324c9..9b6e22308b80 100644 --- a/Filters/Core/Testing/Cxx/CMakeLists.txt +++ b/Filters/Core/Testing/Cxx/CMakeLists.txt @@ -83,6 +83,7 @@ vtk_add_test_cxx(vtkFiltersCoreCxxTests tests TestStaticCleanPolyData.cxx,NO_VALID TestStripper.cxx,NO_VALID TestStructuredGridAppend.cxx,NO_VALID + TestSurfaceNets3DNormalsConsistency.cxx,NO_DATA,NO_VALID TestSynchronizedTemplates2D.cxx,NO_VALID TestSynchronizedTemplates2DRGB.cxx,NO_DATA,NO_VALID TestThreshold.cxx,NO_VALID diff --git a/Filters/Core/Testing/Cxx/TestSurfaceNets3DNormalsConsistency.cxx b/Filters/Core/Testing/Cxx/TestSurfaceNets3DNormalsConsistency.cxx new file mode 100644 index 000000000000..6d255b45abbc --- /dev/null +++ b/Filters/Core/Testing/Cxx/TestSurfaceNets3DNormalsConsistency.cxx @@ -0,0 +1,156 @@ +// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen +// SPDX-License-Identifier: BSD-3-Clause + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +int TestSurfaceNets3DNormalsConsistency(int, char*[]) +{ + // Create labeled image + vtkNew image; + image->SetDimensions(50, 50, 50); + image->AllocateScalars(VTK_INT, 1); + + // Fill with zero + std::fill_n(static_cast(image->GetScalarPointer()), 50 * 50 * 50, 0); + + // Add 4 adjacent cubes (20x20x20) with values 1-4 around the center edge (25,25,:) + for (int z = 15; z < 35; ++z) + { + for (int y = 15; y < 35; ++y) + { + for (int x = 0; x < 50; ++x) + { + int value = 0; + if (x >= 15 && x < 35) + { + if (x < 25 && y < 25) + { + value = 1; + } + else if (x >= 25 && y < 25) + { + value = 2; + } + else if (x < 25) + { + value = 3; + } + else + { + value = 4; + } + } + if (value > 0) + { + *static_cast(image->GetScalarPointer(x, y, z)) = value; + } + } + } + } + + // Run vtkSurfaceNets3D + vtkNew surfacenets; + surfacenets->SetBackgroundLabel(0); + surfacenets->SetInputData(image); + surfacenets->SetInputArrayToProcess( + 0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "ImageScalars"); + surfacenets->SetOutputMeshTypeToTriangles(); + surfacenets->Update(); + + vtkPolyData* surface = surfacenets->GetOutput(); + + // Run vtkPolyDataNormals + vtkNew normalsFilter_original; + normalsFilter_original->SetInputData(surface); + normalsFilter_original->ComputeCellNormalsOn(); + normalsFilter_original->ComputePointNormalsOff(); + normalsFilter_original->AutoOrientNormalsOff(); // Don't change orientation + normalsFilter_original->Update(); + + // Run vtkPolyDataNormals + vtkNew normalsFilter; + normalsFilter->SetInputData(surface); + normalsFilter->ComputeCellNormalsOn(); + normalsFilter->ComputePointNormalsOff(); + normalsFilter->AutoOrientNormalsOn(); // Ensure consistent normal orientation + normalsFilter->Update(); + + // Check normals are consistent + vtkDataArray* normalsBefore = normalsFilter_original->GetOutput()->GetCellData()->GetNormals(); + vtkDataArray* normalsAfter = normalsFilter->GetOutput()->GetCellData()->GetNormals(); + + if (!normalsBefore || !normalsAfter) + { + std::cerr << "Normals not present before or after vtkPolyDataNormals.\n"; + return EXIT_FAILURE; + } + + if (normalsBefore->GetNumberOfTuples() != normalsAfter->GetNumberOfTuples()) + { + std::cerr << "Number of normals mismatch.\n"; + return EXIT_FAILURE; + } + + for (vtkIdType i = 0; i < normalsBefore->GetNumberOfTuples(); ++i) + { + double n0[3], n1[3]; + normalsBefore->GetTuple(i, n0); + normalsAfter->GetTuple(i, n1); + double dot = vtkMath::Dot(n0, n1); + if (dot < 0.99) + { // Allow for small deviation + std::cerr << "Normal mismatch at point " << i << ": dot = " << dot << "\n"; + return EXIT_FAILURE; + } + } + + // Check BoundaryLabels CellData array + vtkIntArray* labels = + vtkIntArray::SafeDownCast(surface->GetCellData()->GetArray("BoundaryLabels")); + if (!labels) + { + std::cerr << "BoundaryLabels array missing from CellData.\n"; + return EXIT_FAILURE; + } + + if (labels->GetNumberOfComponents() != 2) + { + std::cerr << "BoundaryLabels does not have 2 components.\n"; + return EXIT_FAILURE; + } + + for (vtkIdType i = 0; i < labels->GetNumberOfTuples(); ++i) + { + int v[2]; + labels->GetTypedTuple(i, v); + if (v[0] == 0) + { + std::cerr << "BoundaryLabels at cell " << i << " background must be second component: [" + << v[0] << ", " << v[1] << "]\n"; + return EXIT_FAILURE; + } + if (v[0] >= v[1] && v[1] != 0) + { + std::cerr << "BoundaryLabels at cell " << i << " are not sorted: [" << v[0] << ", " << v[1] + << "]\n"; + return EXIT_FAILURE; + } + } + + std::cout << "Test passed: Surface normals consistent and BoundaryLabels correctly formatted.\n"; + return EXIT_SUCCESS; +} -- GitLab From e948892f82188655f9f482cfa57b3839b7ccea3d Mon Sep 17 00:00:00 2001 From: Spiros Tsalikis Date: Tue, 7 Oct 2025 12:49:09 -0400 Subject: [PATCH 5/5] Add changelog --- .../dev/surfacenets3d-fix-consistent-normal-orientation.md | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 Documentation/release/dev/surfacenets3d-fix-consistent-normal-orientation.md diff --git a/Documentation/release/dev/surfacenets3d-fix-consistent-normal-orientation.md b/Documentation/release/dev/surfacenets3d-fix-consistent-normal-orientation.md new file mode 100644 index 000000000000..71c3bce7875c --- /dev/null +++ b/Documentation/release/dev/surfacenets3d-fix-consistent-normal-orientation.md @@ -0,0 +1,3 @@ +## vtkSurfaceNets3D: Fix Consistent Normal Orientation + +`vtkSurfaceNets3D` now generates meshes with consistent normal orientations. -- GitLab