+ //--------------------------------------------------------------------
+ /* 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<class ImageType>
+ 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<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
+ typedef typename SliceType::Pointer SlicePointer;
+ std::vector<SlicePointer> slices;
+ clitk::ExtractSlices<ImageType>(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<SliceType>(slices[i+Ap[dim1]], BG, region);
+ /*
+ typedef itk::ImageRegionIterator<SliceType> 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<ImageType>(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<class ImageType>
+ typename ImageType::Pointer
+ SliceBySliceSetBackgroundFromPoints(const ImageType * input,
+ typename ImageType::PixelType BG,
+ int sliceDim,
+ std::vector<typename ImageType::PointType> & A,
+ bool removeGreaterThanXFlag,
+ bool removeGreaterThanYFlag)
+
+ {
+ // Extract slices
+ typedef typename itk::Image<typename ImageType::PixelType, ImageType::ImageDimension-1> SliceType;
+ typedef typename SliceType::Pointer SlicePointer;
+ std::vector<SlicePointer> slices;
+ clitk::ExtractSlices<ImageType>(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; i<A.size(); i++) {
+ input->TransformPhysicalPointToIndex(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<SliceType>(slices[sliceIndex], BG, region);
+ // Loop
+ }
+
+ // Merge slices
+ typename ImageType::Pointer output;
+ output = clitk::JoinSlices<ImageType>(slices, input, sliceDim);
+ return output;
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template<class ImageType>
+ void
+ FillRegionWithValue(ImageType * input, typename ImageType::PixelType value, typename ImageType::RegionType & region)
+ {
+ typedef itk::ImageRegionIterator<ImageType> IteratorType;
+ IteratorType iter(input, region);
+ iter.GoToBegin();
+ while (!iter.IsAtEnd()) {
+ iter.Set(value);
+ ++iter;
+ }
+ }
+ //--------------------------------------------------------------------
+
+
+ //--------------------------------------------------------------------
+ template<class ImageType>
+ 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; i<ImageType::ImageDimension; i++)
+ max_i[i] = input->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
+ input->TransformIndexToPhysicalPoint(min_i, min);
+ input->TransformIndexToPhysicalPoint(max_i, max);
+ }
+ //--------------------------------------------------------------------
+