X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkSegmentationUtils.txx;h=24ca5194853595a1ac2b815aeeee60533f5412e1;hb=1d3d35c2ac5c6b74d24d9f56d53de6aa990915d1;hp=17603959a99fdf3a7246df225073bb34e1dd3082;hpb=46acae08f9deb8c240a11b6106601e1e72f29478;p=clitk.git diff --git a/itk/clitkSegmentationUtils.txx b/itk/clitkSegmentationUtils.txx index 1760395..24ca519 100644 --- a/itk/clitkSegmentationUtils.txx +++ b/itk/clitkSegmentationUtils.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -32,762 +32,1366 @@ #include #include #include +#include +#include +#include namespace clitk { -//-------------------------------------------------------------------- -template -void ComputeBBFromImageRegion(typename ImageType::Pointer image, - typename ImageType::RegionType region, - typename itk::BoundingBox::Pointer bb) { - typedef typename ImageType::IndexType IndexType; - IndexType firstIndex; - IndexType lastIndex; - for(unsigned int i=0; iGetImageDimension(); i++) { - firstIndex[i] = region.GetIndex()[i]; - lastIndex[i] = firstIndex[i]+region.GetSize()[i]; - } - - typedef itk::BoundingBox BBType; - typedef typename BBType::PointType PointType; - PointType lastPoint; - PointType firstPoint; - image->TransformIndexToPhysicalPoint(firstIndex, firstPoint); - image->TransformIndexToPhysicalPoint(lastIndex, lastPoint); - - bb->SetMaximum(lastPoint); - bb->SetMinimum(firstPoint); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ComputeBBIntersection(typename itk::BoundingBox::Pointer bbo, - typename itk::BoundingBox::Pointer bbi1, - typename itk::BoundingBox::Pointer bbi2) { - - typedef itk::BoundingBox BBType; - typedef typename BBType::PointType PointType; - PointType lastPoint; - PointType firstPoint; - - for(unsigned int i=0; iGetMinimum()[i], - bbi2->GetMinimum()[i]); - lastPoint[i] = std::min(bbi1->GetMaximum()[i], - bbi2->GetMaximum()[i]); - } - - bbo->SetMaximum(lastPoint); - bbo->SetMinimum(firstPoint); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ComputeRegionFromBB(typename ImageType::Pointer image, - const typename itk::BoundingBox::Pointer bb, - typename ImageType::RegionType & region) { - // Types - typedef typename ImageType::IndexType IndexType; - typedef typename ImageType::PointType PointType; - typedef typename ImageType::RegionType RegionType; - typedef typename ImageType::SizeType SizeType; - - // Region starting point - IndexType regionStart; - PointType start = bb->GetMinimum(); - image->TransformPhysicalPointToIndex(start, regionStart); - - // Region size - SizeType regionSize; - PointType maxs = bb->GetMaximum(); - PointType mins = bb->GetMinimum(); - for(unsigned int i=0; iGetSpacing()[i]); - regionSize[i] = lrint((maxs[i] - mins[i])/image->GetSpacing()[i]); - // DD(regionSize[i]); - } - - // Create region - region.SetIndex(regionStart); - region.SetSize(regionSize); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -SetBackground(const ImageType * input, - const TMaskImageType * mask, - typename TMaskImageType::PixelType maskBG, - typename ImageType::PixelType outValue, - bool inPlace) { - typedef SetBackgroundImageFilter - SetBackgroundImageFilterType; - typename SetBackgroundImageFilterType::Pointer setBackgroundFilter - = SetBackgroundImageFilterType::New(); - // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault - setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low - setBackgroundFilter->SetInput(input); - setBackgroundFilter->SetInput2(mask); - setBackgroundFilter->SetMaskValue(maskBG); - setBackgroundFilter->SetOutsideValue(outValue); - setBackgroundFilter->Update(); - return setBackgroundFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -int GetNumberOfConnectedComponentLabels(typename ImageType::Pointer input, - typename ImageType::PixelType BG, - bool isFullyConnected) { - // Connected Component label - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); - connectFilter->SetInput(input); - connectFilter->SetBackgroundValue(BG); - connectFilter->SetFullyConnected(isFullyConnected); - connectFilter->Update(); - - // Return result - return connectFilter->GetObjectCount(); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -/* - Warning : in this cas, we consider outputType like inputType, not - InternalImageType. Be sure it fits. - */ -template -typename ImageType::Pointer -Labelize(const ImageType * input, - typename ImageType::PixelType BG, - bool isFullyConnected, - int minimalComponentSize) { - // InternalImageType for storing large number of component - typedef itk::Image InternalImageType; - - // Connected Component label - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); - // connectFilter->ReleaseDataFlagOn(); - connectFilter->SetInput(input); - connectFilter->SetBackgroundValue(BG); - connectFilter->SetFullyConnected(isFullyConnected); + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + SetBackground(const ImageType * input, + const TMaskImageType * mask, + typename TMaskImageType::PixelType maskBG, + typename ImageType::PixelType outValue, + bool inPlace) { + typedef SetBackgroundImageFilter + SetBackgroundImageFilterType; + typename SetBackgroundImageFilterType::Pointer setBackgroundFilter + = SetBackgroundImageFilterType::New(); + // if (inPlace) setBackgroundFilter->ReleaseDataFlagOn(); // No seg fault + setBackgroundFilter->SetInPlace(inPlace); // This is important to keep memory low + setBackgroundFilter->SetInput(input); + setBackgroundFilter->SetInput2(mask); + setBackgroundFilter->SetMaskValue(maskBG); + setBackgroundFilter->SetOutsideValue(outValue); + setBackgroundFilter->Update(); + return setBackgroundFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + int GetNumberOfConnectedComponentLabels(const ImageType * input, + typename ImageType::PixelType BG, + bool isFullyConnected) { + // Connected Component label + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); + connectFilter->SetInput(input); + connectFilter->SetBackgroundValue(BG); + connectFilter->SetFullyConnected(isFullyConnected); + connectFilter->Update(); - // Sort by size and remove too small area. - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); - // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? - relabelFilter->SetInput(connectFilter->GetOutput()); - relabelFilter->SetMinimumObjectSize(minimalComponentSize); - relabelFilter->Update(); - - // DD(relabelFilter->GetNumberOfObjects()); - // DD(relabelFilter->GetOriginalNumberOfObjects()); - // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]); - - // Return result - typename ImageType::Pointer output = relabelFilter->GetOutput(); - return output; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -/* - Warning : in this cas, we consider outputType like inputType, not - InternalImageType. Be sure it fits. - */ -template -typename ImageType::Pointer -LabelizeAndCountNumberOfObjects(const ImageType * input, - typename ImageType::PixelType BG, - bool isFullyConnected, - int minimalComponentSize, - int & nb) { - // InternalImageType for storing large number of component - typedef itk::Image InternalImageType; + // Return result + return connectFilter->GetObjectCount(); + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + /* + Warning : in this cas, we consider outputType like inputType, not + InternalImageType. Be sure it fits. + */ + template + typename ImageType::Pointer + Labelize(const ImageType * input, + typename ImageType::PixelType BG, + bool isFullyConnected, + int minimalComponentSize) { + // InternalImageType for storing large number of component + typedef itk::Image InternalImageType; - // Connected Component label - typedef itk::ConnectedComponentImageFilter ConnectFilterType; - typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); - // connectFilter->ReleaseDataFlagOn(); - connectFilter->SetInput(input); - connectFilter->SetBackgroundValue(BG); - connectFilter->SetFullyConnected(isFullyConnected); + // Connected Component label + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); + // connectFilter->ReleaseDataFlagOn(); + connectFilter->SetInput(input); + connectFilter->SetBackgroundValue(BG); + connectFilter->SetFullyConnected(isFullyConnected); - // Sort by size and remove too small area. - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); - // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? - relabelFilter->SetInput(connectFilter->GetOutput()); - relabelFilter->SetMinimumObjectSize(minimalComponentSize); - relabelFilter->Update(); - - nb = relabelFilter->GetNumberOfObjects(); - // DD(relabelFilter->GetNumberOfObjects()); - // DD(relabelFilter->GetOriginalNumberOfObjects()); - // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]); - - // Return result - typename ImageType::Pointer output = relabelFilter->GetOutput(); - return output; -} -//-------------------------------------------------------------------- - - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -RemoveLabels(typename ImageType::Pointer input, - typename ImageType::PixelType BG, - std::vector & labelsToRemove) { - typename ImageType::Pointer working_image = input; - for (unsigned int i=0; i SetBackgroundImageFilterType; - typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New(); - setBackgroundFilter->SetInput(input); - setBackgroundFilter->SetInput2(input); - setBackgroundFilter->SetMaskValue(labelsToRemove[i]); - setBackgroundFilter->SetOutsideValue(BG); - setBackgroundFilter->Update(); - working_image = setBackgroundFilter->GetOutput(); + // Sort by size and remove too small area. + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); + // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? + relabelFilter->SetInput(connectFilter->GetOutput()); + relabelFilter->SetMinimumObjectSize(minimalComponentSize); + relabelFilter->Update(); + + // Return result + typename ImageType::Pointer output = relabelFilter->GetOutput(); + return output; } - return working_image; -} -//-------------------------------------------------------------------- + //-------------------------------------------------------------------- -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -KeepLabels(const ImageType * input, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - typename ImageType::PixelType firstKeep, - typename ImageType::PixelType lastKeep, - bool useLastKeep) { - typedef itk::BinaryThresholdImageFilter BinarizeFilterType; - typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); - binarizeFilter->SetInput(input); - binarizeFilter->SetLowerThreshold(firstKeep); - if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep); - binarizeFilter->SetInsideValue(FG); - binarizeFilter->SetOutsideValue(BG); - binarizeFilter->Update(); - return binarizeFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -LabelizeAndSelectLabels(typename ImageType::Pointer input, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool isFullyConnected, - int minimalComponentSize, - LabelizeParameters * param) -{ - typename ImageType::Pointer working_image; - working_image = Labelize(input, BG, isFullyConnected, minimalComponentSize); - working_image = RemoveLabels(working_image, BG, param->GetLabelsToRemove()); - working_image = KeepLabels(working_image, - BG, FG, - param->GetFirstKeep(), - param->GetLastKeep(), - param->GetUseLastKeep()); - return working_image; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -ResizeImageLike(typename ImageType::Pointer input, - typename ImageType::Pointer like, - typename ImageType::PixelType backgroundValue) -{ - typedef CropLikeImageFilter CropFilterType; - typename CropFilterType::Pointer cropFilter = CropFilterType::New(); - cropFilter->SetInput(input); - cropFilter->SetCropLikeImage(like); - cropFilter->SetBackgroundValue(backgroundValue); - cropFilter->Update(); - return cropFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename MaskImageType::Pointer -SliceBySliceRelativePosition(const MaskImageType * input, - const MaskImageType * object, - int direction, - double threshold, - std::string orientation, - bool uniqueConnectedComponent, - double spacing, - bool inverseflag) -{ - typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; - typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); - sliceRelPosFilter->VerboseStepFlagOff(); - sliceRelPosFilter->WriteStepFlagOff(); - sliceRelPosFilter->SetInput(input); - sliceRelPosFilter->SetInputObject(object); - sliceRelPosFilter->SetDirection(direction); - sliceRelPosFilter->SetFuzzyThreshold(threshold); - sliceRelPosFilter->AddOrientationTypeString(orientation); - sliceRelPosFilter->SetResampleBeforeRelativePositionFilter((spacing != -1)); - sliceRelPosFilter->SetIntermediateSpacing(spacing); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent); - sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); - // sliceRelPosFilter->SetAutoCropFlag(true); ?? - sliceRelPosFilter->Update(); - return sliceRelPosFilter->GetOutput(); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -bool -FindExtremaPointInAGivenDirection(const ImageType * input, - typename ImageType::PixelType bg, - int direction, bool opposite, - typename ImageType::PointType & point) -{ - typename ImageType::PointType dummy; - return FindExtremaPointInAGivenDirection(input, bg, direction, - opposite, dummy, 0, point); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -bool -FindExtremaPointInAGivenDirection(const ImageType * input, - typename ImageType::PixelType bg, - int direction, bool opposite, - typename ImageType::PointType refpoint, - double distanceMax, - typename ImageType::PointType & point) -{ + //-------------------------------------------------------------------- /* - loop over input pixels, store the index in the fg that is max - according to the given direction. + Warning : in this cas, we consider outputType like inputType, not + InternalImageType. Be sure it fits. */ - typedef itk::ImageRegionConstIteratorWithIndex IteratorType; - IteratorType iter(input, input->GetLargestPossibleRegion()); - iter.GoToBegin(); - typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex(); - if (opposite) max = max+input->GetLargestPossibleRegion().GetSize(); - bool found=false; - while (!iter.IsAtEnd()) { - if (iter.Get() != bg) { - bool test = iter.GetIndex()[direction] > max[direction]; - if (opposite) test = !test; - if (test) { - typename ImageType::PointType p; - input->TransformIndexToPhysicalPoint(iter.GetIndex(), p); - if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) { - max = iter.GetIndex(); - found = true; + template + typename ImageType::Pointer + LabelizeAndCountNumberOfObjects(const ImageType * input, + typename ImageType::PixelType BG, + bool isFullyConnected, + int minimalComponentSize, + int & nb) { + // InternalImageType for storing large number of component + typedef itk::Image InternalImageType; + + // Connected Component label + typedef itk::ConnectedComponentImageFilter ConnectFilterType; + typename ConnectFilterType::Pointer connectFilter = ConnectFilterType::New(); + // connectFilter->ReleaseDataFlagOn(); + connectFilter->SetInput(input); + connectFilter->SetBackgroundValue(BG); + connectFilter->SetFullyConnected(isFullyConnected); + + // Sort by size and remove too small area. + typedef itk::RelabelComponentImageFilter RelabelFilterType; + typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New(); + // relabelFilter->ReleaseDataFlagOn(); // if yes, fail when ExplosionControlledThresholdConnectedImageFilter ??? + relabelFilter->SetInput(connectFilter->GetOutput()); + relabelFilter->SetMinimumObjectSize(minimalComponentSize); + relabelFilter->Update(); + + nb = relabelFilter->GetNumberOfObjects(); + // DD(relabelFilter->GetOriginalNumberOfObjects()); + // DD(relabelFilter->GetSizeOfObjectsInPhysicalUnits()[0]); + + // Return result + typename ImageType::Pointer output = relabelFilter->GetOutput(); + return output; + } + //-------------------------------------------------------------------- + + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + RemoveLabels(const ImageType * input, + typename ImageType::PixelType BG, + std::vector & labelsToRemove) { + assert(labelsToRemove.size() != 0); + typename ImageType::Pointer working_image;// = input; + for (unsigned int i=0; i SetBackgroundImageFilterType; + typename SetBackgroundImageFilterType::Pointer setBackgroundFilter = SetBackgroundImageFilterType::New(); + setBackgroundFilter->SetInput(input); + setBackgroundFilter->SetInput2(input); + setBackgroundFilter->SetMaskValue(labelsToRemove[i]); + setBackgroundFilter->SetOutsideValue(BG); + setBackgroundFilter->Update(); + working_image = setBackgroundFilter->GetOutput(); + } + return working_image; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + KeepLabels(const ImageType * input, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + typename ImageType::PixelType firstKeep, + typename ImageType::PixelType lastKeep, + bool useLastKeep) { + typedef itk::BinaryThresholdImageFilter BinarizeFilterType; + typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); + binarizeFilter->SetInput(input); + binarizeFilter->SetLowerThreshold(firstKeep); + if (useLastKeep) binarizeFilter->SetUpperThreshold(lastKeep); + binarizeFilter->SetInsideValue(FG); + binarizeFilter->SetOutsideValue(BG); + binarizeFilter->Update(); + return binarizeFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + LabelizeAndSelectLabels(const ImageType * input, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool isFullyConnected, + int minimalComponentSize, + LabelizeParameters * param) + { + typename ImageType::Pointer working_image; + working_image = Labelize(input, BG, isFullyConnected, minimalComponentSize); + if (param->GetLabelsToRemove().size() != 0) + working_image = RemoveLabels(working_image, BG, param->GetLabelsToRemove()); + working_image = KeepLabels(working_image, + BG, FG, + param->GetFirstKeep(), + param->GetLastKeep(), + param->GetUseLastKeep()); + return working_image; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename MaskImageType::Pointer + SliceBySliceRelativePosition(const MaskImageType * input, + const MaskImageType * object, + int direction, + double threshold, + std::string orientation, + bool uniqueConnectedComponent, + double spacing, + bool autocropFlag, + bool singleObjectCCL) + { + typedef clitk::SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->VerboseStepFlagOff(); + sliceRelPosFilter->WriteStepFlagOff(); + sliceRelPosFilter->SetInput(input); + sliceRelPosFilter->SetInputObject(object); + sliceRelPosFilter->SetDirection(direction); + sliceRelPosFilter->SetFuzzyThreshold(threshold); + sliceRelPosFilter->AddOrientationTypeString(orientation); + sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1)); + sliceRelPosFilter->SetIntermediateSpacing(spacing); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent); + sliceRelPosFilter->ObjectCCLSelectionFlagOff(); + sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL); + // sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); + sliceRelPosFilter->SetAutoCropFlag(autocropFlag); + sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); + sliceRelPosFilter->Update(); + return sliceRelPosFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename MaskImageType::Pointer + SliceBySliceRelativePosition(const MaskImageType * input, + const MaskImageType * object, + int direction, + double threshold, + double angle, + bool inverseflag, + bool uniqueConnectedComponent, + double spacing, + bool autocropFlag, + bool singleObjectCCL) + { + typedef clitk::SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->VerboseStepFlagOff(); + sliceRelPosFilter->WriteStepFlagOff(); + sliceRelPosFilter->SetInput(input); + sliceRelPosFilter->SetInputObject(object); + sliceRelPosFilter->SetDirection(direction); + sliceRelPosFilter->SetFuzzyThreshold(threshold); + // sliceRelPosFilter->AddOrientationTypeString(orientation); + sliceRelPosFilter->AddAnglesInRad(angle, 0.0); + sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1)); + sliceRelPosFilter->SetIntermediateSpacing(spacing); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent); + sliceRelPosFilter->ObjectCCLSelectionFlagOff(); + sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL); + sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); + sliceRelPosFilter->SetAutoCropFlag(autocropFlag); + sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); + sliceRelPosFilter->Update(); + return sliceRelPosFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + bool + FindExtremaPointInAGivenDirection(const ImageType * input, + typename ImageType::PixelType bg, + int direction, bool opposite, + typename ImageType::PointType & point) + { + typename ImageType::PointType dummy; + return FindExtremaPointInAGivenDirection(input, bg, direction, + opposite, dummy, 0, point); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + bool + FindExtremaPointInAGivenDirection(const ImageType * input, + typename ImageType::PixelType bg, + int direction, bool opposite, + typename ImageType::PointType refpoint, + double distanceMax, + typename ImageType::PointType & point) + { + /* + loop over input pixels, store the index in the fg that is max + according to the given direction. + */ + typedef itk::ImageRegionConstIteratorWithIndex IteratorType; + IteratorType iter(input, input->GetLargestPossibleRegion()); + iter.GoToBegin(); + typename ImageType::IndexType max = input->GetLargestPossibleRegion().GetIndex(); + if (opposite) max = max+input->GetLargestPossibleRegion().GetSize(); + bool found=false; + while (!iter.IsAtEnd()) { + if (iter.Get() != bg) { + bool test = iter.GetIndex()[direction] > max[direction]; + if (opposite) test = !test; + if (test) { + typename ImageType::PointType p; + input->TransformIndexToPhysicalPoint(iter.GetIndex(), p); + if ((distanceMax==0) || (p.EuclideanDistanceTo(refpoint) < distanceMax)) { + max = iter.GetIndex(); + found = true; + } } } + ++iter; } - ++iter; - } - if (!found) return false; - input->TransformIndexToPhysicalPoint(max, point); - return true; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -CropImageAbove(typename ImageType::Pointer image, - int dim, double min, - bool autoCrop, - typename ImageType::PixelType BG) -{ - return CropImageAlongOneAxis(image, dim, - image->GetOrigin()[dim], - min, - autoCrop, BG); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -CropImageBelow(typename ImageType::Pointer image, - int dim, double max, - bool autoCrop, - typename ImageType::PixelType BG) -{ - typename ImageType::PointType p; - image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+ - image->GetLargestPossibleRegion().GetSize(), p); - return CropImageAlongOneAxis(image, dim, max, p[dim], autoCrop, BG); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -CropImageAlongOneAxis(typename ImageType::Pointer image, - int dim, double min, double max, - bool autoCrop, - typename ImageType::PixelType BG) -{ - // Compute region size - typename ImageType::RegionType region; - typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize(); - typename ImageType::PointType p = image->GetOrigin(); - p[dim] = min; - typename ImageType::IndexType start; - image->TransformPhysicalPointToIndex(p, start); - p[dim] = max; - typename ImageType::IndexType end; - image->TransformPhysicalPointToIndex(p, end); - size[dim] = fabs(end[dim]-start[dim]); - region.SetIndex(start); - region.SetSize(size); + if (!found) return false; + input->TransformIndexToPhysicalPoint(max, point); + return true; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + CropImageRemoveGreaterThan(const ImageType * image, + int dim, double min, bool autoCrop, + typename ImageType::PixelType BG) + { + return CropImageAlongOneAxis(image, dim, + image->GetOrigin()[dim], + min, + autoCrop, BG); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + CropImageRemoveLowerThan(const ImageType * image, + int dim, double max, bool autoCrop, + typename ImageType::PixelType BG) + { + typename ImageType::PointType p; + image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+ + image->GetLargestPossibleRegion().GetSize(), p); + // Add GetSpacing because remove Lower or equal than + // DD(max); + // DD(p); + // DD(max+image->GetSpacing()[dim]); + return CropImageAlongOneAxis(image, dim, max+image->GetSpacing()[dim], p[dim], autoCrop, BG); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + CropImageAlongOneAxis(const ImageType * image, + int dim, double min, double max, + bool autoCrop, typename ImageType::PixelType BG) + { + // Compute region size + typename ImageType::RegionType region; + typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize(); + typename ImageType::PointType p = image->GetOrigin(); + if (min > p[dim]) p[dim] = min; // Check if not outside the image + typename ImageType::IndexType start; + image->TransformPhysicalPointToIndex(p, start); + double m = image->GetOrigin()[dim] + size[dim]*image->GetSpacing()[dim]; + if (max > m) p[dim] = m; // Check if not outside the image + else p[dim] = max; + typename ImageType::IndexType end; + image->TransformPhysicalPointToIndex(p, end); + size[dim] = abs(end[dim]-start[dim]); + region.SetIndex(start); + region.SetSize(size); - // Perform Crop - typedef itk::RegionOfInterestImageFilter CropFilterType; - typename CropFilterType::Pointer cropFilter = CropFilterType::New(); - cropFilter->SetInput(image); - cropFilter->SetRegionOfInterest(region); - cropFilter->Update(); - typename ImageType::Pointer result = cropFilter->GetOutput(); + // Perform Crop + typedef itk::RegionOfInterestImageFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(image); + cropFilter->SetRegionOfInterest(region); + cropFilter->Update(); + typename ImageType::Pointer result = cropFilter->GetOutput(); - // Auto Crop - if (autoCrop) { - result = AutoCrop(result, BG); - } - return result; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -ComputeCentroids(typename ImageType::Pointer image, - typename ImageType::PixelType BG, - std::vector & centroids) -{ - typedef long LabelType; - static const unsigned int Dim = ImageType::ImageDimension; - typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; - typedef itk::LabelMap< LabelObjectType > LabelMapType; - typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; - typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); - typedef itk::ShapeLabelMapFilter ShapeFilterType; - typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); - imageToLabelFilter->SetBackgroundValue(BG); - imageToLabelFilter->SetInput(image); - statFilter->SetInput(imageToLabelFilter->GetOutput()); - statFilter->Update(); - typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); - - centroids.clear(); - typename ImageType::PointType dummy; - centroids.push_back(dummy); // label 0 -> no centroid, use dummy point - for(uint i=1; iGetNumberOfLabelObjects()+1; i++) { - centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid()); - } -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -ExtractSlices(typename ImageType::Pointer image, - int direction, - std::vector::Pointer > & slices) -{ - typedef ExtractSliceFilter ExtractSliceFilterType; - typedef typename ExtractSliceFilterType::SliceType SliceType; - typename ExtractSliceFilterType::Pointer - extractSliceFilter = ExtractSliceFilterType::New(); - extractSliceFilter->SetInput(image); - extractSliceFilter->SetDirection(direction); - extractSliceFilter->Update(); - extractSliceFilter->GetOutputSlices(slices); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -JoinSlices(std::vector::Pointer > & slices, - typename ImageType::Pointer input, - int direction) { - typedef typename itk::Image SliceType; - typedef itk::JoinSeriesImageFilter JoinSeriesFilterType; - typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New(); - joinFilter->SetOrigin(input->GetOrigin()[direction]); - joinFilter->SetSpacing(input->GetSpacing()[direction]); - for(unsigned int i=0; iPushBackInput(slices[i]); - } - joinFilter->Update(); - return joinFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -PointsUtils::Convert2DTo3D(const PointType2D & p, - ImagePointer image, - const int slice, - PointType3D & p3D) -{ - p3D[0] = p[0]; - p3D[1] = p[1]; - p3D[2] = (image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2] - + image->GetOrigin()[2]; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -PointsUtils::Convert2DTo3DList(const MapPoint2DType & map, - ImagePointer image, - VectorPoint3DType & list) -{ - typename MapPoint2DType::const_iterator iter = map.begin(); - while (iter != map.end()) { - PointType3D p; - Convert2DTo3D(iter->second, image, iter->first, p); - list.push_back(p); - ++iter; - } -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -void -WriteListOfLandmarks(std::vector points, - std::string filename) -{ - std::ofstream os; - openFileForWriting(os, filename); - os << "LANDMARKS1" << std::endl; - for(uint i=0; i -typename ImageType::Pointer -Dilate(typename ImageType::Pointer image, - double radiusInMM, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport) -{ - typename ImageType::SizeType r; - for(uint i=0; iGetSpacing()[i]); - return Dilate(image, r, BG, FG, extendSupport); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -Dilate(typename ImageType::Pointer image, - typename ImageType::PointType radiusInMM, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport) -{ - typename ImageType::SizeType r; - for(uint i=0; iGetSpacing()[i]); - return Dilate(image, r, BG, FG, extendSupport); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -typename ImageType::Pointer -Dilate(typename ImageType::Pointer image, - typename ImageType::SizeType radius, - typename ImageType::PixelType BG, - typename ImageType::PixelType FG, - bool extendSupport) -{ - // Create kernel for dilatation - typedef itk::BinaryBallStructuringElement KernelType; - KernelType structuringElement; - structuringElement.SetRadius(radius); - structuringElement.CreateStructuringElement(); - - if (extendSupport) { - typedef itk::ConstantPadImageFilter PadFilterType; - typename PadFilterType::Pointer padFilter = PadFilterType::New(); - padFilter->SetInput(image); - typename ImageType::SizeType lower; - typename ImageType::SizeType upper; - for(uint i=0; i<3; i++) { - lower[i] = upper[i] = 2*(radius[i]+1); + // Auto Crop + if (autoCrop) { + result = AutoCrop(result, BG); } - padFilter->SetPadLowerBound(lower); - padFilter->SetPadUpperBound(upper); - padFilter->Update(); - image = padFilter->GetOutput(); - } - - // Dilate filter - typedef itk::BinaryDilateImageFilter DilateFilterType; - typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New(); - dilateFilter->SetBackgroundValue(BG); - dilateFilter->SetForegroundValue(FG); - dilateFilter->SetBoundaryToForeground(false); - dilateFilter->SetKernel(structuringElement); - dilateFilter->SetInput(image); - dilateFilter->Update(); - return image = dilateFilter->GetOutput(); -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ConvertOption(std::string optionName, uint given, - ValueType * values, VectorType & p, - uint dim, bool required) -{ - if (required && (given == 0)) { - clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or " + return result; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ComputeCentroids(const ImageType * image, + typename ImageType::PixelType BG, + std::vector & centroids) + { + typedef long LabelType; + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; + typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); + typedef itk::ShapeLabelMapFilter ShapeFilterType; + typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); + imageToLabelFilter->SetBackgroundValue(BG); + imageToLabelFilter->SetInput(image); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->Update(); + typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); + + centroids.clear(); + typename ImageType::PointType dummy; + centroids.push_back(dummy); // label 0 -> no centroid, use dummy point for BG + //DS FIXME (not useful ! to change ..) + for(uint i=0; iGetNumberOfLabelObjects(); i++) { + int label = labelMap->GetLabels()[i]; + centroids.push_back(labelMap->GetLabelObject(label)->GetCentroid()); + } + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename itk::LabelMap< itk::ShapeLabelObject >::Pointer + ComputeLabelMap(const ImageType * image, + typename ImageType::PixelType BG, + bool computePerimeterFlag) + { + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; + typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); + typedef itk::ShapeLabelMapFilter ShapeFilterType; + typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); + imageToLabelFilter->SetBackgroundValue(BG); + imageToLabelFilter->SetInput(image); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->SetComputePerimeter(computePerimeterFlag); + statFilter->Update(); + return statFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + ComputeCentroids2(const ImageType * image, + typename ImageType::PixelType BG, + std::vector & centroids) + { + typedef long LabelType; + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; + typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); + typedef itk::ShapeLabelMapFilter ShapeFilterType; + typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); + imageToLabelFilter->SetBackgroundValue(BG); + imageToLabelFilter->SetInput(image); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->Update(); + typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); + + centroids.clear(); + typename ImageType::PointType dummy; + centroids.push_back(dummy); // label 0 -> no centroid, use dummy point + for(uint i=1; iGetNumberOfLabelObjects()+1; i++) { + centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid()); + } + + for(uint i=1; iGetNumberOfLabelObjects()+1; i++) { + DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes()); + DD(labelMap->GetLabelObject(i)->GetBinaryFlatness()); + DD(labelMap->GetLabelObject(i)->GetRoundness ()); + + // search for the point on the boundary alog PA + + } + + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + PointsUtils::Convert2DTo3D(const PointType2D & p2D, + const ImageType * image, + const int slice, + PointType3D & p3D) + { + IndexType3D index3D; + index3D[0] = index3D[1] = 0; + index3D[2] = image->GetLargestPossibleRegion().GetIndex()[2]+slice; + image->TransformIndexToPhysicalPoint(index3D, p3D); + p3D[0] = p2D[0]; + p3D[1] = p2D[1]; + // p3D[2] = p[2];//(image->GetLargestPossibleRegion().GetIndex()[2]+slice)*image->GetSpacing()[2] + // + image->GetOrigin()[2]; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + PointsUtils::Convert2DMapTo3DList(const MapPoint2DType & map, + const ImageType * image, + VectorPoint3DType & list) + { + typename MapPoint2DType::const_iterator iter = map.begin(); + while (iter != map.end()) { + PointType3D p; + Convert2DTo3D(iter->second, image, iter->first, p); + list.push_back(p); + ++iter; + } + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + PointsUtils::Convert2DListTo3DList(const VectorPoint2DType & p2D, + int slice, + const ImageType * image, + VectorPoint3DType & list) + { + for(uint i=0; i + void + WriteListOfLandmarks(std::vector points, + std::string filename) + { + std::ofstream os; + openFileForWriting(os, filename); + os << "LANDMARKS1" << std::endl; + for(uint i=0; i + typename ImageType::Pointer + Dilate(const ImageType * image, double radiusInMM, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport) + { + typename ImageType::SizeType r; + for(uint i=0; iGetSpacing()[i]); + return Dilate(image, r, BG, FG, extendSupport); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Dilate(const ImageType * image, typename ImageType::PointType radiusInMM, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport) + { + typename ImageType::SizeType r; + for(uint i=0; iGetSpacing()[i]); + return Dilate(image, r, BG, FG, extendSupport); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Dilate(const ImageType * image, typename ImageType::SizeType radius, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG, + bool extendSupport) + { + // Create kernel for dilatation + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(radius); + structuringElement.CreateStructuringElement(); + + typename ImageType::Pointer output; + if (extendSupport) { + typedef itk::ConstantPadImageFilter PadFilterType; + typename PadFilterType::Pointer padFilter = PadFilterType::New(); + padFilter->SetInput(image); + typename ImageType::SizeType lower; + typename ImageType::SizeType upper; + for(uint i=0; i<3; i++) { + lower[i] = upper[i] = 2*(radius[i]+1); + } + padFilter->SetPadLowerBound(lower); + padFilter->SetPadUpperBound(upper); + padFilter->Update(); + output = padFilter->GetOutput(); + } + + // Dilate filter + typedef itk::BinaryDilateImageFilter DilateFilterType; + typename DilateFilterType::Pointer dilateFilter = DilateFilterType::New(); + dilateFilter->SetBackgroundValue(BG); + dilateFilter->SetForegroundValue(FG); + dilateFilter->SetBoundaryToForeground(false); + dilateFilter->SetKernel(structuringElement); + if (extendSupport) dilateFilter->SetInput(output); + else dilateFilter->SetInput(image); + dilateFilter->Update(); + return dilateFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Opening(const ImageType * image, typename ImageType::SizeType radius, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG) + { + // Kernel + typedef itk::BinaryBallStructuringElement KernelType; + KernelType structuringElement; + structuringElement.SetRadius(radius); + structuringElement.CreateStructuringElement(); + + // Filter + typedef itk::BinaryMorphologicalOpeningImageFilter OpeningFilterType; + typename OpeningFilterType::Pointer open = OpeningFilterType::New(); + open->SetInput(image); + open->SetBackgroundValue(BG); + open->SetForegroundValue(FG); + open->SetKernel(structuringElement); + open->Update(); + return open->GetOutput(); + } + //-------------------------------------------------------------------- + + + + //-------------------------------------------------------------------- + template + void ConvertOption(std::string optionName, uint given, + ValueType * values, VectorType & p, + uint dim, bool required) + { + if (required && (given == 0)) { + clitkExceptionMacro("The option --" << optionName << " must be set and have 1 or " + << dim << " values."); + } + if (given == 1) { + for(uint i=0; i to accelerate, start with formula, when change sign -> stop and fill -*/ -template -void -SliceBySliceSetBackgroundFromLineSeparation(typename ImageType::Pointer input, - std::vector & lA, - std::vector & lB, - typename ImageType::PixelType BG, - int mainDirection, - double offsetToKeep) -{ + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + /* + http://www.gamedev.net/community/forums/topic.asp?topic_id=542870 + Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute: + (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax) + This will equal zero if the point C is on the line formed by + points A and B, and will have a different sign depending on the + side. Which side this is depends on the orientation of your (x,y) + coordinates, but you can plug test values for A,B and C into this + formula to determine whether negative values are to the left or to + the right. + => to accelerate, start with formula, when change sign -> stop and fill + + offsetToKeep = is used to determine which side of the line we + keep. The point along the mainDirection but 'offsetToKeep' mm away + is kept. - typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; - SliceIteratorType siter = SliceIteratorType(input, - input->GetLargestPossibleRegion()); - siter.SetFirstDirection(0); - siter.SetSecondDirection(1); - siter.GoToBegin(); - int i=0; - typename ImageType::PointType A; - typename ImageType::PointType B; - typename ImageType::PointType C; - while (!siter.IsAtEnd()) { - // Check that the current slice correspond to the current point - input->TransformIndexToPhysicalPoint(siter.GetIndex(), C); - if (C[2] != lA[i][2]) { - // DD(C); - // DD(lA[i]); - } - else { - // Define A,B,C points - A = lA[i]; - B = lB[i]; - C = A; - C[mainDirection] += offsetToKeep; // I know I must keep this point - double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); - bool isPositive = s<0; - while (!siter.IsAtEndOfSlice()) { - while (!siter.IsAtEndOfLine()) { - // Very slow, I know ... but image should be very small - input->TransformIndexToPhysicalPoint(siter.GetIndex(), C); + */ + template + void + SliceBySliceSetBackgroundFromLineSeparation(ImageType * input, + std::vector & lA, + std::vector & lB, + typename ImageType::PixelType BG, + int mainDirection, + double offsetToKeep) + { + assert((mainDirection==0) || (mainDirection==1)); + typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; + SliceIteratorType siter = SliceIteratorType(input, input->GetLargestPossibleRegion()); + siter.SetFirstDirection(0); + siter.SetSecondDirection(1); + siter.GoToBegin(); + uint i=0; + typename ImageType::PointType A; + typename ImageType::PointType B; + typename ImageType::PointType C; + assert(lA.size() == lB.size()); + while ((iTransformIndexToPhysicalPoint(siter.GetIndex(), C); + if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm + } + else { + // Define A,B,C points + A = lA[i]; + B = lB[i]; + C = A; + // Check that the line is not a point (A=B) + bool p = (A[0] == B[0]) && (A[1] == B[1]); + + if (!p) { + C[mainDirection] += offsetToKeep; // I know I must keep this point double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); - if (s == 0) siter.Set(BG); // on the line, we decide to remove - if (isPositive) { - if (s > 0) siter.Set(BG); - } - else { - if (s < 0) siter.Set(BG); - } - ++siter; + bool isPositive = s<0; + while (!siter.IsAtEndOfSlice()) { + while (!siter.IsAtEndOfLine()) { + // Very slow, I know ... but image should be very small + input->TransformIndexToPhysicalPoint(siter.GetIndex(), C); + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + if (s == 0) siter.Set(BG); // on the line, we decide to remove + if (isPositive) { + if (s > 0) siter.Set(BG); + } + else { + if (s < 0) siter.Set(BG); + } + ++siter; + } + siter.NextLine(); + } // end loop slice + } + + ++i; + } // End of current slice + siter.NextSlice(); + } + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + AndNot(ImageType * input, + const ImageType * object, + typename ImageType::PixelType BG) + { + typename ImageType::Pointer o; + bool resized=false; + if (!clitk::HaveSameSizeAndSpacing(input, object)) { + o = clitk::ResizeImageLike(object, input, BG); + resized = true; + } + + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(input); + if (resized) boolFilter->SetInput2(o); + else boolFilter->SetInput2(object); + boolFilter->SetBackgroundValue1(BG); + boolFilter->SetBackgroundValue2(BG); + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + And(ImageType * input, + const ImageType * object, + typename ImageType::PixelType BG) + { + typename ImageType::Pointer o; + bool resized=false; + if (!clitk::HaveSameSizeAndSpacing(input, object)) { + o = clitk::ResizeImageLike(object, input, BG); + resized = true; + } + + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(input); + if (resized) boolFilter->SetInput2(o); + else boolFilter->SetInput2(object); + boolFilter->SetBackgroundValue1(BG); + boolFilter->SetBackgroundValue2(BG); + boolFilter->SetOperationType(BoolFilterType::And); + boolFilter->Update(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + Or(ImageType * input, + const ImageType * object, + typename ImageType::PixelType BG) + { + typename ImageType::Pointer o; + bool resized=false; + if (!clitk::HaveSameSizeAndSpacing(input, object)) { + o = clitk::ResizeImageLike(object, input, BG); + resized = true; + } + + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(input); + if (resized) boolFilter->SetInput2(o); + else boolFilter->SetInput2(object); + boolFilter->SetBackgroundValue1(BG); + boolFilter->SetBackgroundValue2(BG); + boolFilter->SetOperationType(BoolFilterType::Or); + boolFilter->Update(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Binarize(const ImageType * input, + typename ImageType::PixelType lower, + typename ImageType::PixelType upper, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG) + { + typedef itk::BinaryThresholdImageFilter BinaryThresholdFilterType; + typename BinaryThresholdFilterType::Pointer binarizeFilter = BinaryThresholdFilterType::New(); + binarizeFilter->SetInput(input); + binarizeFilter->InPlaceOff(); + binarizeFilter->SetLowerThreshold(lower); + binarizeFilter->SetUpperThreshold(upper); + binarizeFilter->SetInsideValue(FG); + binarizeFilter->SetOutsideValue(BG); + binarizeFilter->Update(); + return binarizeFilter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + GetMinMaxPointPosition(const ImageType * input, + typename ImageType::PointType & min, + typename ImageType::PointType & max) + { + typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex(); + input->TransformIndexToPhysicalPoint(index, min); + index = index+input->GetLargestPossibleRegion().GetSize(); + input->TransformIndexToPhysicalPoint(index, max); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::PointType + FindExtremaPointInAGivenLine(const ImageType * input, + int dimension, + bool inverse, + typename ImageType::PointType p, + typename ImageType::PixelType BG, + double distanceMax) + { + // Which direction ? Increasing or decreasing. + int d=1; + if (inverse) d=-1; + + // Transform to pixel index + typename ImageType::IndexType index; + input->TransformPhysicalPointToIndex(p, index); + + // Loop while inside the mask; + while (input->GetPixel(index) != BG) { + index[dimension] += d; + } + + // Transform back to Physical Units + typename ImageType::PointType result; + input->TransformIndexToPhysicalPoint(index, result); + + // Check that is is not too far away + double distance = p.EuclideanDistanceTo(result); + if (distance > distanceMax) { + result = p; // Get back to initial value + } + + return result; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + bool + IsOnTheSameLineSide(PointType C, PointType A, PointType B, PointType like) + { + // Look at the position of point 'like' according to the AB line + double s = (B[0] - A[0]) * (like[1] - A[1]) - (B[1] - A[1]) * (like[0] - A[0]); + bool negative = s<0; + + // Look the C position + s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + + if (negative && (s<=0)) return true; + if (!negative && (s>=0)) return true; + return false; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + /* Consider an input object, for each slice, find the extrema + position according to a given direction and build a line segment + passing throught this point in a given direction. Output is a + vector of line (from point A to B), for each slice; + */ + template + void + SliceBySliceBuildLineSegmentAccordingToExtremaPosition(const ImageType * input, + typename ImageType::PixelType BG, + int sliceDimension, + int extremaDirection, + bool extremaOppositeFlag, + int lineDirection, + double margin, + std::vector & A, + std::vector & B) + { + // Type of a slice + typedef typename itk::Image SliceType; + + // Build the list of slices + std::vector slices; + clitk::ExtractSlices(input, sliceDimension, slices); + + // Build the list of 2D points + std::map position2D; + for(uint i=0; i(slices[i], BG, + extremaDirection, extremaOppositeFlag, p); + if (found) { + position2D[i] = p; + } + } + + // Convert 2D points in slice into 3D points + clitk::PointsUtils::Convert2DMapTo3DList(position2D, input, A); + + // Create additional point just right to the previous ones, on the + // given lineDirection, in order to create a horizontal/vertical line. + for(uint i=0; i + typename ImageType::Pointer + SliceBySliceKeepMainCCL(const ImageType * input, + typename ImageType::PixelType BG, + typename ImageType::PixelType FG) { + + // Extract slices + const int d = ImageType::ImageDimension-1; + typedef typename itk::Image SliceType; + std::vector slices; + clitk::ExtractSlices(input, d, slices); + + // Labelize and keep the main one + std::vector o; + for(uint i=0; i(slices[i], BG, false, 1)); + o[i] = clitk::KeepLabels(o[i], BG, FG, 1, 1, true); + } + + // Join slices + typename ImageType::Pointer output; + output = clitk::JoinSlices(o, input, d); + return output; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + Clone(const ImageType * input) { + typedef itk::ImageDuplicator DuplicatorType; + typename DuplicatorType::Pointer duplicator = DuplicatorType::New(); + duplicator->SetInputImage(input); + duplicator->Update(); + return duplicator->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + /* Consider an input object, start at A, for each slice (dim1): + - compute the intersection between the AB line and the current slice + - remove what is at lower or greater according to dim2 of this point + - stop at B + */ + template + typename ImageType::Pointer + SliceBySliceSetBackgroundFromSingleLine(const ImageType * input, + typename ImageType::PixelType BG, + typename ImageType::PointType & A, + typename ImageType::PointType & B, + int dim1, int dim2, bool removeLowerPartFlag) + + { + // Extract slices + typedef typename itk::Image SliceType; + typedef typename SliceType::Pointer SlicePointer; + std::vector slices; + clitk::ExtractSlices(input, dim1, slices); + + // Start at slice that contains A, and stop at B + typename ImageType::IndexType Ap; + typename ImageType::IndexType Bp; + input->TransformPhysicalPointToIndex(A, Ap); + input->TransformPhysicalPointToIndex(B, Bp); + + // Determine slice largest region + typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion(); + typename SliceType::SizeType size = region.GetSize(); + typename SliceType::IndexType index = region.GetIndex(); + + // Line slope + double a = (Bp[dim2]-Ap[dim2])/(Bp[dim1]-Ap[dim1]); + double b = Ap[dim2]; + + // Loop from slice A to slice B + for(uint i=0; i<(Bp[dim1]-Ap[dim1]); i++) { + // Compute intersection between line AB and current slice for the dim2 + double p = a*i+b; + // Change region (lower than dim2) + if (removeLowerPartFlag) { + size[dim2] = p-Ap[dim2]; + } + else { + size[dim2] = slices[0]->GetLargestPossibleRegion().GetSize()[dim2]-p; + index[dim2] = p; + } + region.SetSize(size); + region.SetIndex(index); + // Fill region with BG (simple region iterator) + FillRegionWithValue(slices[i+Ap[dim1]], BG, region); + /* + typedef itk::ImageRegionIterator IteratorType; + IteratorType iter(slices[i+Ap[dim1]], region); + iter.GoToBegin(); + while (!iter.IsAtEnd()) { + iter.Set(BG); + ++iter; + } + */ + // Loop + } + + // Merge slices + typename ImageType::Pointer output; + output = clitk::JoinSlices(slices, input, dim1); + return output; + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + /* Consider an input object, slice by slice, use the point A and set + pixel to BG according to their position relatively to A + */ + template + typename ImageType::Pointer + SliceBySliceSetBackgroundFromPoints(const ImageType * input, + typename ImageType::PixelType BG, + int sliceDim, + std::vector & A, + bool removeGreaterThanXFlag, + bool removeGreaterThanYFlag) + + { + // Extract slices + typedef typename itk::Image SliceType; + typedef typename SliceType::Pointer SlicePointer; + std::vector slices; + clitk::ExtractSlices(input, sliceDim, slices); + + // Start at slice that contains A + typename ImageType::IndexType Ap; + + // Determine slice largest region + typename SliceType::RegionType region = slices[0]->GetLargestPossibleRegion(); + typename SliceType::SizeType size = region.GetSize(); + typename SliceType::IndexType index = region.GetIndex(); + + // Loop from slice A to slice B + for(uint i=0; iTransformPhysicalPointToIndex(A[i], Ap); + uint sliceIndex = Ap[2] - input->GetLargestPossibleRegion().GetIndex()[2]; + if ((sliceIndex < 0) || (sliceIndex >= slices.size())) { + continue; // do not consider this slice + } + + // Compute region for BG + if (removeGreaterThanXFlag) { + index[0] = Ap[0]; + size[0] = region.GetSize()[0]-(index[0]-region.GetIndex()[0]); + } + else { + index[0] = region.GetIndex()[0]; + size[0] = Ap[0] - index[0]; + } + + if (removeGreaterThanYFlag) { + index[1] = Ap[1]; + size[1] = region.GetSize()[1]-(index[1]-region.GetIndex()[1]); + } + else { + index[1] = region.GetIndex()[1]; + size[1] = Ap[1] - index[1]; + } + + // Set region + region.SetSize(size); + region.SetIndex(index); + + // Fill region with BG (simple region iterator) + FillRegionWithValue(slices[sliceIndex], BG, region); + // Loop + } + + // Merge slices + typename ImageType::Pointer output; + output = clitk::JoinSlices(slices, input, sliceDim); + return output; + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region) + { + typedef itk::ImageRegionIterator IteratorType; + IteratorType iter(input, region); + iter.GoToBegin(); + while (!iter.IsAtEnd()) { + iter.Set(value); + ++iter; + } + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + GetMinMaxBoundary(ImageType * input, typename ImageType::PointType & min, + typename ImageType::PointType & max) + { + typedef typename ImageType::PointType PointType; + typedef typename ImageType::IndexType IndexType; + IndexType min_i, max_i; + min_i = input->GetLargestPossibleRegion().GetIndex(); + for(uint i=0; iGetLargestPossibleRegion().GetSize()[i] + min_i[i]; + input->TransformIndexToPhysicalPoint(min_i, min); + input->TransformIndexToPhysicalPoint(max_i, max); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename itk::Image::Pointer + DistanceMap(const ImageType * input, typename ImageType::PixelType BG)//, + // typename itk::Image::Pointer dmap) + { + typedef itk::Image FloatImageType; + typedef itk::SignedMaurerDistanceMapImageFilter DistanceMapFilterType; + typename DistanceMapFilterType::Pointer filter = DistanceMapFilterType::New(); + filter->SetInput(input); + filter->SetUseImageSpacing(true); + filter->SquaredDistanceOff(); + filter->SetBackgroundValue(BG); + filter->Update(); + return filter->GetOutput(); + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + void + SliceBySliceBuildLineSegmentAccordingToMinimalDistanceBetweenStructures(const ImageType * S1, + const ImageType * S2, + typename ImageType::PixelType BG, + int sliceDimension, + std::vector & A, + std::vector & B) + { + // Extract slices + typedef typename itk::Image SliceType; + typedef typename SliceType::Pointer SlicePointer; + std::vector slices_s1; + std::vector slices_s2; + clitk::ExtractSlices(S1, sliceDimension, slices_s1); + clitk::ExtractSlices(S2, sliceDimension, slices_s2); + + assert(slices_s1.size() == slices_s2.size()); + + // Prepare dmap + typedef itk::Image FloatImageType; + typedef itk::SignedMaurerDistanceMapImageFilter DistanceMapFilterType; + std::vector dmaps1; + std::vector dmaps2; + typename FloatImageType::Pointer dmap; + + // loop on slices + for(uint i=0; i(slices_s1[i], BG); + dmaps1.push_back(dmap); + writeImage(dmap, "dmap1.mha"); + // Compute dmap for S2 + dmap = clitk::DistanceMap(slices_s2[i], BG); + dmaps2.push_back(dmap); + writeImage(dmap, "dmap2.mha"); + + // Look in S2 for the point the closest to S1 + typename SliceType::PointType p = ComputeClosestPoint(slices_s1[i], dmaps2[i], BG); + typename ImageType::PointType p3D; + clitk::PointsUtils::Convert2DTo3D(p, S1, i, p3D); + A.push_back(p3D); + + // Look in S2 for the point the closest to S1 + p = ComputeClosestPoint(slices_s2[i], dmaps1[i], BG); + clitk::PointsUtils::Convert2DTo3D(p, S2, i, p3D); + B.push_back(p3D); + + } + + // Debug dmap + /* + typedef itk::Image FT; + FT::Pointer f = FT::New(); + typename FT::Pointer d1 = clitk::JoinSlices(dmaps1, S1, 2); + typename FT::Pointer d2 = clitk::JoinSlices(dmaps2, S2, 2); + writeImage(d1, "d1.mha"); + writeImage(d2, "d2.mha"); + */ + } + //-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- + template + typename ImageType::PointType + ComputeClosestPoint(const ImageType * input, + const itk::Image * dmap, + typename ImageType::PixelType & BG) + { + // Loop dmap + S2, if FG, get min + typedef itk::Image FloatImageType; + typedef itk::ImageRegionConstIteratorWithIndex ImageIteratorType; + typedef itk::ImageRegionConstIterator DMapIteratorType; + ImageIteratorType iter1(input, input->GetLargestPossibleRegion()); + DMapIteratorType iter2(dmap, dmap->GetLargestPossibleRegion()); + + iter1.GoToBegin(); + iter2.GoToBegin(); + double dmin = 100000.0; + typename ImageType::IndexType indexmin; + while (!iter1.IsAtEnd()) { + if (iter1.Get() != BG) { + double d = iter2.Get(); + if (dTransformIndexToPhysicalPoint(indexmin, p); + return p; } -} -//-------------------------------------------------------------------- -} + //-------------------------------------------------------------------- + + + + +} // end of namespace +