common/clitkDicomRT_Contour.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_Contour.h"
  21 #include <vtkCellArray.h>
  22
  23 #if GDCM_MAJOR_VERSION == 2
  24 #include "gdcmAttribute.h"
  25 #include "gdcmItem.h"
  26 #endif
  27
  28 //--------------------------------------------------------------------
  29 clitk::DicomRT_Contour::DicomRT_Contour()
  30 {
  31   mMeshIsUpToDate = false;
  32   mNbOfPoints = 0;
  33   mZ = -1;
  34 }
  35 //--------------------------------------------------------------------
  36
  37
  38 //--------------------------------------------------------------------
  39 clitk::DicomRT_Contour::~DicomRT_Contour()
  40 {
  41
  42 }
  43 //--------------------------------------------------------------------
  44
  45
  46
  47 //--------------------------------------------------------------------
  48 void clitk::DicomRT_Contour::Print(std::ostream & os) const
  49 {
  50   DD("TODO : print Contours");
  51 }
  52 //--------------------------------------------------------------------
  53
  54
  55 //--------------------------------------------------------------------
  56 #if GDCM_MAJOR_VERSION == 2
  57 void clitk::DicomRT_Contour::UpdateDicomItem()
  58 {
  59   DD("DicomRT_Contour::UpdateDicomItem");
  60
  61   gdcm::DataSet & nestedds = mItem->GetNestedDataSet();
  62
  63   //NON CONSIDER CONTOUR ITEM NOT SEQ ITEM ?
  64
  65   // Contour type [Contour Geometric Type]
  66   gdcm::Attribute<0x3006,0x0042> contgeotype;
  67   contgeotype.SetFromDataSet( nestedds );
  68
  69   // Number of points [Number of Contour Points]
  70   gdcm::Attribute<0x3006,0x0046> numcontpoints;
  71   numcontpoints.SetFromDataSet( nestedds );
  72   DD(mNbOfPoints);
  73   mNbOfPoints = numcontpoints.GetValue();
  74   DD(mNbOfPoints);
  75
  76   // Contour dicom values from DataPoints
  77   //ComputeDataPointsFromMesh();
  78   uint nb = mData->GetNumberOfPoints();
  79   std::vector<double> points;
  80   points.resize(mNbOfPoints*3);
  81   for(unsigned int i=0; i<nb; i++) {
  82     double * p = mData->GetPoint(i);
  83     points[i*3] = p[0];
  84     points[i*3+1] = p[1];
  85 #if VTK_MAJOR_VERSION <= 5
  86     points[i*3+1] = p[2];
  87 #else
  88     points[i*3+1] = p[2]-0.5;
  89 #endif
  90   }
  91
  92   // Get attribute
  93   gdcm::Attribute<0x3006,0x0050> at;
  94   gdcm::Tag tcontourdata(0x3006,0x0050);
  95   gdcm::DataElement contourdata = nestedds.GetDataElement( tcontourdata );
  96   at.SetFromDataElement( contourdata );
  97
  98   // Set attribute
  99   at.SetValues(&points[0], points.size(), false);
 100   DD(at.GetValues()[0]);
 101
 102   DD("replace");
 103   nestedds.Replace(at.GetAsDataElement());
 104   DD("done");
 105
 106   // Change Number of points [Number of Contour Points]
 107   numcontpoints.SetValue(nb);
 108   nestedds.Replace(numcontpoints.GetAsDataElement());
 109
 110   // Test
 111   gdcm::DataElement aa = nestedds.GetDataElement( tcontourdata );
 112   at.SetFromDataElement( aa );
 113   const double* bb = at.GetValues();
 114   DD(bb[0]);
 115
 116 }
 117 #endif
 118 //--------------------------------------------------------------------
 119
 120
 121 //--------------------------------------------------------------------
 122 #if GDCM_MAJOR_VERSION == 2
 123 bool clitk::DicomRT_Contour::Read(gdcm::Item * item)
 124 {
 125   mItem = item;
 126
 127   const gdcm::DataSet& nestedds2 = item->GetNestedDataSet();
 128
 129   // Contour type [Contour Geometric Type]
 130   gdcm::Attribute<0x3006,0x0042> contgeotype;
 131   contgeotype.SetFromDataSet( nestedds2 );
 132
 133   if (contgeotype.GetValue() != "CLOSED_PLANAR " && contgeotype.GetValue() != "POINT ") { ///WARNING to the space after the name ...
 134     //std::cerr << "Skip this contour : type=" << mType << std::endl;
 135     return false;
 136   }
 137   if (contgeotype.GetValue() == "POINT ") {
 138     std::cerr << "Warning: POINT type not fully supported. (don't use GetMesh() with this!)"
 139       << std::endl;
 140   }
 141
 142   gdcm::Attribute<0x3006,0x0046> numcontpoints;
 143   numcontpoints.SetFromDataSet( nestedds2 );
 144   // Number of points [Number of Contour Points]
 145   mNbOfPoints = numcontpoints.GetValue();
 146   // DD(mNbOfPoints);
 147
 148   gdcm::Attribute<0x3006,0x0050> at;
 149   gdcm::Tag tcontourdata(0x3006,0x0050);
 150   const gdcm::DataElement & contourdata = nestedds2.GetDataElement( tcontourdata );
 151   at.SetFromDataElement( contourdata );
 152   const double* points = at.GetValues();
 153   //  unsigned int npts = at.GetNumberOfValues() / 3;
 154   assert(at.GetNumberOfValues() == static_cast<unsigned int>(mNbOfPoints)*3);
 155
 156   // Organize values
 157   mData = vtkSmartPointer<vtkPoints>::New();
 158   mData->SetDataTypeToDouble();
 159   mData->SetNumberOfPoints(mNbOfPoints);
 160   for(unsigned int i=0; i<mNbOfPoints; i++) {
 161     double p[3];
 162     p[0] = points[i*3];
 163     p[1] = points[i*3+1];
 164 #if VTK_MAJOR_VERSION <= 5
 165     p[2] = points[i*3+2];
 166 #else
 167     p[2] = points[i*3+2]+0.5;
 168 #endif
 169     mData->SetPoint(i, p);
 170     if (mZ == -1) mZ = p[2];
 171     if (p[2] != mZ) {
 172       DD(i);
 173       DD(p[2]);
 174       DD(mZ);
 175       std::cout << "ERROR ! contour not in the same slice" << std::endl;
 176       assert(p[2] == mZ);
 177     }
 178   }
 179
 180   return true;
 181
 182 }
 183 #else
 184 bool clitk::DicomRT_Contour::Read(gdcm::SQItem * item)
 185 {
 186
 187   // Contour type [Contour Geometric Type]
 188   mType = item->GetEntryValue(0x3006,0x0042);
 189   // DD(mType);
 190   if (mType != "CLOSED_PLANAR " && mType != "POINT ") { ///WARNING to the space after the name ...
 191     //std::cerr << "Skip this contour : type=" << mType << std::endl;
 192     return false;
 193   }
 194   if (mType == "POINT ") {
 195     std::cerr << "Warning: POINT type not fully supported. (don't use GetMesh() with this!)"
 196       << std::endl;
 197   }
 198
 199   // Number of points [Number of Contour Points]
 200   mNbOfPoints = parse_value<int>(item->GetEntryValue(0x3006,0x0046));
 201   // DD(mNbOfPoints);
 202
 203   // Read values [Contour Data]
 204   std::vector<float> points = parse_string<float>(item->GetEntryValue(0x3006,0x0050),'\\');
 205   assert(points.size() == static_cast<unsigned int>(mNbOfPoints)*3);
 206
 207   // Organize values
 208   mData = vtkSmartPointer<vtkPoints>::New();
 209   mData->SetDataTypeToDouble();
 210   mData->SetNumberOfPoints(mNbOfPoints);
 211   for(unsigned int i=0; i<mNbOfPoints; i++) {
 212     double p[3];
 213     p[0] = points[i*3];
 214     p[1] = points[i*3+1];
 215 #if VTK_MAJOR_VERSION <= 5
 216     p[2] = points[i*3+2];
 217 #else
 218     p[2] = points[i*3+2]+0.5;
 219 #endif
 220     mData->SetPoint(i, p);
 221     if (mZ == -1) mZ = p[2];
 222     if (p[2] != mZ) {
 223       DD(i);
 224       DD(p[2]);
 225       DD(mZ);
 226       std::cout << "ERROR ! contour not in the same slice" << std::endl;
 227       assert(p[2] == mZ);
 228     }
 229   }
 230
 231   return true;
 232 }
 233 #endif
 234 //--------------------------------------------------------------------
 235
 236
 237 //--------------------------------------------------------------------
 238 vtkPolyData * clitk::DicomRT_Contour::GetMesh()
 239 {
 240   if (!mMeshIsUpToDate) {
 241     ComputeMeshFromDataPoints();
 242   }
 243   return mMesh;
 244 }
 245 //--------------------------------------------------------------------
 246
 247
 248 //--------------------------------------------------------------------
 249 void clitk::DicomRT_Contour::SetMesh(vtkPolyData * mesh)
 250 {
 251   mMesh = mesh;
 252 }
 253 //--------------------------------------------------------------------
 254
 255
 256 //--------------------------------------------------------------------
 257 void clitk::DicomRT_Contour::ComputeMeshFromDataPoints()
 258 {
 259 //  DD("ComputeMesh Contour");
 260   mMesh = vtkSmartPointer<vtkPolyData>::New();
 261   mMesh->Allocate(); //for cell structures
 262   mPoints = vtkSmartPointer<vtkPoints>::New();
 263   mMesh->SetPoints(mPoints);
 264   vtkIdType ids[2];
 265   for (unsigned int idx=0 ; idx<mNbOfPoints ; idx++) {
 266     mMesh->GetPoints()->InsertNextPoint(mData->GetPoint(idx)[0],
 267                                         mData->GetPoint(idx)[1],
 268                                         mData->GetPoint(idx)[2]);
 269     ids[0]=idx;
 270     ids[1]=(ids[0]+1) % mNbOfPoints; //0-1,1-2,...,n-1-0
 271     mMesh->GetLines()->InsertNextCell(2,ids);
 272   }
 273   mMeshIsUpToDate = true;
 274 }
 275 //--------------------------------------------------------------------
 276
 277
 278 //--------------------------------------------------------------------
 279 void clitk::DicomRT_Contour::ComputeDataPointsFromMesh()
 280 {
 281   DD("ComputeDataPointsFromMesh");
 282
 283
 284   /*todo
 285
 286   mMesh = vtkSmartPointer<vtkPolyData>::New();
 287   mMesh->Allocate(); //for cell structures
 288   mPoints = vtkSmartPointer<vtkPoints>::New();
 289   mMesh->SetPoints(mPoints);
 290   vtkIdType ids[2];
 291   for (unsigned int idx=0 ; idx<mNbOfPoints ; idx++) {
 292     mMesh->GetPoints()->InsertNextPoint(mData->GetPoint(idx)[0],
 293                                         mData->GetPoint(idx)[1],
 294                                         mData->GetPoint(idx)[2]);
 295     ids[0]=idx;
 296     ids[1]=(ids[0]+1) % mNbOfPoints; //0-1,1-2,...,n-1-0
 297     mMesh->GetLines()->InsertNextCell(2,ids);
 298   }
 299   mMeshIsUpToDate = true;
 300 */
 301 }
 302 //--------------------------------------------------------------------