From 8f785a7646e65eb73b28e274871b2aafcb34d85c Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 28 Jul 2011 12:58:40 +0200 Subject: [PATCH] first version of Image2DicomRT (dos not work) --- common/clitkDicomRT_Contour.cxx | 4 + common/clitkDicomRT_Contour.h | 12 ++- common/clitkDicomRT_ROI.cxx | 89 +++++++++++++++++-- common/clitkDicomRT_ROI.h | 4 +- common/clitkDicomRT_StructureSet.cxx | 19 +++- common/clitkImage2DicomRTStructFilter.txx | 44 ++++++++- segmentation/clitkExtractLymphStation_3A.txx | 57 ++++++------ segmentation/clitkExtractLymphStation_3P.txx | 4 + segmentation/clitkExtractLymphStation_8.txx | 29 +++--- .../clitkExtractLymphStationsFilter.h | 2 + .../clitkExtractLymphStationsFilter.txx | 10 --- tools/clitkImage2DicomRTStruct.cxx | 3 +- 12 files changed, 206 insertions(+), 71 deletions(-) diff --git a/common/clitkDicomRT_Contour.cxx b/common/clitkDicomRT_Contour.cxx index a4a69dd..8617170 100644 --- a/common/clitkDicomRT_Contour.cxx +++ b/common/clitkDicomRT_Contour.cxx @@ -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 diff --git a/common/clitkDicomRT_Contour.h b/common/clitkDicomRT_Contour.h index 0484da6..91d75d6 100644 --- a/common/clitkDicomRT_Contour.h +++ b/common/clitkDicomRT_Contour.h @@ -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 mData; @@ -61,6 +67,8 @@ protected: bool mMeshIsUpToDate; ///Z location of the contour double mZ; + + gdcm::Item * mItem; private: DicomRT_Contour(); diff --git a/common/clitkDicomRT_ROI.cxx b/common/clitkDicomRT_ROI.cxx index 73d25b9..75fc8e9 100644 --- a/common/clitkDicomRT_ROI.cxx +++ b/common/clitkDicomRT_ROI.cxx @@ -20,6 +20,8 @@ #include "clitkDicomRT_ROI.h" #include #include +#include +#include #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 sqi2 = csq.GetValueAsSQ(); + mContoursSequenceOfItems = csq.GetValueAsSQ(); + gdcm::SmartPointer & 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 append = vtkSmartPointer::New(); for(unsigned int i=0; iGetNestedDataSet(); ds.Replace(de); + // From MESH to CONTOURS + ComputeContoursFromImage(); + // Update contours + DD(mListOfContours.size()); for(uint i=0; iUpdateDicomItem(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 clipper = vtkSmartPointer::New(); + clipper->SetInput(image); + int* extent = image->GetExtent(); + DDV(extent, 6); + // std::vector extend; + + // Prepare the marching squares + vtkSmartPointer squares = vtkSmartPointer::New(); + squares->SetInput(clipper->GetOutput()); + + // Loop on slice + uint n = image->GetDimensions()[2]; + DD(n); + for(uint i=0; iGetOrigin()[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()); + } +} +//-------------------------------------------------------------------- diff --git a/common/clitkDicomRT_ROI.h b/common/clitkDicomRT_ROI.h index 8117352..f0be23d 100644 --- a/common/clitkDicomRT_ROI.h +++ b/common/clitkDicomRT_ROI.h @@ -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 mContoursSequenceOfItems; #endif private: diff --git a/common/clitkDicomRT_StructureSet.cxx b/common/clitkDicomRT_StructureSet.cxx index b10984a..4ac4705 100644 --- a/common/clitkDicomRT_StructureSet.cxx +++ b/common/clitkDicomRT_StructureSet.cxx @@ -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 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(); diff --git a/common/clitkImage2DicomRTStructFilter.txx b/common/clitkImage2DicomRTStructFilter.txx index b4e9212..6ab72e2 100644 --- a/common/clitkImage2DicomRTStructFilter.txx +++ b/common/clitkImage2DicomRTStructFilter.txx @@ -24,6 +24,7 @@ // clitk #include "clitkImage2DicomRTStructFilter.h" #include "clitkImageCommon.h" +#include "vvFromITK.h" // vtk #include @@ -32,6 +33,8 @@ #include #include #include + +// itk #include #include @@ -57,9 +60,44 @@ clitk::Image2DicomRTStructFilter::~Image2DicomRTStructFilter() template void clitk::Image2DicomRTStructFilter::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 mesh = roi->GetMesh(); + DD("done"); + + // Change the mesh (shift by 10); + // const vtkSmartPointer & points = mesh->GetPoints(); + // for(uint i=0; iGetNumberOfVerts (); 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"); } //-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStation_3A.txx b/segmentation/clitkExtractLymphStation_3A.txx index e5529bf..4b1d3c2 100644 --- a/segmentation/clitkExtractLymphStation_3A.txx +++ b/segmentation/clitkExtractLymphStation_3A.txx @@ -18,7 +18,10 @@ clitk::ExtractLymphStationsFilter:: 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(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 ("SVC"); // Trial in 3D -> difficulties superiorly. Stay slice by slice. - /* - typedef clitk::AddRelativePositionConstraintToLabelImageFilter 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(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 ("BrachioCephalicVein"); BrachioCephalicVein = clitk::ResizeImageLike(BrachioCephalicVein, m_Working_Support, GetBackgroundValue()); @@ -318,27 +307,31 @@ ExtractStation_3A_Remove_Structures() for(uint i=0; i(slices_BCV[i], 0, true, 1); + // Compute centroids std::vector centroids; ComputeCentroids(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(slices_BCV[i], slices_BCV[i], - label, - GetBackgroundValue(), true); + else { + label = 1; + } + // "remove" the CCL + slices_BCV[i] = clitk::SetBackground(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(slices[i], slices_BCV[i], GetBackgroundValue()); } // Joint diff --git a/segmentation/clitkExtractLymphStation_3P.txx b/segmentation/clitkExtractLymphStation_3P.txx index 58549df..3478b49 100644 --- a/segmentation/clitkExtractLymphStation_3P.txx +++ b/segmentation/clitkExtractLymphStation_3P.txx @@ -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(m_ListOfStations["3P"], "seg/Station3P.mhd"); GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd"); WriteAFDB(); + StopSubStep(); } //-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStation_8.txx b/segmentation/clitkExtractLymphStation_8.txx index 7402a62..7d5a15c 100644 --- a/segmentation/clitkExtractLymphStation_8.txx +++ b/segmentation/clitkExtractLymphStation_8.txx @@ -24,19 +24,22 @@ void clitk::ExtractLymphStationsFilter:: 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(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(m_ListOfStations["8"], "seg/Station8.mhd"); + GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd"); + WriteAFDB(); + StopSubStep(); } //-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStationsFilter.h b/segmentation/clitkExtractLymphStationsFilter.h index 8bbc301..7c2376a 100644 --- a/segmentation/clitkExtractLymphStationsFilter.h +++ b/segmentation/clitkExtractLymphStationsFilter.h @@ -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(); diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx index 904bfc6..520038c 100644 --- a/segmentation/clitkExtractLymphStationsFilter.txx +++ b/segmentation/clitkExtractLymphStationsFilter.txx @@ -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 diff --git a/tools/clitkImage2DicomRTStruct.cxx b/tools/clitkImage2DicomRTStruct.cxx index e48272c..d71491b 100644 --- a/tools/clitkImage2DicomRTStruct.cxx +++ b/tools/clitkImage2DicomRTStruct.cxx @@ -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; } +//-------------------------------------------------------------------- -- 2.47.1