+ std::vector<typename MaskSliceType::Pointer> bronchi_slices;
+ clitk::ExtractSlices<MaskImageType>(crop_trachea, 2, bronchi_slices);
+
+ // List of midpoints
+ std::vector<MaskImagePointType> midpoints;
+
+ // Add mid points below carina, from foot to head
+ for(uint i=0; i<m_LeftMostInRightBronchus.size(); i++) {
+ MaskImagePointType p;
+ p[0] = (m_LeftMostInRightBronchus[i][0]+m_RightMostInLeftBronchus[i][0])/2.0;
+ p[1] = (m_LeftMostInRightBronchus[i][1]+m_RightMostInLeftBronchus[i][1])/2.0;
+ p[2] = m_LeftMostInRightBronchus[i][2];
+ midpoints.push_back(p);
+ }
+
+ // Crop the trachea above the carina. Find largest Z position
+ MaskImageIndexType p = Trachea->GetLargestPossibleRegion().GetIndex();
+ for(uint i=0; i<3; i++) {
+ p[i] += Trachea->GetLargestPossibleRegion().GetSize()[i];
+ }
+ MaskImagePointType q;
+ Trachea->TransformIndexToPhysicalPoint(p, q);
+ double maxZ = q[2];
+ double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2);
+ MaskImagePointer m_above_carina =
+ clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2,
+ m_CarinaZ,
+ maxZ, true,
+ GetBackgroundValue());
+ writeImage<MaskImageType>(m_above_carina, "above.mhd");
+
+ // Extract all the slices
+ std::vector<typename MaskSliceType::Pointer> trachea_slices;
+ clitk::ExtractSlices<MaskImageType>(m_above_carina, 2, trachea_slices);
+
+ // Find centroid of the trachea in each slice
+ std::vector<MaskImagePointType> points;
+ typedef typename MaskSliceType::PointType SlicePointType;
+ SlicePointType previous;
+ // Start from patient top (head)
+ for(uint i=0; i<trachea_slices.size(); i++) {
+ // Labelize
+ trachea_slices[i] = Labelize<MaskSliceType>(trachea_slices[i], GetBackgroundValue(), true, 10);
+ // Get centroid
+ std::vector<SlicePointType> c;
+ clitk::ComputeCentroids<MaskSliceType>(trachea_slices[i], GetBackgroundValue(), c);
+ // Keep first one (first connected component label) and convert to 3D
+ MaskImagePointType p;
+ p[2] = m_above_carina->GetOrigin()[2] + i*m_above_carina->GetSpacing()[2];
+ p[0] = c[1][0];
+ p[1] = c[1][1];
+ midpoints.push_back(p);
+ }
+
+ // DEBUG POINTS
+ std::ofstream osp;
+ openFileForWriting(osp, "mp.txt");
+ osp << "LANDMARKS1" << std::endl;
+ for(uint i=0; i<midpoints.size(); i++) {
+ osp << i << " " << midpoints[i][0] << " "
+ << midpoints[i][1] << " " << midpoints[i][2] << " 0 0 " << std::endl;
+ }
+ osp.close();
+
+ // Create and allocate left and right part of the support (we create
+ // images because we will then crop them)
+ m_LeftSupport = clitk::NewImageLike<MaskImageType>(m_Working_Support);
+ m_RightSupport = clitk::NewImageLike<MaskImageType>(m_Working_Support);
+
+ // Loop on current support, slice by slice
+ typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
+ SliceIteratorType iter = SliceIteratorType(m_Working_Support,
+ m_Working_Support->GetLargestPossibleRegion());
+ SliceIteratorType iterL = SliceIteratorType(m_LeftSupport,
+ m_LeftSupport->GetLargestPossibleRegion());
+ SliceIteratorType iterR = SliceIteratorType(m_RightSupport,
+ m_RightSupport->GetLargestPossibleRegion());
+ iter.SetFirstDirection(0);
+ iter.SetSecondDirection(1);
+ iter.GoToBegin();
+ iterL.SetFirstDirection(0);
+ iterL.SetSecondDirection(1);
+ iterL.GoToBegin();
+ iterR.SetFirstDirection(0);
+ iterR.SetSecondDirection(1);
+ iterR.GoToBegin();
+
+ int slice=0;
+ // Assert starting of image has the same Z than first midpoints
+ while (midpoints[slice][2] != m_Working_Support->GetOrigin()[2]) {
+ slice++;
+ if ((uint)slice >= midpoints.size()) {
+ clitkExceptionMacro("Bug while searching for first midpoint to use");
+ }
+ }
+
+ // Start loop
+ while (!iter.IsAtEnd()) {
+ ImageIndexType index;
+ m_Working_Support->TransformPhysicalPointToIndex(midpoints[slice], index);
+ while (!iter.IsAtEndOfSlice()) {
+ // convert into index
+ while (!iter.IsAtEndOfLine()) {
+ // if indexcourant <index this is Right part. Left otherwise
+ if (iter.GetIndex()[0]<index[0]) {
+ iterR.Set(iter.Get());
+ iterL.Set(GetBackgroundValue());
+ }
+ else {
+ iterL.Set(iter.Get());
+ iterR.Set(GetBackgroundValue());
+ }
+ ++iter;
+ ++iterL;
+ ++iterR;
+ }
+ iter.NextLine();
+ iterL.NextLine();
+ iterR.NextLine();
+
+ }
+ iter.NextSlice();
+ iterL.NextSlice();
+ iterR.NextSlice();
+ ++slice;
+ }
+
+ // Auto crop !
+ m_LeftSupport = clitk::AutoCrop<MaskImageType>(m_LeftSupport, GetBackgroundValue());
+ m_RightSupport = clitk::AutoCrop<MaskImageType>(m_RightSupport, GetBackgroundValue());
+ writeImage<MaskImageType>(m_LeftSupport, "lsac.mhd");
+ writeImage<MaskImageType>(m_RightSupport, "rsac.mhd");
+
+ StopCurrentStep<MaskImageType>(m_RightSupport);
+
+ // -------------------------------------------------------------------------
+ // Constraint at right from the SVC. Must be after MidTrachea limits
+ // (if not, bronchus can be split, midposition is false)
+ StartNewStep("[Station 4R] R limits with SVC ");
+ MaskImagePointer svc = GetAFDB()->template GetImage<MaskImageType>("SVC");
+ typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
+ typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
+ /* relPosFilter->SetVerboseStep(false);
+ relPosFilter->SetInput(m_RightSupport); // only right here ...
+ relPosFilter->SetInputObject(svc);
+ relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
+ relPosFilter->SetIntermediateSpacing(2);
+ relPosFilter->SetFuzzyThreshold(0.3);
+ relPosFilter->Update();
+ m_RightSupport = relPosFilter->GetOutput();
+ */