diff --git a/Modules/DICOM/src/mitkGantryTiltInformation.cpp b/Modules/DICOM/src/mitkGantryTiltInformation.cpp index b6efd46967..f268ddf8c8 100644 --- a/Modules/DICOM/src/mitkGantryTiltInformation.cpp +++ b/Modules/DICOM/src/mitkGantryTiltInformation.cpp @@ -1,254 +1,255 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "mitkGantryTiltInformation.h" #include "mitkDICOMTag.h" #include "mitkLogMacros.h" mitk::GantryTiltInformation::GantryTiltInformation() : m_ShiftUp(0.0) , m_ShiftRight(0.0) , m_ShiftNormal(0.0) , m_ITKAssumedSliceSpacing(0.0) , m_NumberOfSlicesApart(0) { } #define doublepoint(x) \ Point3Dd x; \ x[0] = x ## f[0]; \ x[1] = x ## f[1]; \ x[2] = x ## f[2]; #define doublevector(x) \ Vector3Dd x; \ x[0] = x ## f[0]; \ x[1] = x ## f[1]; \ x[2] = x ## f[2]; mitk::GantryTiltInformation::GantryTiltInformation( const Point3D& origin1f, const Point3D& origin2f, const Vector3D& rightf, const Vector3D& upf, unsigned int numberOfSlicesApart) : m_ShiftUp(0.0) , m_ShiftRight(0.0) , m_ShiftNormal(0.0) +, m_ITKAssumedSliceSpacing(0.0) , m_NumberOfSlicesApart(numberOfSlicesApart) { assert(numberOfSlicesApart); doublepoint(origin1); doublepoint(origin2); doublevector(right); doublevector(up); // determine if slice 1 (imagePosition1 and imageOrientation1) and slice 2 can be in one orthogonal slice stack: // calculate a line from origin 1, directed along the normal of slice (calculated as the cross product of orientation 1) // check if this line passes through origin 2 /* Determine if line (imagePosition2 + l * normal) contains imagePosition1. Done by calculating the distance of imagePosition1 from line (imagePosition2 + l *normal) E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2 / |pointAlongNormal - origin2| ^ 2 ( x meaning the cross product ) */ Vector3Dd normal = itk::CrossProduct(right, up); Point3Dd pointAlongNormal = origin2 + normal; double numerator = itk::CrossProduct( pointAlongNormal - origin2 , origin2 - origin1 ).GetSquaredNorm(); double denominator = (pointAlongNormal - origin2).GetSquaredNorm(); double distance = sqrt(numerator / denominator); if ( distance > 0.001 ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt { MITK_DEBUG << " Series seems to contain a tilted (or sheared) geometry"; MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance; MITK_DEBUG << " ==> storing this shift for later analysis:"; MITK_DEBUG << " v origin1: " << origin1; MITK_DEBUG << " v origin2: " << origin2; MITK_DEBUG << " v right: " << right; MITK_DEBUG << " v up: " << up; MITK_DEBUG << " v normal: " << normal; Point3Dd projectionRight = projectPointOnLine( origin1, origin2, right ); Point3Dd projectionNormal = projectPointOnLine( origin1, origin2, normal ); m_ShiftRight = (projectionRight - origin2).GetNorm(); m_ShiftNormal = (projectionNormal - origin2).GetNorm(); /* now also check to which side the image is shifted. Calculation e.g. from http://mathworld.wolfram.com/Point-PlaneDistance.html */ Point3Dd testPoint = origin1; Vector3Dd planeNormal = up; double signedDistance = ( planeNormal[0] * testPoint[0] + planeNormal[1] * testPoint[1] + planeNormal[2] * testPoint[2] - ( planeNormal[0] * origin2[0] + planeNormal[1] * origin2[1] + planeNormal[2] * origin2[2] ) ) / sqrt( planeNormal[0] * planeNormal[0] + planeNormal[1] * planeNormal[1] + planeNormal[2] * planeNormal[2] ); m_ShiftUp = signedDistance; m_ITKAssumedSliceSpacing = (origin2 - origin1).GetNorm(); // How do we now this is assumed? See header documentation for ITK code references //double itkAssumedSliceSpacing = sqrt( m_ShiftUp * m_ShiftUp + m_ShiftNormal * m_ShiftNormal ); MITK_DEBUG << " calculated from slices " << m_NumberOfSlicesApart << " slices apart"; MITK_DEBUG << " shift normal: " << m_ShiftNormal; MITK_DEBUG << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing; MITK_DEBUG << " shift up: " << m_ShiftUp; MITK_DEBUG << " shift right: " << m_ShiftRight; MITK_DEBUG << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; } } mitk::GantryTiltInformation mitk::GantryTiltInformation ::MakeFromTagValues( const std::string& origin1String, const std::string& origin2String, const std::string& orientationString, unsigned int numberOfSlicesApart) { Vector3D right; right.Fill(0.0); Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point bool orientationConversion(false); DICOMStringToOrientationVectors( orientationString, right, up, orientationConversion ); if (orientationConversion && !origin1String.empty() && !origin2String.empty() ) { bool firstOriginConversion(false); bool lastOriginConversion(false); Point3D firstOrigin = DICOMStringToPoint3D( origin1String, firstOriginConversion ); Point3D lastOrigin = DICOMStringToPoint3D( origin2String, lastOriginConversion ); if (firstOriginConversion && lastOriginConversion) { return GantryTiltInformation( firstOrigin, lastOrigin, right, up, numberOfSlicesApart ); } } std::stringstream ss; ss << "Invalid tag values when constructing tilt information from origin1 '" << origin1String << "', origin2 '" << origin2String << "', and orientation '" << orientationString << "'"; throw std::invalid_argument(ss.str()); } void mitk::GantryTiltInformation ::Print(std::ostream& os) const { os << " calculated from slices " << m_NumberOfSlicesApart << " slices apart" << std::endl; os << " shift normal: " << m_ShiftNormal << std::endl; os << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing << std::endl; os << " shift up: " << m_ShiftUp << std::endl; os << " shift right: " << m_ShiftRight << std::endl; os << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535 << std::endl; } mitk::Point3D mitk::GantryTiltInformation::projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ) { /** See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) */ Vector3Dd lineOriginToP = p - lineOrigin; double innerProduct = lineOriginToP * lineDirection; double factor = innerProduct / lineDirection.GetSquaredNorm(); Point3Dd projection = lineOrigin + factor * lineDirection; return projection; } double mitk::GantryTiltInformation::GetTiltCorrectedAdditionalSize(unsigned int imageSizeZ) const { return fabs(m_ShiftUp / static_cast(m_NumberOfSlicesApart) * static_cast(imageSizeZ-1)); } double mitk::GantryTiltInformation::GetTiltAngleInDegrees() const { return atan( fabs(m_ShiftUp) / m_ShiftNormal ) * 180.0 / 3.1415926535; } double mitk::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const { // so many mm need to be shifted per slice! return m_ShiftUp / static_cast(m_NumberOfSlicesApart); } double mitk::GantryTiltInformation::GetRealZSpacing() const { return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); } bool mitk::GantryTiltInformation::IsSheared() const { return m_NumberOfSlicesApart && ( fabs(m_ShiftRight) > 0.001 || fabs(m_ShiftUp) > 0.001); } bool mitk::GantryTiltInformation::IsRegularGantryTilt() const { return m_NumberOfSlicesApart && ( fabs(m_ShiftRight) < 0.001 && fabs(m_ShiftUp) > 0.001); } diff --git a/Modules/Multilabel/mitkLabelSetImage.cpp b/Modules/Multilabel/mitkLabelSetImage.cpp index 8ebc7838fb..9f3c76e61d 100644 --- a/Modules/Multilabel/mitkLabelSetImage.cpp +++ b/Modules/Multilabel/mitkLabelSetImage.cpp @@ -1,1012 +1,1022 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkLabelSetImage.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkImagePixelReadAccessor.h" #include "mitkImagePixelWriteAccessor.h" #include "mitkInteractionConst.h" #include "mitkLookupTableProperty.h" #include "mitkPadImageFilter.h" #include "mitkRenderingManager.h" #include "mitkDICOMSegmentationPropertyHelper.h" #include "mitkDICOMQIPropertyHelper.h" #include #include #include #include #include #include //#include #include template void SetToZero(itk::Image *source) { source->FillBuffer(0); } template void CreateLabelMaskProcessing(mitk::Image *layerImage, mitk::Image *mask, mitk::LabelSet::PixelType index) { mitk::ImagePixelReadAccessor readAccessor(layerImage); mitk::ImagePixelWriteAccessor writeAccessor(mask); std::size_t numberOfPixels = 1; for (int dim = 0; dim < static_cast(VImageDimension); ++dim) numberOfPixels *= static_cast(readAccessor.GetDimension(dim)); auto src = readAccessor.GetData(); auto dest = writeAccessor.GetData(); for (std::size_t i = 0; i < numberOfPixels; ++i) { if (index == *(src + i)) *(dest + i) = 1; } } mitk::LabelSetImage::LabelSetImage() : mitk::Image(), m_ActiveLayer(0), m_activeLayerInvalid(false), m_ExteriorLabel(nullptr) { // Iniitlaize Background Label mitk::Color color; color.Set(0, 0, 0); m_ExteriorLabel = mitk::Label::New(); m_ExteriorLabel->SetColor(color); m_ExteriorLabel->SetName("Exterior"); m_ExteriorLabel->SetOpacity(0.0); m_ExteriorLabel->SetLocked(false); m_ExteriorLabel->SetValue(0); // Add some DICOM Tags as properties to segmentation image DICOMSegmentationPropertyHelper::DeriveDICOMSegmentationProperties(this); } mitk::LabelSetImage::LabelSetImage(const mitk::LabelSetImage &other) : Image(other), m_ActiveLayer(other.GetActiveLayer()), m_activeLayerInvalid(false), m_ExteriorLabel(other.GetExteriorLabel()->Clone()) { for (unsigned int i = 0; i < other.GetNumberOfLayers(); i++) { // Clone LabelSet data mitk::LabelSet::Pointer lsClone = other.GetLabelSet(i)->Clone(); // add modified event listener to LabelSet (listen to LabelSet changes) itk::SimpleMemberCommand::Pointer command = itk::SimpleMemberCommand::New(); command->SetCallbackFunction(this, &mitk::LabelSetImage::OnLabelSetModified); lsClone->AddObserver(itk::ModifiedEvent(), command); m_LabelSetContainer.push_back(lsClone); // clone layer Image data mitk::Image::Pointer liClone = other.GetLayerImage(i)->Clone(); m_LayerContainer.push_back(liClone); } // Add some DICOM Tags as properties to segmentation image DICOMSegmentationPropertyHelper::DeriveDICOMSegmentationProperties(this); } void mitk::LabelSetImage::OnLabelSetModified() { Superclass::Modified(); } void mitk::LabelSetImage::SetExteriorLabel(mitk::Label *label) { m_ExteriorLabel = label; } mitk::Label *mitk::LabelSetImage::GetExteriorLabel() { return m_ExteriorLabel; } const mitk::Label *mitk::LabelSetImage::GetExteriorLabel() const { return m_ExteriorLabel; } void mitk::LabelSetImage::Initialize(const mitk::Image *other) { mitk::PixelType pixelType(mitk::MakeScalarPixelType()); if (other->GetDimension() == 2) { const unsigned int dimensions[] = {other->GetDimension(0), other->GetDimension(1), 1}; Superclass::Initialize(pixelType, 3, dimensions); } else { Superclass::Initialize(pixelType, other->GetDimension(), other->GetDimensions()); } auto originalGeometry = other->GetTimeGeometry()->Clone(); this->SetTimeGeometry(originalGeometry); // initialize image memory to zero if (4 == this->GetDimension()) { AccessFixedDimensionByItk(this, SetToZero, 4); } else { AccessByItk(this, SetToZero); } // Transfer some general DICOM properties from the source image to derived image (e.g. Patient information,...) DICOMQIPropertyHelper::DeriveDICOMSourceProperties(other, this); // Add a inital LabelSet ans corresponding image data to the stack AddLayer(); } mitk::LabelSetImage::~LabelSetImage() { m_LabelSetContainer.clear(); } mitk::Image *mitk::LabelSetImage::GetLayerImage(unsigned int layer) { return m_LayerContainer[layer]; } const mitk::Image *mitk::LabelSetImage::GetLayerImage(unsigned int layer) const { return m_LayerContainer[layer]; } unsigned int mitk::LabelSetImage::GetActiveLayer() const { return m_ActiveLayer; } unsigned int mitk::LabelSetImage::GetNumberOfLayers() const { return m_LabelSetContainer.size(); } void mitk::LabelSetImage::RemoveLayer() { int layerToDelete = GetActiveLayer(); // remove all observers from active label set GetLabelSet(layerToDelete)->RemoveAllObservers(); // set the active layer to one below, if exists. if (layerToDelete != 0) { SetActiveLayer(layerToDelete - 1); } else { // we are deleting layer zero, it should not be copied back into the vector m_activeLayerInvalid = true; } // remove labelset and image data m_LabelSetContainer.erase(m_LabelSetContainer.begin() + layerToDelete); m_LayerContainer.erase(m_LayerContainer.begin() + layerToDelete); if (layerToDelete == 0) { this->SetActiveLayer(layerToDelete); } this->Modified(); } unsigned int mitk::LabelSetImage::AddLayer(mitk::LabelSet::Pointer lset) { mitk::Image::Pointer newImage = mitk::Image::New(); newImage->Initialize(this->GetPixelType(), this->GetDimension(), this->GetDimensions(), this->GetImageDescriptor()->GetNumberOfChannels()); newImage->SetTimeGeometry(this->GetTimeGeometry()->Clone()); if (newImage->GetDimension() < 4) { AccessByItk(newImage, SetToZero); } else { AccessFixedDimensionByItk(newImage, SetToZero, 4); } unsigned int newLabelSetId = this->AddLayer(newImage, lset); return newLabelSetId; } unsigned int mitk::LabelSetImage::AddLayer(mitk::Image::Pointer layerImage, mitk::LabelSet::Pointer lset) { unsigned int newLabelSetId = m_LayerContainer.size(); // Add labelset to layer mitk::LabelSet::Pointer ls; if (lset.IsNotNull()) { ls = lset; } else { ls = mitk::LabelSet::New(); ls->AddLabel(GetExteriorLabel()); ls->SetActiveLabel(0 /*Exterior Label*/); } ls->SetLayer(newLabelSetId); // Add exterior Label to label set // mitk::Label::Pointer exteriorLabel = CreateExteriorLabel(); // push a new working image for the new layer m_LayerContainer.push_back(layerImage); // push a new labelset for the new layer m_LabelSetContainer.push_back(ls); // add modified event listener to LabelSet (listen to LabelSet changes) itk::SimpleMemberCommand::Pointer command = itk::SimpleMemberCommand::New(); command->SetCallbackFunction(this, &mitk::LabelSetImage::OnLabelSetModified); ls->AddObserver(itk::ModifiedEvent(), command); SetActiveLayer(newLabelSetId); // MITK_INFO << GetActiveLayer(); this->Modified(); return newLabelSetId; } void mitk::LabelSetImage::AddLabelSetToLayer(const unsigned int layerIdx, const mitk::LabelSet::Pointer labelSet) { if (m_LayerContainer.size() <= layerIdx) { mitkThrow() << "Trying to add labelSet to non-existing layer."; } if (layerIdx < m_LabelSetContainer.size()) { m_LabelSetContainer[layerIdx] = labelSet; } else { while (layerIdx >= m_LabelSetContainer.size()) { mitk::LabelSet::Pointer defaultLabelSet = mitk::LabelSet::New(); defaultLabelSet->AddLabel(GetExteriorLabel()); defaultLabelSet->SetActiveLabel(0 /*Exterior Label*/); defaultLabelSet->SetLayer(m_LabelSetContainer.size()); m_LabelSetContainer.push_back(defaultLabelSet); } m_LabelSetContainer.push_back(labelSet); } } void mitk::LabelSetImage::SetActiveLayer(unsigned int layer) { try { if (4 == this->GetDimension()) { if ((layer != GetActiveLayer() || m_activeLayerInvalid) && (layer < this->GetNumberOfLayers())) { BeforeChangeLayerEvent.Send(); if (m_activeLayerInvalid) { // We should not write the invalid layer back to the vector m_activeLayerInvalid = false; } else { AccessFixedDimensionByItk_n(this, ImageToLayerContainerProcessing, 4, (GetActiveLayer())); } m_ActiveLayer = layer; // only at this place m_ActiveLayer should be manipulated!!! Use Getter and Setter AccessFixedDimensionByItk_n(this, LayerContainerToImageProcessing, 4, (GetActiveLayer())); AfterChangeLayerEvent.Send(); } } else { if ((layer != GetActiveLayer() || m_activeLayerInvalid) && (layer < this->GetNumberOfLayers())) { BeforeChangeLayerEvent.Send(); if (m_activeLayerInvalid) { // We should not write the invalid layer back to the vector m_activeLayerInvalid = false; } else { AccessByItk_1(this, ImageToLayerContainerProcessing, GetActiveLayer()); } m_ActiveLayer = layer; // only at this place m_ActiveLayer should be manipulated!!! Use Getter and Setter AccessByItk_1(this, LayerContainerToImageProcessing, GetActiveLayer()); AfterChangeLayerEvent.Send(); } } } catch (itk::ExceptionObject &e) { mitkThrow() << e.GetDescription(); } this->Modified(); } void mitk::LabelSetImage::Concatenate(mitk::LabelSetImage *other) { const unsigned int *otherDims = other->GetDimensions(); const unsigned int *thisDims = this->GetDimensions(); if ((otherDims[0] != thisDims[0]) || (otherDims[1] != thisDims[1]) || (otherDims[2] != thisDims[2])) mitkThrow() << "Dimensions do not match."; try { int numberOfLayers = other->GetNumberOfLayers(); for (int layer = 0; layer < numberOfLayers; ++layer) { this->SetActiveLayer(layer); AccessByItk_1(this, ConcatenateProcessing, other); mitk::LabelSet *ls = other->GetLabelSet(layer); auto it = ls->IteratorConstBegin(); auto end = ls->IteratorConstEnd(); it++; // skip exterior while (it != end) { GetLabelSet()->AddLabel((it->second)); // AddLabelEvent.Send(); it++; } } } catch (itk::ExceptionObject &e) { mitkThrow() << e.GetDescription(); } this->Modified(); } void mitk::LabelSetImage::ClearBuffer() { try { - AccessByItk(this, ClearBufferProcessing); + if (this->GetDimension() == 4) + { //remark: this extra branch was added, because LabelSetImage instances can be + //dynamic (4D), but AccessByItk by support only supports 2D and 3D. + //The option to change the CMake default dimensions for AccessByItk was + //dropped (for details see discussion in T28756) + AccessFixedDimensionByItk(this, ClearBufferProcessing,4); + } + else + { + AccessByItk(this, ClearBufferProcessing); + } this->Modified(); } catch (itk::ExceptionObject &e) { mitkThrow() << e.GetDescription(); } } bool mitk::LabelSetImage::ExistLabel(PixelType pixelValue) const { bool exist = false; for (unsigned int lidx = 0; lidx < GetNumberOfLayers(); lidx++) exist |= m_LabelSetContainer[lidx]->ExistLabel(pixelValue); return exist; } bool mitk::LabelSetImage::ExistLabel(PixelType pixelValue, unsigned int layer) const { bool exist = m_LabelSetContainer[layer]->ExistLabel(pixelValue); return exist; } bool mitk::LabelSetImage::ExistLabelSet(unsigned int layer) const { return layer < m_LabelSetContainer.size(); } void mitk::LabelSetImage::MergeLabel(PixelType pixelValue, PixelType sourcePixelValue, unsigned int layer) { try { AccessByItk_2(this, MergeLabelProcessing, pixelValue, sourcePixelValue); } catch (itk::ExceptionObject &e) { mitkThrow() << e.GetDescription(); } GetLabelSet(layer)->SetActiveLabel(pixelValue); Modified(); } void mitk::LabelSetImage::MergeLabels(PixelType pixelValue, std::vector& vectorOfSourcePixelValues, unsigned int layer) { try { for (unsigned int idx = 0; idx < vectorOfSourcePixelValues.size(); idx++) { AccessByItk_2(this, MergeLabelProcessing, pixelValue, vectorOfSourcePixelValues[idx]); } } catch (itk::ExceptionObject &e) { mitkThrow() << e.GetDescription(); } GetLabelSet(layer)->SetActiveLabel(pixelValue); Modified(); } void mitk::LabelSetImage::RemoveLabels(std::vector &VectorOfLabelPixelValues, unsigned int layer) { for (unsigned int idx = 0; idx < VectorOfLabelPixelValues.size(); idx++) { GetLabelSet(layer)->RemoveLabel(VectorOfLabelPixelValues[idx]); EraseLabel(VectorOfLabelPixelValues[idx], layer); } } void mitk::LabelSetImage::EraseLabels(std::vector &VectorOfLabelPixelValues, unsigned int layer) { for (unsigned int i = 0; i < VectorOfLabelPixelValues.size(); i++) { this->EraseLabel(VectorOfLabelPixelValues[i], layer); } } void mitk::LabelSetImage::EraseLabel(PixelType pixelValue, unsigned int layer) { try { if (4 == this->GetDimension()) { AccessFixedDimensionByItk_2(this, EraseLabelProcessing, 4, pixelValue, layer); } else { AccessByItk_2(this, EraseLabelProcessing, pixelValue, layer); } } catch (const itk::ExceptionObject &e) { mitkThrow() << e.GetDescription(); } Modified(); } mitk::Label *mitk::LabelSetImage::GetActiveLabel(unsigned int layer) { if (m_LabelSetContainer.size() <= layer) return nullptr; else return m_LabelSetContainer[layer]->GetActiveLabel(); } const mitk::Label* mitk::LabelSetImage::GetActiveLabel(unsigned int layer) const { if (m_LabelSetContainer.size() <= layer) return nullptr; else return m_LabelSetContainer[layer]->GetActiveLabel(); } mitk::Label *mitk::LabelSetImage::GetLabel(PixelType pixelValue, unsigned int layer) const { if (m_LabelSetContainer.size() <= layer) return nullptr; else return m_LabelSetContainer[layer]->GetLabel(pixelValue); } mitk::LabelSet *mitk::LabelSetImage::GetLabelSet(unsigned int layer) { if (m_LabelSetContainer.size() <= layer) return nullptr; else return m_LabelSetContainer[layer].GetPointer(); } const mitk::LabelSet *mitk::LabelSetImage::GetLabelSet(unsigned int layer) const { if (m_LabelSetContainer.size() <= layer) return nullptr; else return m_LabelSetContainer[layer].GetPointer(); } mitk::LabelSet *mitk::LabelSetImage::GetActiveLabelSet() { if (m_LabelSetContainer.size() == 0) return nullptr; else return m_LabelSetContainer[GetActiveLayer()].GetPointer(); } const mitk::LabelSet* mitk::LabelSetImage::GetActiveLabelSet() const { if (m_LabelSetContainer.size() == 0) return nullptr; else return m_LabelSetContainer[GetActiveLayer()].GetPointer(); } void mitk::LabelSetImage::UpdateCenterOfMass(PixelType pixelValue, unsigned int layer) { if (4 == this->GetDimension()) { AccessFixedDimensionByItk_2(this, CalculateCenterOfMassProcessing, 4, pixelValue, layer); } else { AccessByItk_2(this, CalculateCenterOfMassProcessing, pixelValue, layer); } } unsigned int mitk::LabelSetImage::GetNumberOfLabels(unsigned int layer) const { return m_LabelSetContainer[layer]->GetNumberOfLabels(); } unsigned int mitk::LabelSetImage::GetTotalNumberOfLabels() const { unsigned int totalLabels(0); auto layerIter = m_LabelSetContainer.begin(); for (; layerIter != m_LabelSetContainer.end(); ++layerIter) totalLabels += (*layerIter)->GetNumberOfLabels(); return totalLabels; } void mitk::LabelSetImage::MaskStamp(mitk::Image *mask, bool forceOverwrite) { try { mitk::PadImageFilter::Pointer padImageFilter = mitk::PadImageFilter::New(); padImageFilter->SetInput(0, mask); padImageFilter->SetInput(1, this); padImageFilter->SetPadConstant(0); padImageFilter->SetBinaryFilter(false); padImageFilter->SetLowerThreshold(0); padImageFilter->SetUpperThreshold(1); padImageFilter->Update(); mitk::Image::Pointer paddedMask = padImageFilter->GetOutput(); if (paddedMask.IsNull()) return; AccessByItk_2(this, MaskStampProcessing, paddedMask, forceOverwrite); } catch (...) { mitkThrow() << "Could not stamp the provided mask on the selected label."; } } mitk::Image::Pointer mitk::LabelSetImage::CreateLabelMask(PixelType index, bool useActiveLayer, unsigned int layer) { auto previousActiveLayer = this->GetActiveLayer(); auto mask = mitk::Image::New(); try { mask->Initialize(this); auto byteSize = sizeof(LabelSetImage::PixelType); for (unsigned int dim = 0; dim < mask->GetDimension(); ++dim) byteSize *= mask->GetDimension(dim); { ImageWriteAccessor accessor(mask); memset(accessor.GetData(), 0, byteSize); } if (!useActiveLayer) this->SetActiveLayer(layer); if (4 == this->GetDimension()) { ::CreateLabelMaskProcessing<4>(this, mask, index); } else if (3 == this->GetDimension()) { ::CreateLabelMaskProcessing(this, mask, index); } else { mitkThrow(); } } catch (...) { if (!useActiveLayer) this->SetActiveLayer(previousActiveLayer); mitkThrow() << "Could not create a mask out of the selected label."; } if (!useActiveLayer) this->SetActiveLayer(previousActiveLayer); return mask; } void mitk::LabelSetImage::InitializeByLabeledImage(mitk::Image::Pointer image) { if (image.IsNull() || image->IsEmpty() || !image->IsInitialized()) mitkThrow() << "Invalid labeled image."; try { this->Initialize(image); unsigned int byteSize = sizeof(LabelSetImage::PixelType); for (unsigned int dim = 0; dim < image->GetDimension(); ++dim) { byteSize *= image->GetDimension(dim); } mitk::ImageWriteAccessor *accessor = new mitk::ImageWriteAccessor(static_cast(this)); memset(accessor->GetData(), 0, byteSize); delete accessor; auto geometry = image->GetTimeGeometry()->Clone(); this->SetTimeGeometry(geometry); if (image->GetDimension() == 3) { AccessTwoImagesFixedDimensionByItk(this, image, InitializeByLabeledImageProcessing, 3); } else if (image->GetDimension() == 4) { AccessTwoImagesFixedDimensionByItk(this, image, InitializeByLabeledImageProcessing, 4); } else { mitkThrow() << image->GetDimension() << "-dimensional label set images not yet supported"; } } catch (...) { mitkThrow() << "Could not intialize by provided labeled image."; } this->Modified(); } template void mitk::LabelSetImage::InitializeByLabeledImageProcessing(LabelSetImageType *labelSetImage, ImageType *image) { typedef itk::ImageRegionConstIteratorWithIndex SourceIteratorType; typedef itk::ImageRegionIterator TargetIteratorType; TargetIteratorType targetIter(labelSetImage, labelSetImage->GetRequestedRegion()); targetIter.GoToBegin(); SourceIteratorType sourceIter(image, image->GetRequestedRegion()); sourceIter.GoToBegin(); while (!sourceIter.IsAtEnd()) { auto sourceValue = static_cast(sourceIter.Get()); targetIter.Set(sourceValue); if (!this->ExistLabel(sourceValue)) { std::stringstream name; name << "object-" << sourceValue; double rgba[4]; m_LabelSetContainer[this->GetActiveLayer()]->GetLookupTable()->GetTableValue(sourceValue, rgba); mitk::Color color; color.SetRed(rgba[0]); color.SetGreen(rgba[1]); color.SetBlue(rgba[2]); auto label = mitk::Label::New(); label->SetName(name.str().c_str()); label->SetColor(color); label->SetOpacity(rgba[3]); label->SetValue(sourceValue); this->GetLabelSet()->AddLabel(label); if (GetActiveLabelSet()->GetNumberOfLabels() >= mitk::Label::MAX_LABEL_VALUE || sourceValue >= mitk::Label::MAX_LABEL_VALUE) this->AddLayer(); } ++sourceIter; ++targetIter; } } template void mitk::LabelSetImage::MaskStampProcessing(ImageType *itkImage, mitk::Image *mask, bool forceOverwrite) { typename ImageType::Pointer itkMask; mitk::CastToItkImage(mask, itkMask); typedef itk::ImageRegionConstIterator SourceIteratorType; typedef itk::ImageRegionIterator TargetIteratorType; SourceIteratorType sourceIter(itkMask, itkMask->GetLargestPossibleRegion()); sourceIter.GoToBegin(); TargetIteratorType targetIter(itkImage, itkImage->GetLargestPossibleRegion()); targetIter.GoToBegin(); int activeLabel = this->GetActiveLabel(GetActiveLayer())->GetValue(); while (!sourceIter.IsAtEnd()) { PixelType sourceValue = sourceIter.Get(); PixelType targetValue = targetIter.Get(); if ((sourceValue != 0) && (forceOverwrite || !this->GetLabel(targetValue)->GetLocked())) // skip exterior and locked labels { targetIter.Set(activeLabel); } ++sourceIter; ++targetIter; } this->Modified(); } template void mitk::LabelSetImage::CalculateCenterOfMassProcessing(ImageType *itkImage, PixelType pixelValue, unsigned int layer) { // for now, we just retrieve the voxel in the middle typedef itk::ImageRegionConstIterator IteratorType; IteratorType iter(itkImage, itkImage->GetLargestPossibleRegion()); iter.GoToBegin(); std::vector indexVector; while (!iter.IsAtEnd()) { // TODO fix comparison warning more effective if (iter.Get() == pixelValue) { indexVector.push_back(iter.GetIndex()); } ++iter; } mitk::Point3D pos; pos.Fill(0.0); if (!indexVector.empty()) { typename itk::ImageRegionConstIteratorWithIndex::IndexType centerIndex; centerIndex = indexVector.at(indexVector.size() / 2); if (centerIndex.GetIndexDimension() == 3) { pos[0] = centerIndex[0]; pos[1] = centerIndex[1]; pos[2] = centerIndex[2]; } else return; } GetLabelSet(layer)->GetLabel(pixelValue)->SetCenterOfMassIndex(pos); this->GetSlicedGeometry()->IndexToWorld(pos, pos); // TODO: TimeGeometry? GetLabelSet(layer)->GetLabel(pixelValue)->SetCenterOfMassCoordinates(pos); } template void mitk::LabelSetImage::ClearBufferProcessing(ImageType *itkImage) { itkImage->FillBuffer(0); } // todo: concatenate all layers and not just the active one template void mitk::LabelSetImage::ConcatenateProcessing(ImageType *itkTarget, mitk::LabelSetImage *other) { typename ImageType::Pointer itkSource = ImageType::New(); mitk::CastToItkImage(other, itkSource); typedef itk::ImageRegionConstIterator ConstIteratorType; typedef itk::ImageRegionIterator IteratorType; ConstIteratorType sourceIter(itkSource, itkSource->GetLargestPossibleRegion()); IteratorType targetIter(itkTarget, itkTarget->GetLargestPossibleRegion()); int numberOfTargetLabels = this->GetNumberOfLabels(GetActiveLayer()) - 1; // skip exterior sourceIter.GoToBegin(); targetIter.GoToBegin(); while (!sourceIter.IsAtEnd()) { PixelType sourceValue = sourceIter.Get(); PixelType targetValue = targetIter.Get(); if ((sourceValue != 0) && !this->GetLabel(targetValue)->GetLocked()) // skip exterior and locked labels { targetIter.Set(sourceValue + numberOfTargetLabels); } ++sourceIter; ++targetIter; } } template void mitk::LabelSetImage::LayerContainerToImageProcessing(itk::Image *target, unsigned int layer) { typedef itk::Image ImageType; typename ImageType::Pointer itkSource; // mitk::CastToItkImage(m_LayerContainer[layer], itkSource); itkSource = ImageToItkImage(m_LayerContainer[layer]); typedef itk::ImageRegionConstIterator SourceIteratorType; typedef itk::ImageRegionIterator TargetIteratorType; SourceIteratorType sourceIter(itkSource, itkSource->GetLargestPossibleRegion()); sourceIter.GoToBegin(); TargetIteratorType targetIter(target, target->GetLargestPossibleRegion()); targetIter.GoToBegin(); while (!sourceIter.IsAtEnd()) { targetIter.Set(sourceIter.Get()); ++sourceIter; ++targetIter; } } template void mitk::LabelSetImage::ImageToLayerContainerProcessing(itk::Image *source, unsigned int layer) const { typedef itk::Image ImageType; typename ImageType::Pointer itkTarget; // mitk::CastToItkImage(m_LayerContainer[layer], itkTarget); itkTarget = ImageToItkImage(m_LayerContainer[layer]); typedef itk::ImageRegionConstIterator SourceIteratorType; typedef itk::ImageRegionIterator TargetIteratorType; SourceIteratorType sourceIter(source, source->GetLargestPossibleRegion()); sourceIter.GoToBegin(); TargetIteratorType targetIter(itkTarget, itkTarget->GetLargestPossibleRegion()); targetIter.GoToBegin(); while (!sourceIter.IsAtEnd()) { targetIter.Set(sourceIter.Get()); ++sourceIter; ++targetIter; } } template void mitk::LabelSetImage::EraseLabelProcessing(ImageType *itkImage, PixelType pixelValue, unsigned int /*layer*/) { typedef itk::ImageRegionIterator IteratorType; IteratorType iter(itkImage, itkImage->GetLargestPossibleRegion()); iter.GoToBegin(); while (!iter.IsAtEnd()) { PixelType value = iter.Get(); if (value == pixelValue) { iter.Set(0); } ++iter; } } template void mitk::LabelSetImage::MergeLabelProcessing(ImageType *itkImage, PixelType pixelValue, PixelType index) { typedef itk::ImageRegionIterator IteratorType; IteratorType iter(itkImage, itkImage->GetLargestPossibleRegion()); iter.GoToBegin(); while (!iter.IsAtEnd()) { if (iter.Get() == index) { iter.Set(pixelValue); } ++iter; } } bool mitk::Equal(const mitk::LabelSetImage &leftHandSide, const mitk::LabelSetImage &rightHandSide, ScalarType eps, bool verbose) { bool returnValue = true; /* LabelSetImage members */ MITK_INFO(verbose) << "--- LabelSetImage Equal ---"; // number layers returnValue = leftHandSide.GetNumberOfLayers() == rightHandSide.GetNumberOfLayers(); if (!returnValue) { MITK_INFO(verbose) << "Number of layers not equal."; return false; } // total number labels returnValue = leftHandSide.GetTotalNumberOfLabels() == rightHandSide.GetTotalNumberOfLabels(); if (!returnValue) { MITK_INFO(verbose) << "Total number of labels not equal."; return false; } // active layer returnValue = leftHandSide.GetActiveLayer() == rightHandSide.GetActiveLayer(); if (!returnValue) { MITK_INFO(verbose) << "Active layer not equal."; return false; } if (4 == leftHandSide.GetDimension()) { MITK_INFO(verbose) << "Can not compare image data for 4D images - skipping check."; } else { // working image data returnValue = mitk::Equal((const mitk::Image &)leftHandSide, (const mitk::Image &)rightHandSide, eps, verbose); if (!returnValue) { MITK_INFO(verbose) << "Working image data not equal."; return false; } } for (unsigned int layerIndex = 0; layerIndex < leftHandSide.GetNumberOfLayers(); layerIndex++) { if (4 == leftHandSide.GetDimension()) { MITK_INFO(verbose) << "Can not compare image data for 4D images - skipping check."; } else { // layer image data returnValue = mitk::Equal(*leftHandSide.GetLayerImage(layerIndex), *rightHandSide.GetLayerImage(layerIndex), eps, verbose); if (!returnValue) { MITK_INFO(verbose) << "Layer image data not equal."; return false; } } // layer labelset data returnValue = mitk::Equal(*leftHandSide.GetLabelSet(layerIndex), *rightHandSide.GetLabelSet(layerIndex), eps, verbose); if (!returnValue) { MITK_INFO(verbose) << "Layer labelset data not equal."; return false; } } return returnValue; }