common/clitkDicomRT_StructureSet.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_StructureSet.h"
  21 #include <vtksys/SystemTools.hxx>
  22 #include "gdcmFile.h"
  23
  24 //--------------------------------------------------------------------
  25 clitk::DicomRT_StructureSet::DicomRT_StructureSet()
  26 {
  27   mStudyID = "NoStudyID";
  28   mStudyTime = "NoStudyTime";
  29   mStudyDate = "NoStudyDate";
  30   mLabel = "NoLabel";
  31   mName = "NoName";
  32   mDate = "NoDate";
  33   mTime = "NoTime";
  34   mFile = NULL;
  35 }
  36 //--------------------------------------------------------------------
  37
  38
  39 //--------------------------------------------------------------------
  40 clitk::DicomRT_StructureSet::~DicomRT_StructureSet()
  41 {
  42 }
  43 //--------------------------------------------------------------------
  44
  45
  46 //--------------------------------------------------------------------
  47 const std::string & clitk::DicomRT_StructureSet::GetStudyID() const
  48 {
  49   return mStudyID;
  50 }
  51 //--------------------------------------------------------------------
  52
  53
  54 //--------------------------------------------------------------------
  55 const std::string & clitk::DicomRT_StructureSet::GetStudyTime() const
  56 {
  57   return mStudyTime;
  58 }
  59 //--------------------------------------------------------------------
  60
  61
  62 //--------------------------------------------------------------------
  63 const std::string & clitk::DicomRT_StructureSet::GetStudyDate() const
  64 {
  65   return mStudyDate;
  66 }
  67 //--------------------------------------------------------------------
  68
  69
  70 //--------------------------------------------------------------------
  71 const std::string & clitk::DicomRT_StructureSet::GetLabel() const
  72 {
  73   return mLabel;
  74 }
  75 //--------------------------------------------------------------------
  76
  77
  78 //--------------------------------------------------------------------
  79 const std::string & clitk::DicomRT_StructureSet::GetName() const
  80 {
  81   return mName;
  82 }
  83 //--------------------------------------------------------------------
  84
  85
  86 //--------------------------------------------------------------------
  87 const std::string & clitk::DicomRT_StructureSet::GetDate() const
  88 {
  89   return mDate;
  90 }
  91 //--------------------------------------------------------------------
  92
  93
  94 //--------------------------------------------------------------------
  95 void clitk::DicomRT_StructureSet::SetTransformMatrix(vtkMatrix4x4* matrix)
  96 {
  97   mTransformMatrix = matrix;
  98 }
  99 //--------------------------------------------------------------------
 100
 101
 102 //--------------------------------------------------------------------
 103 const std::string & clitk::DicomRT_StructureSet::GetTime() const
 104 {
 105   return mTime;
 106 }
 107 //--------------------------------------------------------------------
 108
 109
 110 //--------------------------------------------------------------------
 111 // const std::vector<clitk::DicomRT_ROI::Pointer> & clitk::DicomRT_StructureSet::GetListOfROI() const
 112 // {
 113 //   return mListOfROI;
 114 // }
 115 //--------------------------------------------------------------------
 116
 117
 118 //--------------------------------------------------------------------
 119 clitk::DicomRT_ROI* clitk::DicomRT_StructureSet::GetROIFromROINumber(int n)
 120 {
 121   if (mROIs.find(n) == mROIs.end()) {
 122     std::cerr << "No ROI number " << n << std::endl;
 123     return NULL;
 124   }
 125   return mROIs[n];
 126 }
 127 //--------------------------------------------------------------------
 128
 129 //--------------------------------------------------------------------
 130 clitk::DicomRT_ROI* clitk::DicomRT_StructureSet::GetROIFromROIName(const std::string& name)
 131 {
 132   std::map<int, std::string>::iterator it = mMapOfROIName.begin();
 133   int number = -1;
 134   while (it != mMapOfROIName.end() && number == -1) {
 135     if (it->second == name)
 136       number = it->first;
 137     else
 138       it++;
 139   }
 140
 141   if (number == -1) {
 142     std::cerr << "No ROI name " << name << std::endl;
 143     return NULL;
 144   }
 145
 146   return mROIs[number];
 147 }
 148 //--------------------------------------------------------------------
 149 /*
 150 // RP: 08/02/2013
 151 // RegEx version shall be available when C++x11 supports it propely
 152 //
 153 //--------------------------------------------------------------------
 154 clitk::DicomRT_ROI* clitk::DicomRT_StructureSet::GetROIFromROINameRegEx(const std::string& regEx)
 155 {
 156   std::map<int, std::string>::iterator it = mMapOfROIName.begin();
 157   int number = -1;
 158
 159   while (it != mMapOfROIName.end() && number == -1) {
 160     if (std::tr1::regex_match (it->second, std::tr1::regex(regEx)))
 161       number = it->first;
 162     else
 163       it++;
 164   }
 165
 166   if (number == -1) {
 167     std::cerr << "No ROI name " << number << std::endl;
 168     return NULL;
 169   }
 170
 171   return mROIs[number];
 172 }
 173 //--------------------------------------------------------------------
 174 */
 175 //--------------------------------------------------------------------
 176 clitk::DicomRT_ROI* clitk::DicomRT_StructureSet::GetROIFromROINameSubstr(const std::string& s)
 177 {
 178   std::map<int, std::string>::iterator it = mMapOfROIName.begin();
 179   int number = -1;
 180
 181   while (it != mMapOfROIName.end() && number == -1) {
 182     if (it->second.find(s) != std::string::npos)
 183       number = it->first;
 184     else
 185       it++;
 186   }
 187
 188   if (number == -1) {
 189     std::cerr << "No ROI name " << s << std::endl;
 190     return NULL;
 191   }
 192
 193   return mROIs[number];
 194 }
 195 //--------------------------------------------------------------------
 196
 197 //--------------------------------------------------------------------
 198 clitk::DicomRT_StructureSet::ROIMapContainer *
 199 clitk::DicomRT_StructureSet::GetROIsFromROINameSubstr(const std::string& s)
 200 {
 201   static ROIMapContainer rois;
 202   rois.clear();
 203
 204   ROIMapContainer::iterator it = mROIs.begin();
 205   int number = -1;
 206
 207   while (it != mROIs.end()) {
 208     if (it->second->GetName().find(s) != std::string::npos) {
 209       number = it->first;
 210       rois[number] = it->second;
 211     }
 212     it++;
 213   }
 214
 215   if (number == -1) {
 216     std::cerr << "No ROI name " << s << std::endl;
 217     return NULL;
 218   }
 219
 220   return &rois;
 221
 222 }
 223 //--------------------------------------------------------------------
 224
 225 //--------------------------------------------------------------------
 226 void clitk::DicomRT_StructureSet::Print(std::ostream & os) const
 227 {
 228   os << "Study ID      = " << mStudyID << std::endl
 229      << "Study Date    = " << mStudyDate << std::endl
 230      << "Study Time    = " << mStudyTime << std::endl
 231      << "Struct Label  = " << mLabel << std::endl
 232      << "Struct Name   = " << mName << std::endl
 233      << "Struct Time   = " << mTime << std::endl
 234      << "Number of ROI = " << mROIs.size() << std::endl;
 235   for(ROIConstIteratorType iter = mROIs.begin(); iter != mROIs.end(); iter++) {
 236     iter->second->Print(os);
 237   }
 238 }
 239 //--------------------------------------------------------------------
 240
 241
 242 #if GDCM_MAJOR_VERSION >= 2
 243 //--------------------------------------------------------------------
 244 int clitk::DicomRT_StructureSet::ReadROINumber(const gdcm::Item & item)
 245 {
 246   // 0x3006,0x0022 = [ROI Number]
 247   const gdcm::DataSet & nestedds = item.GetNestedDataSet();
 248   gdcm::Attribute<0x3006,0x0022> roinumber;
 249   roinumber.SetFromDataSet( nestedds );
 250   return roinumber.GetValue();
 251 }
 252 //--------------------------------------------------------------------
 253 #endif
 254
 255 //--------------------------------------------------------------------
 256 void clitk::DicomRT_StructureSet::Write(const std::string & filename)
 257 {
 258 #if GDCM_MAJOR_VERSION >= 2
 259
 260   // Assert that the gdcm file is still open (we can write only if it was readed)
 261   if (mFile == NULL) {
 262     //assert(mFile != NULL);
 263     FATAL("Sorry, I can write DICOM only if it was read first from a file with 'Read' function");
 264   }
 265
 266   // Loop and update each ROI
 267   int i=0;
 268   for(ROIIteratorType iter = mROIs.begin(); iter != mROIs.end(); iter++) {
 269     iter->second->UpdateDicomItem();
 270     i++;
 271   }
 272
 273   // Write [ Structure Set ROI Sequence ] = 0x3006,0x0020
 274   gdcm::DataSet & ds = mFile->GetDataSet();
 275   gdcm::Tag tssroisq(0x3006,0x0020);
 276   const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
 277   gdcm::DataElement de(ssroisq);
 278   de.SetValue(*mROIInfoSequenceOfItems);
 279   ds.Replace(de);
 280
 281   // Write [ ROI Contour Sequence ] = 0x3006,0x0039
 282   DD("ici");
 283   gdcm::Tag troicsq(0x3006,0x0039);
 284   const gdcm::DataElement &roicsq = ds.GetDataElement( troicsq );
 285   gdcm::DataElement de2(roicsq);
 286   de2.SetValue(*mROIContoursSequenceOfItems);
 287   ds.Replace(de2);
 288
 289   //DEBUG
 290   gdcm::DataSet & a = mROIContoursSequenceOfItems->GetItem(1).GetNestedDataSet();
 291   gdcm::Tag tcsq(0x3006,0x0040);
 292   const gdcm::DataElement& csq = a.GetDataElement( tcsq );
 293   gdcm::SmartPointer<gdcm::SequenceOfItems> sqi2 = csq.GetValueAsSQ();
 294   gdcm::Item & j = sqi2->GetItem(1);
 295   gdcm::DataSet & b = j.GetNestedDataSet();
 296   gdcm::Attribute<0x3006,0x0050> at;
 297   gdcm::Tag tcontourdata(0x3006,0x0050);
 298   gdcm::DataElement contourdata = b.GetDataElement( tcontourdata );
 299   at.SetFromDataElement( contourdata );
 300   const double* points = at.GetValues();
 301   DD(points[0]);
 302
 303
 304   // Write dicom
 305   gdcm::Writer writer;
 306   //writer.CheckFileMetaInformationOff();
 307   writer.SetFileName(filename.c_str());
 308   writer.SetFile(*mFile);
 309   DD("before write");
 310   writer.Write();
 311   DD("End write");
 312 #else
 313   FATAL("Sorry not compatible with GDCM1, use GDCM2");
 314 #endif
 315 }
 316 //--------------------------------------------------------------------
 317
 318
 319 //--------------------------------------------------------------------
 320 void clitk::DicomRT_StructureSet::Read(const std::string & filename)
 321 {
 322
 323 //Try to avoid to use extern GDCM library
 324
 325   //check the RS file is available before conversion
 326   gdcm::Reader RTreader;
 327   RTreader.SetFileName( filename.c_str() );
 328   if( !RTreader.Read() )
 329   {
 330     std::cout << "Problem reading file: " << filename << std::endl;
 331     return;
 332   }
 333
 334   const gdcm::DataSet& ds = RTreader.GetFile().GetDataSet();
 335
 336   gdcm::MediaStorage ms;
 337   ms.SetFromFile( RTreader.GetFile() );
 338
 339   // (3006,0020) SQ (Sequence with explicit length #=4)      # 370, 1 StructureSetROISequence
 340   gdcm::Tag tssroisq(0x3006,0x0020);
 341   if( !ds.FindDataElement( tssroisq ) )
 342   {
 343     std::cout << "Problem locating 0x3006,0x0020 - Is this a valid RT Struct file?" << std::endl;
 344     return;
 345   }
 346   gdcm::Tag troicsq(0x3006,0x0039);
 347   if( !ds.FindDataElement( troicsq ) )
 348   {
 349     std::cout << "Problem locating 0x3006,0x0039 - Is this a valid RT Struct file?" << std::endl;
 350     return;
 351   }
 352
 353   const gdcm::DataElement &roicsq = ds.GetDataElement( troicsq );
 354
 355   gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = roicsq.GetValueAsSQ();
 356   if( !sqi || !sqi->GetNumberOfItems() )
 357   {
 358     return;
 359   }
 360   const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
 361   gdcm::SmartPointer<gdcm::SequenceOfItems> ssqi = ssroisq.GetValueAsSQ();
 362   if( !ssqi || !ssqi->GetNumberOfItems() )
 363   {
 364     return;
 365   }
 366
 367   // Read global info
 368   gdcm::Attribute<0x20,0x10> studyid;
 369   studyid.SetFromDataSet( ds );
 370   gdcm::Attribute<0x8,0x20> studytime;
 371   studytime.SetFromDataSet( ds );
 372   gdcm::Attribute<0x8,0x30> studydate;
 373   studydate.SetFromDataSet( ds );
 374   gdcm::Attribute<0x3006,0x02> label;
 375   label.SetFromDataSet( ds );
 376   gdcm::Attribute<0x3006,0x04> atname;
 377   atname.SetFromDataSet( ds );
 378   gdcm::Attribute<0x3006,0x09> time;
 379   time.SetFromDataSet( ds );
 380
 381   mStudyID   = studyid.GetValue();
 382   mStudyTime = studytime.GetValue();
 383   mStudyDate = studydate.GetValue();
 384   mLabel     = label.GetValue();
 385   mName      = atname.GetValue();
 386   mTime      = time.GetValue();
 387
 388   // Temporary store the list of items
 389   std::map<int, gdcm::Item*> mMapOfROIInfo;
 390   std::map<int, gdcm::Item*> mMapOfROIContours;
 391
 392
 393   //----------------------------------
 394   // Read all ROI Names and number
 395   // 0x3006,0x0020 = [ Structure Set ROI Sequence ]
 396   //gdcm::Tag tssroisq(0x3006,0x0020);
 397   //const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq );
 398   mROIInfoSequenceOfItems = ssroisq.GetValueAsSQ();
 399   gdcm::SmartPointer<gdcm::SequenceOfItems> & roi_seq = mROIInfoSequenceOfItems;
 400   assert(roi_seq); // TODO error message
 401   for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx)
 402     {
 403     gdcm::Item & item = roi_seq->GetItem( ridx + 1); // Item starts at 1
 404     const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 405
 406     gdcm::Attribute<0x3006,0x26> roiname;
 407     roiname.SetFromDataSet( nestedds );
 408     std::string name = roiname.GetValue(); // 0x3006,0x0026 = [ROI Name]
 409
 410     // 0x3006,0x0022 = [ROI Number]
 411     int nb = ReadROINumber(item);
 412
 413     // Store the item
 414     mMapOfROIInfo[nb] = &item;
 415
 416     // Check if such a number already exist
 417     if (mMapOfROIName.find(nb) != mMapOfROIName.end()) {
 418       std::cerr << "WARNING. A Roi already exist with the number "
 419         << nb << ". I replace." << std::endl;
 420     }
 421     // Add in map
 422     mMapOfROIName[nb] = name;
 423     }
 424
 425   //----------------------------------
 426   // Read all ROI item
 427   // 0x3006,0x0039 = [ ROI Contour Sequence ]
 428   //gdcm::Tag troicsq(0x3006,0x0039);
 429   //const gdcm::DataElement &roicsq = ds.GetDataElement( troicsq );
 430   gdcm::SmartPointer<gdcm::SequenceOfItems> roi_contour_seq = roicsq.GetValueAsSQ();
 431   mROIContoursSequenceOfItems = roi_contour_seq;
 432   assert(roi_contour_seq); // TODO error message
 433   for(unsigned int ridx = 0; ridx < roi_contour_seq->GetNumberOfItems(); ++ridx) {
 434     gdcm::Item & item = roi_contour_seq->GetItem( ridx + 1); // Item starts at 1
 435     // ROI number [Referenced ROI Number]
 436     const gdcm::DataSet& nestedds = item.GetNestedDataSet();
 437     gdcm::Attribute<0x3006,0x0084> referencedroinumber;
 438     referencedroinumber.SetFromDataSet( nestedds );
 439     int nb = referencedroinumber.GetValue();
 440     // Store the item
 441     mMapOfROIContours[nb] = &item;
 442   }
 443
 444   //----------------------------------
 445   // Create the ROIs
 446   for(std::map<int, gdcm::Item*>::iterator i = mMapOfROIInfo.begin(); i != mMapOfROIInfo.end(); i++) {
 447     int nb = i->first;//ReadROINumber(i);//mROIIndex[i];
 448     // Create the roi
 449     mROIs[nb] = DicomRT_ROI::New();
 450     mROIs[nb]->SetTransformMatrix(mTransformMatrix);
 451     mROIs[nb]->Read(mMapOfROIInfo[nb], mMapOfROIContours[nb]);
 452   }
 453
 454   return;
 455
 456 }
 457 //--------------------------------------------------------------------
 458
 459
 460 //--------------------------------------------------------------------
 461 bool clitk::DicomRT_StructureSet::IsDicomRTStruct(const std::string & filename)
 462 {
 463   // Open DICOM
 464 #if GDCM_MAJOR_VERSION >= 2
 465   // Read gdcm file
 466   mReader = new gdcm::Reader;
 467   mReader->SetFileName(filename.c_str());
 468   mReader->Read();
 469   mFile = &(mReader->GetFile());
 470   const gdcm::DataSet & ds = mFile->GetDataSet();
 471
 472   // Check file type
 473   //Verify if the file is a RT-Structure-Set dicom file
 474   gdcm::MediaStorage ms;
 475   ms.SetFromFile(*mFile);
 476   if( ms != gdcm::MediaStorage::RTStructureSetStorage ) return false;
 477
 478   gdcm::Attribute<0x8,0x60> modality;
 479   modality.SetFromDataSet( ds );
 480   if( modality.GetValue() != "RTSTRUCT" ) return false;
 481
 482   return true;
 483
 484   //----------------------------------------------------------------------------------------
 485 #else
 486   mFile = new gdcm::File;
 487   mFile->SetFileName(filename.c_str());
 488   mFile->SetMaxSizeLoadEntry(16384); // Needed ...
 489   mFile->SetLoadMode(gdcm::LD_NOSHADOW); // don't load shadow tags (in order to save memory)
 490   mFile->Load();
 491
 492   // Check file type
 493   //Verify if the file is a RT-Structure-Set dicom file
 494   if (!gdcm::Util::DicomStringEqual(mFile->GetEntryValue(0x0008,0x0016),"1.2.840.10008.5.1.4.1.1.481.3"))
 495     return false;
 496   if (!gdcm::Util::DicomStringEqual(mFile->GetEntryValue(0x0008,0x0060),"RTSTRUCT"))
 497     return false;
 498
 499   return true;
 500
 501 #endif
 502 }
 503 //--------------------------------------------------------------------
 504
 505
 506 //--------------------------------------------------------------------
 507 int clitk::DicomRT_StructureSet::AddBinaryImageAsNewROI(vvImage * im, std::string n)
 508 {
 509   // Search max ROI number
 510   int max = -1;
 511   for(ROIConstIteratorType iter = mROIs.begin(); iter != mROIs.end(); iter++) {
 512     //  for(unsigned int i=0; i<mListOfROI.size(); i++) {
 513     clitk::DicomRT_ROI::Pointer roi = iter->second;
 514     if (roi->GetROINumber() > max)
 515       max = roi->GetROINumber();
 516   }
 517   ++max;
 518
 519   // Compute name
 520   std::ostringstream oss;
 521   oss << vtksys::SystemTools::GetFilenameName(vtksys::SystemTools::GetFilenameWithoutLastExtension(n));
 522   mMapOfROIName[max] = oss.str();
 523
 524   // Set color
 525   std::vector<double> color;
 526   color.push_back(1);
 527   color.push_back(0);
 528   color.push_back(0);
 529
 530   // Create ROI
 531   DicomRT_ROI::Pointer roi = DicomRT_ROI::New();
 532   roi->SetFromBinaryImage(im, max, oss.str(), color, n);
 533   mROIs[max] = roi;
 534   return max;
 535 }
 536 //--------------------------------------------------------------------
 537
 538