/*========================================================================= * * 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 * * https://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. * *=========================================================================*/ #ifndef itkBinaryImageToLevelSetImageAdaptor_hxx #define itkBinaryImageToLevelSetImageAdaptor_hxx #include "itkConnectedImageNeighborhoodShape.h" #include "itkSignedMaurerDistanceMapImageFilter.h" namespace itk { template BinaryImageToLevelSetImageAdaptor>::BinaryImageToLevelSetImageAdaptor() { this->m_SignedDistanceTransformFilter = SignedMaurerDistanceMapImageFilter::New(); } template BinaryImageToLevelSetImageAdaptor>::~BinaryImageToLevelSetImageAdaptor() = default; template void BinaryImageToLevelSetImageAdaptor>::Initialize() { if (this->m_InputImage.IsNull()) { itkGenericExceptionMacro("m_InputImage is nullptr"); } if (m_SignedDistanceTransformFilter.IsNull()) { itkGenericExceptionMacro("m_SignedDistanceTransformFilter is nullptr"); } m_SignedDistanceTransformFilter->SetInput(this->m_InputImage); m_SignedDistanceTransformFilter->Update(); auto tempImage = LevelSetImageType::New(); tempImage->Graft(m_SignedDistanceTransformFilter->GetOutput()); this->m_LevelSet = LevelSetType::New(); this->m_LevelSet->SetImage(tempImage); } template BinaryImageToLevelSetImageAdaptor>:: BinaryImageToLevelSetImageAdaptor() = default; template BinaryImageToLevelSetImageAdaptor>:: ~BinaryImageToLevelSetImageAdaptor() = default; template void BinaryImageToLevelSetImageAdaptor>::Initialize() { if (this->m_InputImage.IsNull()) { itkGenericExceptionMacro("m_InputImage is nullptr"); } this->m_LabelMap = LevelSetLabelMapType::New(); this->m_LabelMap->SetBackgroundValue(LevelSetType::PlusThreeLayer()); this->m_LabelMap->CopyInformation(this->m_InputImage); this->m_InternalImage = InternalImageType::New(); this->m_InternalImage->CopyInformation(this->m_InputImage); this->m_InternalImage->SetRegions(this->m_InputImage->GetBufferedRegion()); this->m_InternalImage->Allocate(); this->m_InternalImage->FillBuffer(LevelSetType::PlusThreeLayer()); LevelSetLabelObjectPointer innerPart = LevelSetLabelObjectType::New(); innerPart->SetLabel(LevelSetType::MinusThreeLayer()); // Precondition labelmap and phi InputIteratorType inputIt(this->m_InputImage, this->m_InputImage->GetLargestPossibleRegion()); inputIt.GoToBegin(); InternalIteratorType internalIt(this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); internalIt.GoToBegin(); while (!inputIt.IsAtEnd()) { if (inputIt.Get() != InputImagePixelType{}) { innerPart->AddIndex(inputIt.GetIndex()); internalIt.Set(LevelSetType::MinusThreeLayer()); } ++internalIt; ++inputIt; } innerPart->Optimize(); this->m_LabelMap->AddLabelObject(innerPart); FindActiveLayer(); FindPlusOneMinusOneLayer(); PropagateToOuterLayers(LevelSetType::MinusOneLayer(), LevelSetType::MinusTwoLayer(), LevelSetType::MinusThreeLayer()); PropagateToOuterLayers(LevelSetType::PlusOneLayer(), LevelSetType::PlusTwoLayer(), LevelSetType::PlusThreeLayer()); this->m_LabelMap->Optimize(); this->m_LevelSet->SetLabelMap(this->m_LabelMap); // release the memory this->m_InternalImage = nullptr; } template void BinaryImageToLevelSetImageAdaptor>:: PropagateToOuterLayers(LayerIdType layerToBeScanned, LayerIdType outputLayer, LayerIdType testValue) { const LevelSetLayerType layerPlus1 = this->m_LevelSet->GetLayer(layerToBeScanned); LevelSetLayerType & layerPlus2 = this->m_LevelSet->GetLayer(outputLayer); const auto plus2 = static_cast(outputLayer); typename NeighborhoodIteratorType::RadiusType radius; radius.Fill(1); ZeroFluxNeumannBoundaryCondition im_nbc; NeighborhoodIteratorType neighIt(radius, this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); neighIt.OverrideBoundaryCondition(&im_nbc); neighIt.ActivateOffsets(GenerateConnectedImageNeighborhoodShapeOffsets()); // iterate on the layer to be scanned auto nodeIt = layerPlus1.begin(); auto nodeEnd = layerPlus1.end(); while (nodeIt != nodeEnd) { LevelSetInputType idx = nodeIt->first; neighIt.SetLocation(idx); for (typename NeighborhoodIteratorType::Iterator it = neighIt.Begin(); !it.IsAtEnd(); ++it) { // check if in the neighborhood there are values equal to testValue if (it.Get() == testValue) { LevelSetInputType tempIndex = neighIt.GetIndex(it.GetNeighborhoodOffset()); layerPlus2.insert(LayerPairType(tempIndex, plus2)); } } ++nodeIt; } LevelSetLabelObjectPointer ObjectPlus2 = LevelSetLabelObjectType::New(); ObjectPlus2->SetLabel(static_cast(outputLayer)); nodeIt = layerPlus2.begin(); nodeEnd = layerPlus2.end(); while (nodeIt != nodeEnd) { if (testValue != this->m_LabelMap->GetBackgroundValue()) { this->m_LabelMap->GetLabelObject(testValue)->RemoveIndex(nodeIt->first); } ObjectPlus2->AddIndex(nodeIt->first); this->m_InternalImage->SetPixel(nodeIt->first, outputLayer); ++nodeIt; } ObjectPlus2->Optimize(); this->m_LabelMap->AddLabelObject(ObjectPlus2); } template void BinaryImageToLevelSetImageAdaptor>::FindActiveLayer() { LevelSetLabelObjectPointer labelObject = this->m_LabelMap->GetLabelObject(LevelSetType::MinusThreeLayer()); const LevelSetOutputType zero{}; LevelSetLayerType & layer0 = this->m_LevelSet->GetLayer(LevelSetType::ZeroLayer()); typename NeighborhoodIteratorType::RadiusType radius; radius.Fill(1); ZeroFluxNeumannBoundaryCondition im_nbc; NeighborhoodIteratorType neighIt(radius, this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); neighIt.OverrideBoundaryCondition(&im_nbc); neighIt.ActivateOffsets(GenerateConnectedImageNeighborhoodShapeOffsets()); typename LevelSetLabelObjectType::ConstIndexIterator lineIt(labelObject); lineIt.GoToBegin(); while (!lineIt.IsAtEnd()) { const LevelSetInputType & idx = lineIt.GetIndex(); neighIt.SetLocation(idx); for (typename NeighborhoodIteratorType::Iterator it = neighIt.Begin(); !it.IsAtEnd(); ++it) { if (it.Get() == LevelSetType::PlusThreeLayer()) { layer0.insert(LayerPairType(idx, zero)); break; } } ++lineIt; } if (!layer0.empty()) { LevelSetLabelObjectPointer ZeroSet = LevelSetLabelObjectType::New(); ZeroSet->SetLabel(LevelSetType::ZeroLayer()); LevelSetLayerConstIterator nodeIt = layer0.begin(); LevelSetLayerConstIterator nodeEnd = layer0.end(); while (nodeIt != nodeEnd) { this->m_LabelMap->GetLabelObject(LevelSetType::MinusThreeLayer())->RemoveIndex(nodeIt->first); ZeroSet->AddIndex(nodeIt->first); this->m_InternalImage->SetPixel(nodeIt->first, LevelSetType::ZeroLayer()); ++nodeIt; } ZeroSet->Optimize(); this->m_LabelMap->AddLabelObject(ZeroSet); } } template void BinaryImageToLevelSetImageAdaptor>:: FindPlusOneMinusOneLayer() { const LevelSetOutputType minus1 = -NumericTraits::OneValue(); const LevelSetOutputType plus1 = NumericTraits::OneValue(); const LevelSetLayerType layer0 = this->m_LevelSet->GetLayer(LevelSetType::ZeroLayer()); LevelSetLayerType & layerMinus1 = this->m_LevelSet->GetLayer(LevelSetType::MinusOneLayer()); LevelSetLayerType & layerPlus1 = this->m_LevelSet->GetLayer(LevelSetType::PlusOneLayer()); typename NeighborhoodIteratorType::RadiusType radius; radius.Fill(1); ZeroFluxNeumannBoundaryCondition im_nbc; NeighborhoodIteratorType neighIt(radius, this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); neighIt.OverrideBoundaryCondition(&im_nbc); neighIt.ActivateOffsets(GenerateConnectedImageNeighborhoodShapeOffsets()); auto nodeIt = layer0.begin(); auto nodeEnd = layer0.end(); while (nodeIt != nodeEnd) { LevelSetInputType idx = nodeIt->first; neighIt.SetLocation(idx); for (typename NeighborhoodIteratorType::Iterator it = neighIt.Begin(); !it.IsAtEnd(); ++it) { if (it.Get() == LevelSetType::PlusThreeLayer()) { LevelSetInputType tempIndex = neighIt.GetIndex(it.GetNeighborhoodOffset()); layerPlus1.insert(LayerPairType(tempIndex, plus1)); } if (it.Get() == LevelSetType::MinusThreeLayer()) { LevelSetInputType tempIndex = neighIt.GetIndex(it.GetNeighborhoodOffset()); layerMinus1.insert(LayerPairType(tempIndex, minus1)); } } ++nodeIt; } LevelSetLabelObjectPointer ObjectMinus1 = LevelSetLabelObjectType::New(); ObjectMinus1->SetLabel(LevelSetType::MinusOneLayer()); nodeIt = layerMinus1.begin(); nodeEnd = layerMinus1.end(); while (nodeIt != nodeEnd) { this->m_LabelMap->GetLabelObject(LevelSetType::MinusThreeLayer())->RemoveIndex(nodeIt->first); ObjectMinus1->AddIndex(nodeIt->first); this->m_InternalImage->SetPixel(nodeIt->first, LevelSetType::MinusOneLayer()); ++nodeIt; } ObjectMinus1->Optimize(); this->m_LabelMap->AddLabelObject(ObjectMinus1); LevelSetLabelObjectPointer ObjectPlus1 = LevelSetLabelObjectType::New(); ObjectPlus1->SetLabel(LevelSetType::PlusOneLayer()); nodeIt = layerPlus1.begin(); nodeEnd = layerPlus1.end(); while (nodeIt != nodeEnd) { ObjectPlus1->AddIndex(nodeIt->first); this->m_InternalImage->SetPixel(nodeIt->first, LevelSetType::PlusOneLayer()); ++nodeIt; } ObjectPlus1->Optimize(); this->m_LabelMap->AddLabelObject(ObjectPlus1); } //////////////////////////////////////////////////////////////////////////////// template BinaryImageToLevelSetImageAdaptor>::BinaryImageToLevelSetImageAdaptor() = default; template BinaryImageToLevelSetImageAdaptor>:: ~BinaryImageToLevelSetImageAdaptor() = default; template void BinaryImageToLevelSetImageAdaptor>::Initialize() { if (this->m_InputImage.IsNull()) { itkGenericExceptionMacro("m_InputImage is nullptr"); } this->m_LabelMap = LevelSetLabelMapType::New(); this->m_LabelMap->SetBackgroundValue(LevelSetType::PlusThreeLayer()); this->m_LabelMap->CopyInformation(this->m_InputImage); this->m_InternalImage = InternalImageType::New(); this->m_InternalImage->CopyInformation(this->m_InputImage); this->m_InternalImage->SetRegions(this->m_InputImage->GetBufferedRegion()); this->m_InternalImage->Allocate(); this->m_InternalImage->FillBuffer(LevelSetType::PlusThreeLayer()); LevelSetLabelObjectPointer innerPart = LevelSetLabelObjectType::New(); innerPart->SetLabel(LevelSetType::MinusThreeLayer()); // Precondition labelmap and phi InputIteratorType iIt(this->m_InputImage, this->m_InputImage->GetLargestPossibleRegion()); iIt.GoToBegin(); InternalIteratorType labelIt(this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); labelIt.GoToBegin(); while (!iIt.IsAtEnd()) { if (iIt.Get() != InputImagePixelType{}) { innerPart->AddIndex(iIt.GetIndex()); labelIt.Set(LevelSetType::MinusThreeLayer()); } ++labelIt; ++iIt; } innerPart->Optimize(); this->m_LabelMap->AddLabelObject(innerPart); FindActiveLayer(); this->m_LevelSet->SetLabelMap(this->m_LabelMap); this->m_InternalImage = nullptr; } template void BinaryImageToLevelSetImageAdaptor>::FindActiveLayer() { LevelSetLabelObjectPointer labelObject = this->m_LabelMap->GetLabelObject(LevelSetType::MinusThreeLayer()); LevelSetLayerType & layerMinus1 = this->m_LevelSet->GetLayer(LevelSetType::MinusOneLayer()); LevelSetLayerType & layerPlus1 = this->m_LevelSet->GetLayer(LevelSetType::PlusOneLayer()); typename NeighborhoodIteratorType::RadiusType radius; radius.Fill(1); ZeroFluxNeumannBoundaryCondition im_nbc; NeighborhoodIteratorType neighIt(radius, this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); neighIt.OverrideBoundaryCondition(&im_nbc); neighIt.ActivateOffsets(GenerateConnectedImageNeighborhoodShapeOffsets()); typename LevelSetLabelObjectType::ConstIndexIterator lineIt(labelObject); lineIt.GoToBegin(); while (!lineIt.IsAtEnd()) { const LevelSetInputType & idx = lineIt.GetIndex(); neighIt.SetLocation(idx); bool boundary = false; for (typename NeighborhoodIteratorType::Iterator it = neighIt.Begin(); !it.IsAtEnd(); ++it) { if (it.Get() == LevelSetType::PlusThreeLayer()) { LevelSetInputType tempIndex = neighIt.GetIndex(it.GetNeighborhoodOffset()); layerPlus1.insert(LayerPairType(tempIndex, LevelSetOutputType(1.0))); boundary = true; } } if (boundary) { layerMinus1.insert(LayerPairType(idx, LevelSetOutputType(-1.0))); } ++lineIt; } LevelSetLabelObjectPointer ObjectMinus1 = LevelSetLabelObjectType::New(); ObjectMinus1->SetLabel(LevelSetType::MinusOneLayer()); LevelSetLayerConstIterator nodeIt = layerMinus1.begin(); LevelSetLayerConstIterator nodeEnd = layerMinus1.end(); while (nodeIt != nodeEnd) { this->m_LabelMap->GetLabelObject(LevelSetType::MinusThreeLayer())->RemoveIndex(nodeIt->first); ObjectMinus1->AddIndex(nodeIt->first); ++nodeIt; } ObjectMinus1->Optimize(); this->m_LabelMap->AddLabelObject(ObjectMinus1); LevelSetLabelObjectPointer ObjectPlus1 = LevelSetLabelObjectType::New(); ObjectPlus1->SetLabel(LevelSetType::PlusOneLayer()); nodeIt = layerPlus1.begin(); nodeEnd = layerPlus1.end(); while (nodeIt != nodeEnd) { ObjectPlus1->AddIndex(nodeIt->first); ++nodeIt; } ObjectPlus1->Optimize(); this->m_LabelMap->AddLabelObject(ObjectPlus1); } //////////////////////////////////////////////////////////////////////////////// template BinaryImageToLevelSetImageAdaptor>:: BinaryImageToLevelSetImageAdaptor() : Superclass() {} template BinaryImageToLevelSetImageAdaptor>:: ~BinaryImageToLevelSetImageAdaptor() = default; template void BinaryImageToLevelSetImageAdaptor>::Initialize() { if (this->m_InputImage.IsNull()) { itkGenericExceptionMacro("m_InputImage is nullptr"); } this->m_LabelMap = LevelSetLabelMapType::New(); this->m_LabelMap->SetBackgroundValue(LevelSetType::PlusOneLayer()); this->m_LabelMap->CopyInformation(this->m_InputImage); this->m_InternalImage = InternalImageType::New(); this->m_InternalImage->CopyInformation(this->m_InputImage); this->m_InternalImage->SetRegions(this->m_InputImage->GetBufferedRegion()); this->m_InternalImage->Allocate(); this->m_InternalImage->FillBuffer(LevelSetType::PlusOneLayer()); LevelSetLabelObjectPointer innerPart = LevelSetLabelObjectType::New(); innerPart->SetLabel(LevelSetType::MinusOneLayer()); // Precondition labelmap and phi InputIteratorType inputIt(this->m_InputImage, this->m_InputImage->GetLargestPossibleRegion()); inputIt.GoToBegin(); InternalIteratorType internalIt(this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); internalIt.GoToBegin(); while (!inputIt.IsAtEnd()) { if (inputIt.Get() != InputImagePixelType{}) { innerPart->AddIndex(inputIt.GetIndex()); internalIt.Set(LevelSetType::MinusOneLayer()); } ++internalIt; ++inputIt; } innerPart->Optimize(); this->m_LabelMap->AddLabelObject(innerPart); this->FindActiveLayer(); this->CreateMinimalInterface(); this->m_LevelSet->SetLabelMap(this->m_LabelMap); this->m_InternalImage = nullptr; } template void BinaryImageToLevelSetImageAdaptor>::FindActiveLayer() { LevelSetLabelObjectPointer labelObject = this->m_LabelMap->GetLabelObject(LevelSetType::MinusOneLayer()); LevelSetLayerType & layer = this->m_LevelSet->GetLayer(LevelSetType::ZeroLayer()); typename NeighborhoodIteratorType::RadiusType radius; radius.Fill(1); ZeroFluxNeumannBoundaryCondition im_nbc; NeighborhoodIteratorType neighIt(radius, this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); neighIt.OverrideBoundaryCondition(&im_nbc); neighIt.ActivateOffsets(GenerateConnectedImageNeighborhoodShapeOffsets()); typename LevelSetLabelObjectType::ConstIndexIterator lineIt(labelObject); lineIt.GoToBegin(); while (!lineIt.IsAtEnd()) { const LevelSetInputType & idx = lineIt.GetIndex(); neighIt.SetLocation(idx); bool ZeroSet = false; for (typename NeighborhoodIteratorType::Iterator it = neighIt.Begin(); !it.IsAtEnd(); ++it) { if (it.Get() == LevelSetType::PlusOneLayer()) { ZeroSet = true; } } if (ZeroSet) { layer.insert(LayerPairType(idx, LevelSetOutputType{})); this->m_InternalImage->SetPixel(idx, LevelSetType::ZeroLayer()); } ++lineIt; } LevelSetLabelObjectPointer ObjectZero = LevelSetLabelObjectType::New(); ObjectZero->SetLabel(LevelSetType::ZeroLayer()); auto nodeIt = layer.begin(); auto nodeEnd = layer.end(); while (nodeIt != nodeEnd) { labelObject->RemoveIndex(nodeIt->first); ObjectZero->AddIndex(nodeIt->first); ++nodeIt; } ObjectZero->Optimize(); this->m_LabelMap->AddLabelObject(ObjectZero); } template void BinaryImageToLevelSetImageAdaptor>::CreateMinimalInterface() { LevelSetLayerType & list_0 = this->m_LevelSet->GetLayer(LevelSetType::ZeroLayer()); ZeroFluxNeumannBoundaryCondition sp_nbc; typename NeighborhoodIteratorType::RadiusType radius; radius.Fill(1); NeighborhoodIteratorType neighIt(radius, this->m_InternalImage, this->m_InternalImage->GetLargestPossibleRegion()); neighIt.OverrideBoundaryCondition(&sp_nbc); neighIt.ActivateOffsets(GenerateConnectedImageNeighborhoodShapeOffsets()); auto nodeIt = list_0.begin(); auto nodeEnd = list_0.end(); while (nodeIt != nodeEnd) { LevelSetInputType currentIdx = nodeIt->first; neighIt.SetLocation(currentIdx); bool hasPositiveLayerNeighbor = false; bool hasNegativeLayerNeighbor = false; for (typename NeighborhoodIteratorType::Iterator i = neighIt.Begin(); !i.IsAtEnd(); ++i) { LayerIdType tempValue = i.Get(); if (tempValue != LayerIdType{}) { if (tempValue == LevelSetType::MinusOneLayer()) { hasNegativeLayerNeighbor = true; if (hasPositiveLayerNeighbor) { break; } } else // ( tempValue == NumericTraits< LevelSetOutputType >::OneValue() ) { hasPositiveLayerNeighbor = true; if (hasNegativeLayerNeighbor) { break; } } } } if (hasNegativeLayerNeighbor && !hasPositiveLayerNeighbor) { auto tempIt = nodeIt; ++nodeIt; list_0.erase(tempIt); this->m_LabelMap->GetLabelObject(LevelSetType::ZeroLayer())->RemoveIndex(currentIdx); this->m_LabelMap->GetLabelObject(LevelSetType::MinusOneLayer())->AddIndex(currentIdx); } else { if (hasPositiveLayerNeighbor && !hasNegativeLayerNeighbor) { auto tempIt = nodeIt; ++nodeIt; list_0.erase(tempIt); this->m_LabelMap->GetLabelObject(LevelSetType::ZeroLayer())->RemoveIndex(currentIdx); } else { ++nodeIt; } } } } } // namespace itk #endif // itkBinaryImageToLevelSetImageAdaptor_hxx