Skip to content

Commit

Permalink
BUG: Handle Non-Manifold Geometry
Browse files Browse the repository at this point in the history
itk::QuadEdgeMesh enforces manifold geometry, and therefore cases
where vertices appear at the same physical location must be handled
carefully.  A bugfix was recently applied which handles the vast
majority (252/256) of cases.  However, in the method fails in the
four possible cases where two background pixels appear at opposite
corners of a 2x2x2 region, with foreground pixels elsewhere.
This patch handles these four cases manually.  Testing and
documentation are updated as appropriate.
  • Loading branch information
sudomakeinstall committed Jun 2, 2020
1 parent 978cdbd commit 5cdccb1
Show file tree
Hide file tree
Showing 3 changed files with 132 additions and 20 deletions.
51 changes: 31 additions & 20 deletions include/itkCuberilleImageToMeshFilter.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -682,15 +682,22 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
{

// The commented code below iterates through all possible binary 2x2x2 regions,
// runs connected components on the result. The total number of foreground
// and runs connected components on the result. The total number of foreground
// connected components is equal to the number of vertices which must be
// replicated--these vertices are stored in the VertexMap. One less than the
// number of the connected component is equal to the index of the corresponding
// vertex in the VertexMap. While it would be easily possible to generate this at
// runtime, it could theoretically cause a significant performance hit in a
// case where many instances of this class must be instantiated. Therefore,
// the values are hardcoded below. The commented code was used to generate
// the hardcoded values.
// vertex in the VertexMap.
//
// This method works for all but four of the 256 possible cases. The failing
// cases occur when there are two background pixels at opposite corners of the
// 2x2x2 region, and all remaining pixels are foreground. In these cases,
// there is only one connected component, but two vertices must be added. These
// cases are handled separately.
//
// While it would be possible to generate this at runtime, it could theoretically
// cause a significant performance hit in a case where many instances of this class
// must be instantiated. Therefore, the values are hardcoded below. The commented
// code was used to generate the hardcoded values.
//
// NOTE: There was previously a bug in itk::ConnectedComponentImageFilter which
// caused the hardcoded values below to be incorrect, though this has since been
Expand All @@ -702,7 +709,7 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
// using TImage = itk::Image<unsigned char, 3>;
// using TConnected = itk::ConnectedComponentImageFilter< TImage, TImage >;
//
// for (size_t mask = 0; mask < pow(2, 8); ++mask) {
// for (size_t mask = 0; mask < std::pow(2, 8); ++mask) {
//
// std::bitset<8> bitmask(mask);
//
Expand All @@ -719,7 +726,7 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
// image->Allocate();
// image->FillBuffer( 0 );
//
// for (size_t index = 0; index < pow(2, 3); ++index) {
// for (size_t index = 0; index < std::pow(2, 3); ++index) {
// std::bitset<3> bitindex(index);
// image->SetPixel( {{bitindex[0], bitindex[1], bitindex[2]}}, bitmask[index] );
// }
Expand All @@ -731,7 +738,7 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
//
// TLabels labels;
//
// for (size_t index = 0; index < pow(2, 3); ++index) {
// for (size_t index = 0; index < std::pow(2, 3); ++index) {
// std::bitset<3> bitindex(index);
// const auto component = connected->GetOutput()->GetPixel({{
// bitindex[0],
Expand All @@ -745,18 +752,22 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
//
// }
//
// for (size_t i = 0; i < pow(2,8); ++i) {
// // Manually handle the corner cases, discussed above.
//
// this->m_LabelsArray[126] = { -1, 0, 0, 1, 0, 1, 1, -1 };
// this->m_LabelsArray[189] = { 0, -1, 1, 0, 1, 0, -1, 1 };
// this->m_LabelsArray[219] = { 0, 1, -1, 0, 1, -1, 0, 1 };
// this->m_LabelsArray[231] = { 1, 0, 0, -1, -1, 1, 1, 0 };
//
// for (size_t i = 0; i < std::pow(2,8); ++i) {
// std::cout << "this->m_LabelsArray[" << i << "] ";
// if (i < 10) std::cout << ' ';
// if (i < 100) std::cout << ' ';
// std::cout << "= {";
// for (size_t j = 0; j < pow(2,3); ++j) {
// for (size_t j = 0; j < std::pow(2,3); ++j) {
// const auto k = static_cast<int>(this->m_LabelsArray[i][j]);
// if (k >= 0) std::cout << ' ';
// std::cout << k;
// std::cout << ' ' << k;
// if (j != 7) std::cout << ',';
// }
// std::cout << "};\n";
// std::cout << " };\n";
// }

this->m_LabelsArray[0] = { -1, -1, -1, -1, -1, -1, -1, -1 };
Expand Down Expand Up @@ -885,7 +896,7 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
this->m_LabelsArray[123] = { 0, 0, -1, 0, 0, 0, 0, -1 };
this->m_LabelsArray[124] = { -1, -1, 0, 0, 0, 0, 0, -1 };
this->m_LabelsArray[125] = { 0, -1, 0, 0, 0, 0, 0, -1 };
this->m_LabelsArray[126] = { -1, 0, 0, 0, 0, 0, 0, -1 };
this->m_LabelsArray[126] = { -1, 0, 0, 1, 0, 1, 1, -1 };
this->m_LabelsArray[127] = { 0, 0, 0, 0, 0, 0, 0, -1 };
this->m_LabelsArray[128] = { -1, -1, -1, -1, -1, -1, -1, 0 };
this->m_LabelsArray[129] = { 0, -1, -1, -1, -1, -1, -1, 1 };
Expand Down Expand Up @@ -948,7 +959,7 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
this->m_LabelsArray[186] = { -1, 0, -1, 0, 0, 0, -1, 0 };
this->m_LabelsArray[187] = { 0, 0, -1, 0, 0, 0, -1, 0 };
this->m_LabelsArray[188] = { -1, -1, 0, 0, 0, 0, -1, 0 };
this->m_LabelsArray[189] = { 0, -1, 0, 0, 0, 0, -1, 0 };
this->m_LabelsArray[189] = { 0, -1, 1, 0, 1, 0, -1, 1 };
this->m_LabelsArray[190] = { -1, 0, 0, 0, 0, 0, -1, 0 };
this->m_LabelsArray[191] = { 0, 0, 0, 0, 0, 0, -1, 0 };
this->m_LabelsArray[192] = { -1, -1, -1, -1, -1, -1, 0, 0 };
Expand Down Expand Up @@ -978,7 +989,7 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
this->m_LabelsArray[216] = { -1, -1, -1, 0, 0, -1, 0, 0 };
this->m_LabelsArray[217] = { 0, -1, -1, 0, 0, -1, 0, 0 };
this->m_LabelsArray[218] = { -1, 0, -1, 0, 0, -1, 0, 0 };
this->m_LabelsArray[219] = { 0, 0, -1, 0, 0, -1, 0, 0 };
this->m_LabelsArray[219] = { 0, 1, -1, 0, 1, -1, 0, 1 };
this->m_LabelsArray[220] = { -1, -1, 0, 0, 0, -1, 0, 0 };
this->m_LabelsArray[221] = { 0, -1, 0, 0, 0, -1, 0, 0 };
this->m_LabelsArray[222] = { -1, 0, 0, 0, 0, -1, 0, 0 };
Expand All @@ -990,7 +1001,7 @@ CuberilleImageToMeshFilter<TInputImage, TOutputMesh, TInterpolator>::CalculateLa
this->m_LabelsArray[228] = { -1, -1, 0, -1, -1, 0, 0, 0 };
this->m_LabelsArray[229] = { 0, -1, 0, -1, -1, 0, 0, 0 };
this->m_LabelsArray[230] = { -1, 0, 0, -1, -1, 0, 0, 0 };
this->m_LabelsArray[231] = { 0, 0, 0, -1, -1, 0, 0, 0 };
this->m_LabelsArray[231] = { 1, 0, 0, -1, -1, 1, 1, 0 };
this->m_LabelsArray[232] = { -1, -1, -1, 0, -1, 0, 0, 0 };
this->m_LabelsArray[233] = { 0, -1, -1, 1, -1, 1, 1, 1 };
this->m_LabelsArray[234] = { -1, 0, -1, 0, -1, 0, 0, 0 };
Expand Down
4 changes: 4 additions & 0 deletions test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ set(CuberilleTests
CuberilleTest01.cxx
CuberilleTest02.cxx
CuberilleTest03.cxx
CuberilleTest04.cxx
)

CreateTestDriver(Cuberille "${Cuberille-Test_LIBRARIES}" "${CuberilleTests}")
Expand Down Expand Up @@ -355,3 +356,6 @@ itk_add_test(NAME CuberilleTestDirectionMatrix

itk_add_test(NAME CuberilleTestCellData
COMMAND ${itk-module}TestDriver CuberilleTest03 )

itk_add_test(NAME CuberilleTestNonManifoldGeometry
COMMAND ${itk-module}TestDriver CuberilleTest04 )
97 changes: 97 additions & 0 deletions test/CuberilleTest04.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/

#include <itkTestingMacros.h>
#include <itkImage.h>
#include <itkConstantPadImageFilter.h>
#include <itkQuadEdgeMesh.h>
#include <itkCuberilleImageToMeshFilter.h>

// In the case where a 2x2x2 region contains two background pixels at
// opposite corners, and foreground pixels elsewhere, it is necessary
// that there be two vertices at the center. Previously, this was
// handled incorrectly with a single vertex, causing non-manifold
// geometry. This test covers the four cases in which this occurs.

int
CuberilleTest04(int itkNotUsed(argc), char * itkNotUsed(argv)[])
{

using TPixel = unsigned char;
using TCoordinate = double;
const unsigned int Dimension = 3;
using TImage = itk::Image<TPixel, Dimension>;
using TMesh = itk::QuadEdgeMesh<TCoordinate, Dimension>;
using TPad = itk::ConstantPadImageFilter<TImage, TImage>;
using TExtract = itk::CuberilleImageToMeshFilter<TImage, TMesh>;

// 126 01111110
// 189 10111101
// 219 11011011
// 231 11100111

std::array<size_t, 4> masks{ 126, 189, 219, 231 };

for (const auto & mask : masks)
{

std::bitset<8> bitmask(mask);

std::cout << mask << ' ' << bitmask << std::endl;

const auto image = TImage::New();
TImage::SizeType size;
size.Fill(2);

TImage::IndexType origin;
origin.Fill(0);

TImage::RegionType region(origin, size);

image->SetRegions(region);
image->Allocate();
image->FillBuffer(0);

for (size_t index = 0; index < std::pow(2, 3); ++index)
{
std::bitset<3> bitindex(index);
image->SetPixel({ { bitindex[0], bitindex[1], bitindex[2] } }, bitmask[index]);
}

typename TImage::SizeType padding;
padding.Fill(1);

const auto pad = TPad::New();
pad->SetInput(image);
pad->SetPadUpperBound(padding);
pad->SetPadLowerBound(padding);
pad->SetConstant(static_cast<TPixel>(0));

const auto extract = TExtract::New();
extract->SetInput(pad->GetOutput());
extract->GenerateTriangleFacesOn();
extract->ProjectVerticesToIsoSurfaceOff();
extract->SavePixelAsCellDataOn();
extract->Update();

ITK_TEST_EXPECT_EQUAL(extract->GetOutput()->GetNumberOfPoints(), 26);
ITK_TEST_EXPECT_EQUAL(extract->GetOutput()->GetNumberOfCells(), 48);
}

return EXIT_SUCCESS;
}

0 comments on commit 5cdccb1

Please sign in to comment.