]> Creatis software - clitk.git/commitdiff
first version of Image2DicomRT (dos not work)
authorDavid Sarrut <david.sarrut@creatis.insa-lyon.fr>
Thu, 28 Jul 2011 10:58:40 +0000 (12:58 +0200)
committerDavid Sarrut <david.sarrut@creatis.insa-lyon.fr>
Thu, 28 Jul 2011 10:58:40 +0000 (12:58 +0200)
12 files changed:
common/clitkDicomRT_Contour.cxx
common/clitkDicomRT_Contour.h
common/clitkDicomRT_ROI.cxx
common/clitkDicomRT_ROI.h
common/clitkDicomRT_StructureSet.cxx
common/clitkImage2DicomRTStructFilter.txx
segmentation/clitkExtractLymphStation_3A.txx
segmentation/clitkExtractLymphStation_3P.txx
segmentation/clitkExtractLymphStation_8.txx
segmentation/clitkExtractLymphStationsFilter.h
segmentation/clitkExtractLymphStationsFilter.txx
tools/clitkImage2DicomRTStruct.cxx

index a4a69dd600df11350d337434e36937e0316d5895..861717002ea65fa76b0a19ae12d62d19ea42500e 100644 (file)
@@ -145,6 +145,10 @@ bool clitk::DicomRT_Contour::Read(gdcm::Item * item)
   at.SetFromDataElement( contourdata );
   const double* points = at.GetValues();
 <<<<<<< Updated upstream
+<<<<<<< variant A
+>>>>>>> variant B
+  //  unsigned int npts = at.GetNumberOfValues() / 3;
+======= end
 =======
   //  unsigned int npts = at.GetNumberOfValues() / 3;
 >>>>>>> Stashed changes
index 0484da6c3c0301517089e37874f3f30c0f46b938..91d75d6ffb411be64962f2091527a1e169428372 100644 (file)
@@ -42,17 +42,23 @@ public:
   itkNewMacro(Self);
 
   void Print(std::ostream & os = std::cout) const;
+
 #if GDCM_MAJOR_VERSION == 2
-  bool Read(gdcm::Item const & item);
+  bool Read(gdcm::Item * item);
+  void UpdateDicomItem();
 #else
   bool Read(gdcm::SQItem * item);
 #endif
+
   vtkPolyData * GetMesh();
+  void SetMesh(vtkPolyData * mesh);
   vtkPoints * GetPoints() {return mData;}
   double GetZ() const {return mZ;}
   
+  
 protected:
-  void ComputeMesh();
+  void ComputeMeshFromDataPoints();
+  void ComputeDataPointsFromMesh();
   unsigned int mNbOfPoints;
   std::string mType;
   vtkSmartPointer<vtkPoints> mData;
@@ -61,6 +67,8 @@ protected:
   bool mMeshIsUpToDate;
   ///Z location of the contour
   double mZ;
+  
+  gdcm::Item * mItem;
 
 private:
   DicomRT_Contour();
index 73d25b99cf038934b5479b61edca6237e74e9075..75fc8e97fdbc517bdc61da7dc5b60bd672d9ffe9 100644 (file)
@@ -20,6 +20,8 @@
 #include "clitkDicomRT_ROI.h"
 #include <vtkSmartPointer.h>
 #include <vtkAppendPolyData.h>
+#include <vtkImageClip.h>
+#include <vtkMarchingSquares.h>
 
 #if GDCM_MAJOR_VERSION == 2
 #include "gdcmAttribute.h"
@@ -31,6 +33,7 @@ clitk::DicomRT_ROI::DicomRT_ROI()
 {
   mName = "NoName";
   mNumber = -1;
+  mImage = NULL;
   mColor.resize(3);
   mColor[0] = mColor[1] = mColor[2] = 0;
   mMeshIsUpToDate = false;
@@ -185,7 +188,8 @@ bool clitk::DicomRT_ROI::Read(gdcm::Item * itemInfo, gdcm::Item * itemContour)
       return false;
     }
   const gdcm::DataElement& csq = nestedds.GetDataElement( tcsq );
-  gdcm::SmartPointer<gdcm::SequenceOfItems> sqi2 = csq.GetValueAsSQ();
+  mContoursSequenceOfItems = csq.GetValueAsSQ();
+  gdcm::SmartPointer<gdcm::SequenceOfItems> & sqi2 = mContoursSequenceOfItems;
   if( !sqi2 || !sqi2->GetNumberOfItems() )
     {
     }
@@ -193,9 +197,9 @@ bool clitk::DicomRT_ROI::Read(gdcm::Item * itemInfo, gdcm::Item * itemContour)
 
   for(unsigned int i = 0; i < nitems; ++i)
     {
-      const gdcm::Item & j = sqi2->GetItem(i+1); // Item start at #1
+      gdcm::Item & j = sqi2->GetItem(i+1); // Item start at #1
       DicomRT_Contour::Pointer c = DicomRT_Contour::New();
-      bool b = c->Read(j);
+      bool b = c->Read(&j);
       if (b) {
         mListOfContours.push_back(c);
       }
@@ -249,7 +253,7 @@ void clitk::DicomRT_ROI::SetImage(vvImage * image)
 vtkPolyData * clitk::DicomRT_ROI::GetMesh()
 {
   if (!mMeshIsUpToDate) {
-    ComputeMesh();
+    ComputeMeshFromContour();
   }
   return mMesh;
 }
@@ -265,7 +269,7 @@ clitk::DicomRT_Contour * clitk::DicomRT_ROI::GetContour(int n)
 
 
 //--------------------------------------------------------------------
-void clitk::DicomRT_ROI::ComputeMesh()
+void clitk::DicomRT_ROI::ComputeMeshFromContour()
 {
   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
   for(unsigned int i=0; i<mListOfContours.size(); i++) {
@@ -298,12 +302,37 @@ void clitk::DicomRT_ROI::UpdateDicomItem()
   gdcm::DataSet & ds = mItemInfo->GetNestedDataSet();  
   ds.Replace(de);
 
+  // From MESH to CONTOURS
+  ComputeContoursFromImage();
+
   // Update contours
+  DD(mListOfContours.size());
   for(uint i=0; i<mListOfContours.size(); i++) {
     DD(i);
-    DicomRT_Contour::Pointer roi = mListOfContours[i];
-    roi->UpdateDicomItem(mItemContour);
+    DicomRT_Contour::Pointer contour = mListOfContours[i];
+    contour->UpdateDicomItem();//mItemContour);
   }
+
+  // Nb of contours
+  unsigned int nitems = mContoursSequenceOfItems->GetNumberOfItems();
+  DD(nitems);
+
+  // Write [Contour Sequence] = 0x3006,0x0040)
+  gdcm::DataSet & dsc = mItemContour->GetNestedDataSet();
+  gdcm::Tag tcsq(0x3006,0x0040);
+  const gdcm::DataElement& csq = dsc.GetDataElement( tcsq );
+  gdcm::DataElement dec(csq);
+  dec.SetValue(*mContoursSequenceOfItems);
+  dsc.Replace(dec);
+
+  gdcm::DataSet & a = mContoursSequenceOfItems->GetItem(1).GetNestedDataSet();
+  gdcm::Attribute<0x3006,0x0050> at;
+  gdcm::Tag tcontourdata(0x3006,0x0050);
+  gdcm::DataElement contourdata = a.GetDataElement( tcontourdata );
+  at.SetFromDataElement( contourdata );
+  const double* points = at.GetValues();
+  DD(points[0]);
+
 }
 //--------------------------------------------------------------------
 #endif
@@ -340,3 +369,49 @@ vvImage * clitk::DicomRT_ROI::GetImage() const
   return mImage;
 }
 //--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+void clitk::DicomRT_ROI::ComputeContoursFromImage()
+{
+  DD("ComputeMeshFromImage");
+
+  // Check that an image is loaded
+  if (!mImage) return;
+
+  // Only consider 3D here
+  if (mImage->GetNumberOfDimensions() != 3) {
+    FATAL("DicomRT_ROI::ComputeMeshFromImage only work with 3D images");
+  }
+
+  // Get the VTK image
+  vtkImageData * image = mImage->GetVTKImages()[0];
+  
+  // Get initial extend for the clipping
+  vtkSmartPointer<vtkImageClip> clipper = vtkSmartPointer<vtkImageClip>::New();
+  clipper->SetInput(image);
+  int* extent = image->GetExtent();
+  DDV(extent, 6);
+  //  std::vector<int> extend;
+
+  // Prepare the marching squares
+  vtkSmartPointer<vtkMarchingSquares> squares = vtkSmartPointer<vtkMarchingSquares>::New();
+  squares->SetInput(clipper->GetOutput());
+
+  // Loop on slice
+  uint n = image->GetDimensions()[2];
+  DD(n);
+  for(uint i=0; i<n; i++) {
+    // Clip to the current slice
+    extent[4] = extent[5] = image->GetOrigin()[2]+i*image->GetSpacing()[2];
+    DDV(extent, 6);
+    clipper->SetOutputWholeExtent(extent[0],extent[1],extent[2],
+                                extent[3],extent[4],extent[5]);
+    
+
+    squares->Update();
+    DD(squares->GetNumberOfContours());
+    mListOfContours[i]->SetMesh(squares->GetOutput());
+  }
+}
+//--------------------------------------------------------------------
index 811735218d3327e982d873d71db1d301cc942dc0..f0be23d74e47777647dc6714cf512415ec2193a9 100644 (file)
@@ -60,7 +60,8 @@ public:
   DicomRT_Contour* GetContour(int n);
 
   // Compute a vtk mesh from the dicom contours
-  void ComputeMesh();
+  void ComputeMeshFromContour();
+  void ComputeContoursFromImage();
   
   // Indicate if the mesh is uptodate according to the dicom
   void SetDicomUptodateFlag(bool b) { m_DicomUptodateFlag = b; }
@@ -91,6 +92,7 @@ protected:
 #if GDCM_MAJOR_VERSION == 2
   gdcm::Item * mItemInfo;
   gdcm::Item * mItemContour;
+  gdcm::SmartPointer<gdcm::SequenceOfItems> mContoursSequenceOfItems;
 #endif
 
 private:
index b10984ad56df49bf0c07786f5a1466ddcc749974..4ac4705ef73a673b0487980e569aac56ba3f01ef 100644 (file)
@@ -180,8 +180,23 @@ void clitk::DicomRT_StructureSet::Write(const std::string & filename)
   const gdcm::DataElement &roicsq = ds.GetDataElement( troicsq );
   gdcm::DataElement de2(roicsq);
   de2.SetValue(*mROIContoursSequenceOfItems);
-  ds.Replace(de);
-  
+  ds.Replace(de2);
+
+  //DEBUG
+  gdcm::DataSet & a = mROIContoursSequenceOfItems->GetItem(1).GetNestedDataSet();
+  gdcm::Tag tcsq(0x3006,0x0040);
+  const gdcm::DataElement& csq = a.GetDataElement( tcsq );
+  gdcm::SmartPointer<gdcm::SequenceOfItems> sqi2 = csq.GetValueAsSQ();
+  gdcm::Item & j = sqi2->GetItem(1);
+  gdcm::DataSet & b = j.GetNestedDataSet();
+  gdcm::Attribute<0x3006,0x0050> at;
+  gdcm::Tag tcontourdata(0x3006,0x0050);
+  gdcm::DataElement contourdata = b.GetDataElement( tcontourdata );
+  at.SetFromDataElement( contourdata );
+  const double* points = at.GetValues();
+  DD(points[0]);
+
+
   // Write dicom
   gdcm::Writer writer;
   //writer.CheckFileMetaInformationOff();
index b4e921269a1ff0dc7d874176a9c5b967ac711502..6ab72e2ed6916e48a6ce8978d7016508baea85e8 100644 (file)
@@ -24,6 +24,7 @@
 // clitk
 #include "clitkImage2DicomRTStructFilter.h"
 #include "clitkImageCommon.h"
+#include "vvFromITK.h"
 
 // vtk
 #include <vtkPolyDataToImageStencil.h>
@@ -32,6 +33,8 @@
 #include <vtkLinearExtrusionFilter.h>
 #include <vtkMetaImageWriter.h>
 #include <vtkImageData.h>
+
+// itk
 #include <itkImage.h>
 #include <itkVTKImageToImageFilter.h>
 
@@ -57,9 +60,44 @@ clitk::Image2DicomRTStructFilter<PixelType>::~Image2DicomRTStructFilter()
 template<class PixelType>
 void clitk::Image2DicomRTStructFilter<PixelType>::Update() 
 {
-  DD("Image2DicomRTStructFilter::Update");
-
-
+  DD("Image2DicomRTStructFilter::GenerateData");
+
+  // Read DicomRTStruct
+  std::string filename = "RS.zzQAnotmt_french01_.dcm";
+  clitk::DicomRT_StructureSet::Pointer structset = clitk::DicomRT_StructureSet::New();
+  structset->Read(filename);
+  
+  DD(structset->GetName());
+  clitk::DicomRT_ROI * roi = structset->GetROIFromROINumber(1); // Aorta
+  DD(roi->GetName());
+  DD(roi->GetROINumber());
+
+  // Add an image to the roi
+  vvImage::Pointer im = vvImageFromITK<3, PixelType>(m_Input);
+  roi->SetImage(im);
+
+  // Get one contour
+  DD("Compute Mesh");
+  roi->ComputeMeshFromImage();
+  vtkSmartPointer<vtkPolyData> mesh = roi->GetMesh();
+  DD("done");
+  
+  // Change the mesh (shift by 10);
+  // const vtkSmartPointer<vtkPoints> & points = mesh->GetPoints();
+  // for(uint i=0; i<mesh->GetNumberOfVerts (); i++) {
+  //   DD(i);
+  //   double * p = points->GetPoint(i);
+  //   p[0] += 30;
+  //   points->SetPoint(i, p);
+  // }
+  roi->SetName("TOTO");
+  roi->SetDicomUptodateFlag(false); // indicate that dicom info must be updated from the mesh.
+    
+  // Convert to dicom ?
+  DD("TODO");
+  
+  // Write
+  structset->Write("toto.dcm");
 }
 //--------------------------------------------------------------------
 
index e5529bf4c1852a86e4975e0a33ec2adffdedd94c..4b1d3c2906b13dac04cd915b881f66673f2f446b 100644 (file)
@@ -18,7 +18,10 @@ clitk::ExtractLymphStationsFilter<TImageType>::
 ExtractStation_3A()
 {
   if (!CheckForStation("3A")) return;
-    
+
+  StartNewStep("Station 3A");
+  StartSubStep();   
+
   // Get the current support 
   StartNewStep("[Station 3A] Get the current 3A suppport");
   m_Working_Support = m_ListOfSupports["S3A"];
@@ -28,7 +31,6 @@ ExtractStation_3A()
   ExtractStation_3A_AntPost_S5();
   ExtractStation_3A_AntPost_S6();
   ExtractStation_3A_AntPost_Superiorly();
-
   ExtractStation_3A_Remove_Structures();
 
   Remove_Structures("3A", "Aorta");
@@ -38,8 +40,10 @@ ExtractStation_3A()
   Remove_Structures("3A", "CommonCarotidArteryLeft");
   Remove_Structures("3A", "CommonCarotidArteryRight");
   Remove_Structures("3A", "BrachioCephalicArtery");
-  //  Remove_Structures("3A", "BrachioCephalicVein"); ?
 
+  ExtractStation_3A_PostToBones();
+  
+  
 
   // ExtractStation_3A_Ant_Limits(); --> No, already in support; to remove
   // ExtractStation_3A_Post_Limits(); --> No, more complex, see Vessels etc
@@ -48,6 +52,7 @@ ExtractStation_3A()
   writeImage<MaskImageType>(m_ListOfStations["3A"], "seg/Station3A.mhd");
   GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd"); 
   WriteAFDB(); 
+  StopSubStep();
 }
 //--------------------------------------------------------------------
 
@@ -112,24 +117,6 @@ ExtractStation_3A_AntPost_S5()
   MaskImagePointer SVC = GetAFDB()->template GetImage <MaskImageType>("SVC");
 
   // Trial in 3D -> difficulties superiorly. Stay slice by slice.
-  /*
-  typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->VerboseStepFlagOff();
-  relPosFilter->WriteStepFlagOff();
-  relPosFilter->SetBackgroundValue(GetBackgroundValue());
-  relPosFilter->SetInput(m_Working_Support); 
-  relPosFilter->SetInputObject(SVC); 
-  relPosFilter->AddOrientationTypeString("PostTo");
-  relPosFilter->SetInverseOrientationFlag(true);
-  relPosFilter->SetIntermediateSpacing(4);
-  relPosFilter->SetIntermediateSpacingFlag(false);
-  relPosFilter->SetFuzzyThreshold(0.5);
-  //  relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
-  relPosFilter->Update();
-  m_Working_Support = relPosFilter->GetOutput();
-  */
-  
   // Slice by slice not post to SVC. Use initial spacing
   m_Working_Support = 
     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, SVC, 2, 
@@ -308,6 +295,8 @@ ExtractStation_3A_Remove_Structures()
   // resize like support, extract slices
   // while single CCL -> remove
   // when two remove only the most post
+  MaskImagePointer BrachioCephalicVein = 
+    GetAFDB()->template GetImage <MaskImageType>("BrachioCephalicVein");
   BrachioCephalicVein = clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, 
                                                               m_Working_Support, 
                                                               GetBackgroundValue());
@@ -318,27 +307,31 @@ ExtractStation_3A_Remove_Structures()
   for(uint i=0; i<slices.size(); i++) {
     // Labelize slices_BCV
     slices_BCV[i] = Labelize<MaskSliceType>(slices_BCV[i], 0, true, 1);
+
     // Compute centroids
     std::vector<typename MaskSliceType::PointType> centroids;
     ComputeCentroids<MaskSliceType>(slices_BCV[i], GetBackgroundValue(), centroids);
-    if (centroids.size() > 1) {
+
+    // If several centroid, select the one most anterior
+    if (centroids.size() > 2) {
       // Only keep the one most post
       typename MaskSliceType::PixelType label;
       if (centroids[1][1] > centroids[2][1]) {
-        label = 1;
-      }
-      else {
         label = 2;
       }
-      
-      HERE
-
-      slices_BCV[i] = clitk::SetBackground<MaskSliceType>(slices_BCV[i], slices_BCV[i], 
-                                                          label, 
-                                                          GetBackgroundValue(), true);
+      else {
+        label = 1;
+      }      
+      // "remove" the CCL 
+      slices_BCV[i] = clitk::SetBackground<MaskSliceType, MaskSliceType>(slices_BCV[i], 
+                                                                         slices_BCV[i], 
+                                                                         label, 
+                                                                         GetBackgroundValue(), 
+                                                                         true);
     }
+    
     // Remove from the support
-    slices[i] = clitk::AndNot(slices[i], slices_BCV[i], GetBackgroundValue());
+    clitk::AndNot<MaskSliceType>(slices[i], slices_BCV[i], GetBackgroundValue());
   }
   
   // Joint
index 58549dfcfce17ac5919fd4d51335b64d01ccfa0d..3478b490959b166d75fa33b169cbdacc0026fbdc 100644 (file)
@@ -17,6 +17,9 @@ ExtractStation_3P()
 {
   if (!CheckForStation("3P")) return;
 
+  StartNewStep("Station 3P");
+  StartSubStep();  
+
   // Get the current support 
   StartNewStep("[Station 3P] Get the current 3P suppport");
   m_Working_Support = m_ListOfSupports["S3P"];
@@ -38,6 +41,7 @@ ExtractStation_3P()
   writeImage<MaskImageType>(m_ListOfStations["3P"], "seg/Station3P.mhd");
   GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd"); 
   WriteAFDB();   
+  StopSubStep();
 }
 //--------------------------------------------------------------------
 
index 7402a626bebd7e86bd06c92c43f843b5291f930e..7d5a15c7f97928597979364049062296780980b3 100644 (file)
@@ -24,19 +24,22 @@ void
 clitk::ExtractLymphStationsFilter<TImageType>::
 ExtractStation_8() 
 {
-  if (CheckForStation("8")) {
-    ExtractStation_8_SI_Limits();         // OK, validated
-    ExtractStation_8_Ant_Limits();        // OK, validated
-    ExtractStation_8_Left_Sup_Limits();   // OK, validated
-    ExtractStation_8_Left_Inf_Limits();   // OK, validated
-    ExtractStation_8_Single_CCL_Limits(); // OK, validated
-    ExtractStation_8_Remove_Structures(); // OK, validated
-
-    // Store image filenames into AFDB 
-    writeImage<MaskImageType>(m_ListOfStations["8"], "seg/Station8.mhd");
-    GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd");  
-    WriteAFDB();
-  }
+  if (!CheckForStation("8")) return;
+
+  StartNewStep("Station 8");
+  StartSubStep();   
+  ExtractStation_8_SI_Limits();         // OK, validated
+  ExtractStation_8_Ant_Limits();        // OK, validated
+  ExtractStation_8_Left_Sup_Limits();   // OK, validated
+  ExtractStation_8_Left_Inf_Limits();   // OK, validated
+  ExtractStation_8_Single_CCL_Limits(); // OK, validated
+  ExtractStation_8_Remove_Structures(); // OK, validated
+  
+  // Store image filenames into AFDB 
+  writeImage<MaskImageType>(m_ListOfStations["8"], "seg/Station8.mhd");
+  GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd");  
+  WriteAFDB();
+  StopSubStep();
 }
 //--------------------------------------------------------------------
 
index 8bbc3013f424cda711e3b500c4d4a28678e66507..7c2376aa8bfe56471aac85f90fc9bf981ba5c965 100644 (file)
@@ -200,6 +200,8 @@ namespace clitk {
     void ExtractStation_3A_AntPost_S5();
     void ExtractStation_3A_AntPost_S6();
     void ExtractStation_3A_AntPost_Superiorly();
+    void ExtractStation_3A_Remove_Structures();
+    
     void ExtractStation_3A_SI_Limits();
     void ExtractStation_3A_Ant_Limits();
     void ExtractStation_3A_Post_Limits();
index 904bfc62011b684a9835d1c00f65b73553703df9..520038c89785c58abc9eab9b8a35c8bb1c4bbaae 100644 (file)
@@ -118,23 +118,13 @@ GenerateOutputInformation() {
   }
 
   // Extract Station8
-  StartNewStep("Station 8");
-  StartSubStep(); 
   ExtractStation_8();
-  StopSubStep();
 
   // Extract Station3P
-  StartNewStep("Station 3P");
-  StartSubStep(); 
   ExtractStation_3P();
-  StopSubStep();
 
   // Extract Station3A
-  StartNewStep("Station 3A");
-  StartSubStep(); 
   ExtractStation_3A();
-  StopSubStep();
-
 
   // HERE
 
index e48272c779ca68d47b06bb6380cedb8bb07c151e..d71491ba8d12ea5e6ef26a6f45316a4b8fb44a73 100644 (file)
@@ -39,8 +39,9 @@ int main(int argc, char * argv[]) {
   
   // Write result
   clitk::DicomRT_StructureSet::Pointer s = filter.GetDicomRTStruct();
-  s->Write(args_info.output_arg);
+  //  s->Write(args_info.output_arg);
 
   // This is the end my friend 
   return 0;
 }
+//--------------------------------------------------------------------