common/clitkDicomRT_ROI.cxx

   1 /*=========================================================================
   2   Program:         vv http://www.creatis.insa-lyon.fr/rio/vv
   3   Main authors :   XX XX XX
   4
   5   Authors belongs to:
   6   - University of LYON           http://www.universite-lyon.fr/
   7   - Léon Bérard cancer center    http://www.centreleonberard.fr
   8   - CREATIS CNRS laboratory      http://www.creatis.insa-lyon.fr
   9
  10   This software is distributed WITHOUT ANY WARRANTY; without even
  11   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12   PURPOSE.  See the copyright notices for more information.
  13
  14   It is distributed under dual licence
  15   - BSD       http://www.opensource.org/licenses/bsd-license.php
  16   - CeCILL-B  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17
  18   =========================================================================*/
  19
  20 #include "clitkDicomRT_ROI.h"
  21 #include <vtkSmartPointer.h>
  22 #include <vtkAppendPolyData.h>
  23 #include <vtkImageClip.h>
  24 #include <vtkMarchingSquares.h>
  25 #include <vtkPolyDataWriter.h>
  26
  27 #if GDCM_MAJOR_VERSION == 2
  28 #include "gdcmAttribute.h"
  29 #include "gdcmItem.h"
  30 #endif
  31
  32 //--------------------------------------------------------------------
  33 clitk::DicomRT_ROI::DicomRT_ROI()
  34 {
  35   mName = "NoName";
  36   mNumber = -1;
  37   mImage = NULL;
  38   mColor.resize(3);
  39   mColor[0] = mColor[1] = mColor[2] = 0;
  40   mMeshIsUpToDate = false;
  41   mBackgroundValue = 0;
  42   mForegroundValue = 1;
  43   SetDicomUptodateFlag(false);
  44 }
  45 //--------------------------------------------------------------------
  46
  47
  48 //--------------------------------------------------------------------
  49 clitk::DicomRT_ROI::~DicomRT_ROI()
  50 {
  51 }
  52 //--------------------------------------------------------------------
  53
  54
  55 //--------------------------------------------------------------------
  56 void clitk::DicomRT_ROI::SetDisplayColor(double r, double v, double b)
  57 {
  58   mColor.resize(3);
  59   mColor[0] = r;
  60   mColor[1] = v;
  61   mColor[2] = b;
  62 }
  63 //--------------------------------------------------------------------
  64
  65
  66 //--------------------------------------------------------------------
  67 int clitk::DicomRT_ROI::GetROINumber() const
  68 {
  69   return mNumber;
  70 }
  71 //--------------------------------------------------------------------
  72
  73
  74 //--------------------------------------------------------------------
  75 const std::string & clitk::DicomRT_ROI::GetName() const
  76 {
  77   return mName;
  78 }
  79 //--------------------------------------------------------------------
  80
  81
  82 //--------------------------------------------------------------------
  83 const std::string & clitk::DicomRT_ROI::GetFilename() const
  84 {
  85   return mFilename;
  86 }
  87 //--------------------------------------------------------------------
  88
  89
  90 //--------------------------------------------------------------------
  91 const std::vector<double> & clitk::DicomRT_ROI::GetDisplayColor() const
  92 {
  93   return mColor;
  94 }
  95 //--------------------------------------------------------------------
  96
  97
  98 //--------------------------------------------------------------------
  99 void clitk::DicomRT_ROI::Print(std::ostream & os) const
 100 {
 101   os << "ROI " << mNumber << "\t" << mName
 102      << "\t(" << mColor[0] << " " << mColor[1] << " " << mColor[2] << ")"
 103      << "\t Contours = " << mListOfContours.size() << std::endl;
 104 }
 105 //--------------------------------------------------------------------
 106
 107
 108 //--------------------------------------------------------------------
 109 void clitk::DicomRT_ROI::SetBackgroundValueLabelImage(double bg)
 110 {
 111   mBackgroundValue = bg;
 112 }
 113 //--------------------------------------------------------------------
 114
 115
 116 //--------------------------------------------------------------------
 117 double clitk::DicomRT_ROI::GetBackgroundValueLabelImage() const
 118 {
 119   return mBackgroundValue;
 120 }
 121 //--------------------------------------------------------------------
 122
 123
 124 //--------------------------------------------------------------------
 125 void clitk::DicomRT_ROI::SetForegroundValueLabelImage(double bg)
 126 {
 127   mForegroundValue = bg;
 128 }
 129 //--------------------------------------------------------------------
 130
 131
 132 //--------------------------------------------------------------------
 133 double clitk::DicomRT_ROI::GetForegroundValueLabelImage() const
 134 {
 135   return mForegroundValue;
 136 }
 137 //--------------------------------------------------------------------
 138
 139
 140 //--------------------------------------------------------------------
 141 #if GDCM_MAJOR_VERSION == 2
 142 bool clitk::DicomRT_ROI::Read(gdcm::Item * itemInfo, gdcm::Item * itemContour)
 143 {
 144   // Keep dicom item
 145   mItemInfo = itemInfo;
 146   mItemContour = itemContour;
 147   // DD(mItemInfo);
 148
 149   // ROI number [Referenced ROI Number]
 150   const gdcm::DataSet & nesteddsInfo = mItemInfo->GetNestedDataSet();
 151   gdcm::Attribute<0x3006,0x0022> roinumber;
 152   roinumber.SetFromDataSet( nesteddsInfo );
 153   int nb1 = roinumber.GetValue();
 154
 155   // Check this is the same with the other item
 156   const gdcm::DataSet & nestedds = mItemContour->GetNestedDataSet();
 157   gdcm::Attribute<0x3006,0x0084> referencedroinumber;
 158   referencedroinumber.SetFromDataSet( nestedds );
 159   int nb2 = referencedroinumber.GetValue();
 160
 161   // Must never be different
 162   if (nb1 != nb2) {
 163     DD(nb2);
 164     DD(nb1);
 165     FATAL("nb1 must equal nb2" << std::endl);
 166   }
 167   mNumber = nb1;
 168
 169   // Retrieve ROI Name (in the info item)
 170   gdcm::Attribute<0x3006,0x26> roiname;
 171   roiname.SetFromDataSet( nesteddsInfo );
 172   mName = roiname.GetValue();
 173   // DD(mName);
 174
 175   // ROI Color [ROI Display Color]
 176   gdcm::Attribute<0x3006,0x002a> color = {};
 177   color.SetFromDataSet( nestedds );
 178   assert( color.GetNumberOfValues() == 3 );
 179   mColor[0] = color.GetValue(0);
 180   mColor[1] = color.GetValue(1);
 181   mColor[2] = color.GetValue(2);
 182
 183   // Read contours [Contour Sequence]
 184   gdcm::Tag tcsq(0x3006,0x0040);
 185   if( !nestedds.FindDataElement( tcsq ) )
 186     {
 187       std::cerr << "Warning. Could not read contour for structure <" << mName << ">, number" << mNumber << " ? I ignore it" << std::endl;
 188       SetDicomUptodateFlag(true);
 189       return false;
 190     }
 191   const gdcm::DataElement& csq = nestedds.GetDataElement( tcsq );
 192   mContoursSequenceOfItems = csq.GetValueAsSQ();
 193   gdcm::SmartPointer<gdcm::SequenceOfItems> & sqi2 = mContoursSequenceOfItems;
 194   if( !sqi2 || !sqi2->GetNumberOfItems() )
 195     {
 196     }
 197   unsigned int nitems = sqi2->GetNumberOfItems();
 198
 199   for(unsigned int i = 0; i < nitems; ++i)
 200     {
 201       gdcm::Item & j = sqi2->GetItem(i+1); // Item start at #1
 202       DicomRT_Contour::Pointer c = DicomRT_Contour::New();
 203       bool b = c->Read(&j);
 204       if (b) {
 205         mListOfContours.push_back(c);
 206       }
 207     }
 208   SetDicomUptodateFlag(true);
 209   return true;
 210 }
 211 #else
 212 void clitk::DicomRT_ROI::Read(std::map<int, std::string> & rois, gdcm::SQItem * item)
 213 {
 214   // ROI number [Referenced ROI Number]
 215   mNumber = atoi(item->GetEntryValue(0x3006,0x0084).c_str());
 216
 217   // Retrieve ROI Name
 218   mName = rois[mNumber];
 219
 220   // ROI Color [ROI Display Color]
 221   mColor = clitk::parse_string<double>(item->GetEntryValue(0x3006,0x002a),'\\');
 222
 223   // Read contours [Contour Sequence]
 224   gdcm::SeqEntry * contours=item->GetSeqEntry(0x3006,0x0040);
 225   if (contours) {
 226     int i=0;
 227     for(gdcm::SQItem* j=contours->GetFirstSQItem(); j!=0; j=contours->GetNextSQItem()) {
 228       DicomRT_Contour::Pointer c = DicomRT_Contour::New();
 229       bool b = c->Read(j);
 230       if (b) {
 231         mListOfContours.push_back(c);
 232       }
 233       ++i;
 234     }
 235   }
 236   else {
 237     std::cerr << "Warning. Could not read contour for structure <" << mName << ">, number" << mNumber << " ? I ignore it" << std::endl;
 238   }
 239   SetDicomUptodateFlag(true);
 240 }
 241 #endif
 242 //--------------------------------------------------------------------
 243
 244
 245 //--------------------------------------------------------------------
 246 void clitk::DicomRT_ROI::SetImage(vvImage::Pointer image)
 247 {
 248   mImage = image;
 249 }
 250 //--------------------------------------------------------------------
 251
 252
 253 //--------------------------------------------------------------------
 254 vtkPolyData * clitk::DicomRT_ROI::GetMesh()
 255 {
 256   if (!mMeshIsUpToDate) {
 257     ComputeMeshFromContour();
 258   }
 259   return mMesh;
 260 }
 261 //--------------------------------------------------------------------
 262
 263
 264 //--------------------------------------------------------------------
 265 clitk::DicomRT_Contour * clitk::DicomRT_ROI::GetContour(int n)
 266 {
 267   return mListOfContours[n];
 268 }
 269 //--------------------------------------------------------------------
 270
 271
 272 //--------------------------------------------------------------------
 273 void clitk::DicomRT_ROI::ComputeMeshFromContour()
 274 {
 275   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
 276   for(unsigned int i=0; i<mListOfContours.size(); i++) {
 277     append->AddInput(mListOfContours[i]->GetMesh());
 278   }
 279   append->Update();
 280
 281   mMesh = vtkSmartPointer<vtkPolyData>::New();
 282   mMesh->DeepCopy(append->GetOutput());
 283   mMeshIsUpToDate = true;
 284 }
 285 //--------------------------------------------------------------------
 286
 287
 288 #if GDCM_MAJOR_VERSION == 2
 289
 290 //--------------------------------------------------------------------
 291 void clitk::DicomRT_ROI::UpdateDicomItem()
 292 {
 293   if (GetDicomUptoDateFlag()) return;
 294   DD("ROI::UpdateDicomItem");
 295   DD(GetName());
 296
 297   // From now, only some item can be modified
 298
 299   // Set ROI Name 0x3006,0x26>
 300   gdcm::Attribute<0x3006,0x26> roiname;
 301   roiname.SetValue(GetName());
 302   gdcm::DataElement de = roiname.GetAsDataElement();
 303   gdcm::DataSet & ds = mItemInfo->GetNestedDataSet();
 304   ds.Replace(de);
 305
 306   // From MESH to CONTOURS
 307   ComputeContoursFromImage();
 308
 309   // Update contours
 310   DD(mListOfContours.size());
 311   for(uint i=0; i<mListOfContours.size(); i++) {
 312     DD(i);
 313     DicomRT_Contour::Pointer contour = mListOfContours[i];
 314     contour->UpdateDicomItem();//mItemContour);
 315   }
 316
 317   // Nb of contours
 318   unsigned int nitems = mContoursSequenceOfItems->GetNumberOfItems();
 319   DD(nitems);
 320
 321   // Write [Contour Sequence] = 0x3006,0x0040)
 322   gdcm::DataSet & dsc = mItemContour->GetNestedDataSet();
 323   gdcm::Tag tcsq(0x3006,0x0040);
 324   const gdcm::DataElement& csq = dsc.GetDataElement( tcsq );
 325   gdcm::DataElement dec(csq);
 326   dec.SetValue(*mContoursSequenceOfItems);
 327   dsc.Replace(dec);
 328
 329   gdcm::DataSet & a = mContoursSequenceOfItems->GetItem(1).GetNestedDataSet();
 330   gdcm::Attribute<0x3006,0x0050> at;
 331   gdcm::Tag tcontourdata(0x3006,0x0050);
 332   gdcm::DataElement contourdata = a.GetDataElement( tcontourdata );
 333   at.SetFromDataElement( contourdata );
 334   const double* points = at.GetValues();
 335   DD(points[0]);
 336
 337 }
 338 //--------------------------------------------------------------------
 339 #endif
 340
 341 //--------------------------------------------------------------------
 342 void clitk::DicomRT_ROI::SetFromBinaryImage(vvImage::Pointer image, int n,
 343                                             std::string name,
 344                                             std::vector<double> color,
 345                                             std::string filename)
 346 {
 347
 348   // ROI number [Referenced ROI Number]
 349   mNumber = n;
 350
 351   // ROI Name
 352   mName = name;
 353   mFilename = filename;
 354
 355   // ROI Color [ROI Display Color]
 356   mColor = color;
 357
 358   // No contours [Contour Sequence]
 359   mListOfContours.clear();
 360
 361   // Set image
 362   mImage = image;
 363 }
 364 //--------------------------------------------------------------------
 365
 366
 367 //--------------------------------------------------------------------
 368 vvImage * clitk::DicomRT_ROI::GetImage() const
 369 {
 370   return mImage;
 371 }
 372 //--------------------------------------------------------------------
 373
 374
 375 //--------------------------------------------------------------------
 376 void clitk::DicomRT_ROI::ComputeContoursFromImage()
 377 {
 378   DD("ComputeMeshFromImage");
 379
 380   // Check that an image is loaded
 381   if (!mImage) return;
 382
 383   // Only consider 3D here
 384   if (mImage->GetNumberOfDimensions() != 3) {
 385     FATAL("DicomRT_ROI::ComputeMeshFromImage only work with 3D images");
 386   }
 387
 388   // Get the VTK image
 389   vtkImageData * image = mImage->GetVTKImages()[0];
 390
 391   // Get initial extend for the clipping
 392   vtkSmartPointer<vtkImageClip> clipper = vtkSmartPointer<vtkImageClip>::New();
 393   clipper->SetInput(image);
 394   int* extent = image->GetExtent();
 395   DDV(extent, 6);
 396   //  std::vector<int> extend;
 397
 398
 399   // Loop on slice
 400   uint n = image->GetDimensions()[2];
 401   DD(n);
 402   DD(mListOfContours.size());
 403   mListOfContours.resize(n); /// ???FIXME
 404   DD(mListOfContours.size());
 405   std::vector<vtkSmartPointer<vtkPolyData> > contours;
 406   for(uint i=0; i<n; i++) {
 407     DD(i);
 408
 409     // FIXME     vtkDiscreteMarchingCubes INSTEAD
 410
 411
 412     vtkSmartPointer<vtkMarchingSquares> squares = vtkSmartPointer<vtkMarchingSquares>::New();
 413     squares->SetInput(image);
 414     squares->SetImageRange(extent[0], extent[1], extent[2], extent[3], i, i);
 415     squares->SetValue(1, 1.0);
 416     squares->Update();
 417     DD(squares->GetNumberOfContours());
 418
 419     //clitk::DicomRT_Contour * contour = new clitk::DicomRT_Contour();
 420     //mListOfContours[i]->SetMesh(squares->GetOutput());
 421
 422
 423     vtkSmartPointer<vtkPolyData> m = squares->GetOutput();
 424     contours.push_back(m);
 425
 426     /*
 427     // Clip to the current slice
 428     extent[4] = extent[5] = image->GetOrigin()[2]+i*image->GetSpacing()[2];
 429     DDV(extent, 6);
 430     // Prepare the marching squares
 431     vtkSmartPointer<vtkMarchingSquares> squares = vtkSmartPointer<vtkMarchingSquares>::New();
 432     clipper->SetOutputWholeExtent(extent[0],extent[1],extent[2],
 433                                 extent[3],extent[4],extent[5]);
 434
 435     squares->SetInput(clipper->GetOutput());
 436     squares->Update();
 437     DD(squares->GetNumberOfContours());
 438     mListOfContours[i]->SetMesh(squares->GetOutput());
 439     */
 440   }
 441   DD("done");
 442
 443   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
 444   for(unsigned int i=0; i<n; i++) {
 445     append->AddInput(contours[i]);
 446   }
 447   append->Update();
 448
 449   mMesh = vtkSmartPointer<vtkPolyData>::New();
 450   mMesh->DeepCopy(append->GetOutput());
 451
 452   // Write vtk
 453   vtkPolyDataWriter * w = vtkPolyDataWriter::New();
 454   w->SetInput(mMesh);
 455   w->SetFileName("toto.vtk");
 456   w->Write();
 457
 458   DD("done");
 459 }
 460 //--------------------------------------------------------------------
 461
 462
 463 //--------------------------------------------------------------------
 464 #if CLITK_USE_SYSTEM_GDCM == 1
 465 void clitk::DicomRT_ROI::Read(vtkSmartPointer<vtkGDCMPolyDataReader> & reader, int roiindex)
 466 {
 467   vtkRTStructSetProperties * p = reader->GetRTStructSetProperties();
 468
 469   mName = p->GetStructureSetROIName(roiindex);
 470   mNumber = p->GetStructureSetROINumber(roiindex);
 471   //mColor = //FIXME !!
 472   SetDicomUptodateFlag(true);
 473   // Get the contour
 474   mMesh =  reader->GetOutput(roiindex);
 475   DicomRT_Contour::Pointer c = DicomRT_Contour::New();
 476   c->SetMesh(mMesh); // FIXME no GetZ, not GetPoints
 477   mMeshIsUpToDate = true;
 478   mListOfContours.push_back(c);
 479 }
 480 #endif
 481 //--------------------------------------------------------------------
 482