X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=itk%2FclitkSegmentationUtils.txx;h=75e5ef35f982b2d7e33fff174af5ae7b77c8ad40;hb=5f801bf0b07486889123e941d6913d4369dfc86f;hp=adca988098ad8c78fee5cfb132b890fb0bef0c09;hpb=42e050104aebe7830f2111d096f0dd6673987d1e;p=clitk.git diff --git a/itk/clitkSegmentationUtils.txx b/itk/clitkSegmentationUtils.txx index adca988..75e5ef3 100644 --- a/itk/clitkSegmentationUtils.txx +++ b/itk/clitkSegmentationUtils.txx @@ -38,89 +38,6 @@ namespace clitk { - //-------------------------------------------------------------------- - template - void ComputeBBFromImageRegion(const ImageType * 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(const ImageType * 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]); - } - - // Create region - region.SetIndex(regionStart); - region.SetSize(regionSize); - } - //-------------------------------------------------------------------- - //-------------------------------------------------------------------- template typename ImageType::Pointer @@ -313,19 +230,37 @@ namespace clitk { //-------------------------------------------------------------------- - template - typename ImageType::Pointer - ResizeImageLike(const ImageType * input, - const itk::ImageBase * like, - typename ImageType::PixelType backgroundValue) + 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 CropLikeImageFilter CropFilterType; - typename CropFilterType::Pointer cropFilter = CropFilterType::New(); - cropFilter->SetInput(input); - cropFilter->SetCropLikeImage(like); - cropFilter->SetBackgroundValue(backgroundValue); - cropFilter->Update(); - return cropFilter->GetOutput(); + 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(); } //-------------------------------------------------------------------- @@ -337,13 +272,14 @@ namespace clitk { const MaskImageType * object, int direction, double threshold, - std::string orientation, + double angle, + bool inverseflag, bool uniqueConnectedComponent, double spacing, bool autocropFlag, bool singleObjectCCL) { - typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typedef clitk::SliceBySliceRelativePositionFilter SliceRelPosFilterType; typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); sliceRelPosFilter->VerboseStepFlagOff(); sliceRelPosFilter->WriteStepFlagOff(); @@ -351,13 +287,14 @@ namespace clitk { sliceRelPosFilter->SetInputObject(object); sliceRelPosFilter->SetDirection(direction); sliceRelPosFilter->SetFuzzyThreshold(threshold); - sliceRelPosFilter->AddOrientationTypeString(orientation); + // 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->SetInverseOrientationFlag(inverseflag); sliceRelPosFilter->SetAutoCropFlag(autocropFlag); sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); sliceRelPosFilter->Update(); @@ -417,7 +354,7 @@ namespace clitk { ++iter; } if (!found) return false; - input->TransformIndexToPhysicalPoint(max, point); + input->TransformIndexToPhysicalPoint(max, point); // half of the pixel return true; } //-------------------------------------------------------------------- @@ -445,9 +382,11 @@ namespace clitk { int dim, double max, bool autoCrop, typename ImageType::PixelType BG) { - typename ImageType::PointType p; + typename ImageType::PointType p; + image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+ image->GetLargestPossibleRegion().GetSize(), p); + return CropImageAlongOneAxis(image, dim, max, p[dim], autoCrop, BG); } //-------------------------------------------------------------------- @@ -463,14 +402,24 @@ namespace clitk { // Compute region size typename ImageType::RegionType region; typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize(); - typename ImageType::PointType p = image->GetOrigin(); - p[dim] = min; + + // Starting index + typename ImageType::PointType p = image->GetOrigin(); // not at pixel center ! + if (min > p[dim]) p[dim] = min; // Check if not outside the image typename ImageType::IndexType start; image->TransformPhysicalPointToIndex(p, start); - p[dim] = max; + + // Size of the region + // -1 because last point is size -1 + double m = image->GetOrigin()[dim] + (size[dim]-1)*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]); + size[dim] = abs(end[dim]-start[dim])+1;// +1 because we want to include the point. + + // Set region region.SetIndex(start); region.SetSize(size); @@ -819,12 +768,26 @@ namespace clitk { std::vector & lB, typename ImageType::PixelType BG, int mainDirection, - double offsetToKeep) + double offsetToKeep, + bool keepIfEqual) { assert((mainDirection==0) || (mainDirection==1)); + typename ImageType::PointType offset; + offset[0] = offset[1] = offset[2] = 0.0; + offset[mainDirection] = offsetToKeep; + SliceBySliceSetBackgroundFromLineSeparation_pt(input, lA, lB, BG, offset, keepIfEqual); + } + template + void + SliceBySliceSetBackgroundFromLineSeparation_pt(ImageType * input, + std::vector & lA, + std::vector & lB, + typename ImageType::PixelType BG, + typename ImageType::PointType offsetToKeep, + bool keepIfEqual) + { typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; - SliceIteratorType siter = SliceIteratorType(input, - input->GetLargestPossibleRegion()); + SliceIteratorType siter = SliceIteratorType(input, input->GetLargestPossibleRegion()); siter.SetFirstDirection(0); siter.SetSecondDirection(1); siter.GoToBegin(); @@ -837,31 +800,39 @@ namespace clitk { // Check that the current slice correspond to the current point input->TransformIndexToPhysicalPoint(siter.GetIndex(), C); if ((fabs(C[2] - lA[i][2]))>0.01) { // is !equal with a tolerance of 0.01 mm + // FIXME : if not the same slices, should advance i or slice (not done yet) + // clitkExceptionMacro("Error list of point and slice do not start at the same location"); } 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 + //C[mainDirection] += offsetToKeep; // I know I must keep this point + C[0] += offsetToKeep[0]; + C[1] += offsetToKeep[1]; + //C[2] += offsetToKeep[2]; 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); - 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); + if (siter.Get() != BG) { // do only if not BG + // 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) { + if (!keepIfEqual) 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; } @@ -1380,11 +1351,11 @@ namespace clitk { // Compute dmap for S1 *TO PUT IN FONCTION* dmap = clitk::DistanceMap(slices_s1[i], BG); dmaps1.push_back(dmap); - writeImage(dmap, "dmap1.mha"); + //writeImage(dmap, "dmap1.mha"); // Compute dmap for S2 dmap = clitk::DistanceMap(slices_s2[i], BG); dmaps2.push_back(dmap); - writeImage(dmap, "dmap2.mha"); + //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); @@ -1399,14 +1370,14 @@ namespace clitk { } - // 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"); + // 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"); */ } //-------------------------------------------------------------------- @@ -1430,6 +1401,7 @@ namespace clitk { iter2.GoToBegin(); double dmin = 100000.0; typename ImageType::IndexType indexmin; + indexmin.Fill(0); while (!iter1.IsAtEnd()) { if (iter1.Get() != BG) { double d = iter2.Get(); @@ -1450,6 +1422,30 @@ namespace clitk { //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + template + typename ImageType::Pointer + RemoveNegativeIndexFromRegion(ImageType * input) { + typedef itk::ChangeInformationImageFilter< ImageType > InfoFilterType; + typename InfoFilterType::Pointer indexChangeFilter = InfoFilterType::New(); + indexChangeFilter->ChangeRegionOn(); + // The next line is commented because not exist in itk 3 + // typename InfoFilterType::OutputImageOffsetValueType indexShift[3]; + long indexShift[3]; + typename ImageType::IndexType index = input->GetLargestPossibleRegion().GetIndex(); + for(uint i=0;iGetOrigin()[i] - indexShift[i]*input->GetSpacing()[i]; + indexChangeFilter->SetOutputOffset( indexShift ); + indexChangeFilter->SetInput(input); + indexChangeFilter->SetOutputOrigin(origin); + indexChangeFilter->ChangeOriginOn(); + indexChangeFilter->Update(); + return indexChangeFilter->GetOutput(); + } + //-------------------------------------------------------------------- } // end of namespace