- // 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);
- }
-
- //--------------------------------------------------------------------
- // Step 8: Trial segmentation KMeans
- if (0) {
- StartNewStep("K means");
- // Take input, crop like current mask
- typedef CropLikeImageFilter<ImageType> CropLikeFilterType;
- typename CropLikeFilterType::Pointer cropLikeFilter = CropLikeFilterType::New();
- cropLikeFilter->SetInput(input);
- cropLikeFilter->SetCropLikeImage(output);
- cropLikeFilter->Update();
- ImagePointer working_input = cropLikeFilter->GetOutput();
- writeImage<ImageType>(working_input, "crop-input.mhd");
- // Set bG at -1000
- working_input = clitk::SetBackground<ImageType, MaskImageType>(working_input, output, GetBackgroundValue(), -1000);
- writeImage<ImageType>(working_input, "crop-input2.mhd");
- // Kmeans
- typedef itk::ScalarImageKmeansImageFilter<ImageType> KMeansFilterType;
- typename KMeansFilterType::Pointer kmeansFilter = KMeansFilterType::New();
- kmeansFilter->SetInput(working_input);
- // const unsigned int numberOfInitialClasses = 3;
- // const unsigned int useNonContiguousLabels = 0;
- kmeansFilter->AddClassWithInitialMean(-1000);
- kmeansFilter->AddClassWithInitialMean(30);
- kmeansFilter->AddClassWithInitialMean(-40); // ==> I want this one
- DD("Go!");
- kmeansFilter->Update();
- DD("End");
- typename KMeansFilterType::ParametersType estimatedMeans = kmeansFilter->GetFinalMeans();
- const unsigned int numberOfClasses = estimatedMeans.Size();
- for ( unsigned int i = 0 ; i < numberOfClasses ; ++i ) {
- std::cout << "cluster[" << i << "] ";
- std::cout << " estimated mean : " << estimatedMeans[i] << std::endl;
- }
- MaskImageType::Pointer kmeans = kmeansFilter->GetOutput();
- kmeans = clitk::SetBackground<MaskImageType, MaskImageType>(kmeans, kmeans,
- 1, GetBackgroundValue());
- writeImage<MaskImageType>(kmeans, "kmeans.mhd");
- // Get final results, and remove from current mask
- boolFilter = BoolFilterType::New();
- boolFilter->InPlaceOn();
- boolFilter->SetInput1(output);
- boolFilter->SetInput2(kmeans);
- boolFilter->SetOperationType(BoolFilterType::And);
- boolFilter->Update();
- output = boolFilter->GetOutput();
- writeImage<MaskImageType>(output, "out-kmean.mhd");
- this->template StopCurrentStep<MaskImageType>(output);
-
- // TODO -> FillMASK ?
- // comment speed ? mask ? 2 class ?
-
-
- //TODO
- // Confidence connected ?
-
- }
-
- //--------------------------------------------------------------------
- // Step 8: Lower limits from lung (need separate lung ?)
- if (1) {
- // StartNewStep("Trial : minus segmented struct");
- // MaskImagePointer heart = GetAFDB()->template GetImage <MaskImageType>("heart");
- // boolFilter = BoolFilterType::New();
- // boolFilter->InPlaceOn();
- // boolFilter->SetInput1(output);
- // boolFilter->SetInput2(heart);
- // boolFilter->SetOperationType(BoolFilterType::AndNot);
- // boolFilter->Update();
- // output = boolFilter->GetOutput(); // not needed because InPlace
-
- // Not below the heart
- // relPosFilter = RelPosFilterType::New();
- // relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
- // relPosFilter->VerboseStepOff();
- // relPosFilter->WriteStepOff();
- // relPosFilter->SetInput(output);
- // relPosFilter->SetInputObject(heart);
- // relPosFilter->SetOrientationType(RelPosFilterType::SupTo);
- // relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing());
- // relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold3());
- // relPosFilter->Update();
- // output = relPosFilter->GetOutput();
- }