- //--------------------------------------------------------------------
- // Step 7 : Slice by Slice to optimize posterior part
- // Warning slice does not necesseraly correspond between 'output' and 'bones'
- typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
- typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
- typedef typename ExtractSliceFilterType::SliceType SliceType;
- std::vector<typename SliceType::Pointer> mSlices;
- if (GetUseBones()) {
- StartNewStep("Rafine posterior part according to vertebral body");
- extractSliceFilter->SetInput(bones_post);
- extractSliceFilter->SetDirection(2);
- extractSliceFilter->Update();
- std::vector<double> mVertebralAntPositionBySlice;
- extractSliceFilter->GetOutputSlices(mSlices);
- for(unsigned int i=0; i<mSlices.size(); i++) {
- mSlices[i] = Labelize<SliceType>(mSlices[i], 0, true, 10);
- mSlices[i] = KeepLabels<SliceType>(mSlices[i],
- GetBackgroundValue(),
- GetForegroundValue(), 1, 2, true); // keep two first
- // Find most anterior point (start of the vertebral)
- typename itk::ImageRegionIteratorWithIndex<SliceType>
- iter(mSlices[i], mSlices[i]->GetLargestPossibleRegion());
- iter.GoToBegin();
- bool stop = false;
- while (!stop) {
- if (iter.Get() != GetBackgroundValue())
- stop = true; // not foreground because we keep two main label
- ++iter;
- if (iter.IsAtEnd()) stop = true;
- }
- if (!iter.IsAtEnd()) {
- typename SliceType::PointType p;
- mSlices[i]->TransformIndexToPhysicalPoint(iter.GetIndex(),p);
- mVertebralAntPositionBySlice.push_back(p[1]);
- }
- else {
- mVertebralAntPositionBySlice.push_back(bones_post->GetOrigin()[1]+(bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1]));
- DD(mVertebralAntPositionBySlice.back());
- DD("ERROR ?? NO FG in bones here ?");
- }
- }
-
- // Cut Post position slice by slice
- {
- MaskImageRegionType region;
- MaskImageSizeType size;
- MaskImageIndexType start;
- size[2] = 1;
- start[0] = output->GetLargestPossibleRegion().GetIndex()[0];
- for(unsigned int i=0; i<mSlices.size(); i++) {
- // Compute index
- MaskImagePointType point;
- point[0] = 0;
-
- //TODO 10 mm OPTION
-
- point[1] = mVertebralAntPositionBySlice[i]+GetDistanceMaxToAnteriorPartOfTheSpine();// ADD ONE CM
- point[2] = bones_post->GetOrigin()[2]+(bones_post->GetLargestPossibleRegion().GetIndex()[2]+i)*bones_post->GetSpacing()[2];
- MaskImageIndexType index;
- output->TransformPhysicalPointToIndex(point, index);
- // Compute region
- start[2] = index[2];
- start[1] = output->GetLargestPossibleRegion().GetIndex()[1]+index[1];
- size[0] = output->GetLargestPossibleRegion().GetSize()[0];
- size[1] = output->GetLargestPossibleRegion().GetSize()[1]-start[1];
- region.SetSize(size);
- region.SetIndex(start);
- // Fill Region
- if (output->GetLargestPossibleRegion().IsInside(start)) {
- itk::ImageRegionIteratorWithIndex<MaskImageType> it(output, region);
- it.GoToBegin();
- while (!it.IsAtEnd()) {
- it.Set(GetBackgroundValue());
- ++it;
- }
- }
- }
- }
- this->template StopCurrentStep<MaskImageType>(output);
- }