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     points[i*3+1] = p[2]-0.5;
  86   }
  87
  88   // Get attribute
  89   gdcm::Attribute<0x3006,0x0050> at;
  90   gdcm::Tag tcontourdata(0x3006,0x0050);
  91   gdcm::DataElement contourdata = nestedds.GetDataElement( tcontourdata );
  92   at.SetFromDataElement( contourdata );
  93
  94   // Set attribute
  95   at.SetValues(&points[0], points.size());
  96   DD(at.GetValues()[0]);
  97
  98   DD("replace");
  99   nestedds.Replace(at.GetAsDataElement());
 100   DD("done");
 101
 102   // Change Number of points [Number of Contour Points]
 103   numcontpoints.SetValue(nb);
 104   nestedds.Replace(numcontpoints.GetAsDataElement());
 105
 106   // Test
 107   gdcm::DataElement aa = nestedds.GetDataElement( tcontourdata );
 108   at.SetFromDataElement( aa );
 109   const double* bb = at.GetValues();
 110   DD(bb[0]);
 111
 112 }
 113 #endif
 114 //--------------------------------------------------------------------
 115
 116
 117 //--------------------------------------------------------------------
 118 #if GDCM_MAJOR_VERSION >= 2
 119 bool clitk::DicomRT_Contour::Read(gdcm::Item * item)
 120 {
 121   mItem = item;
 122
 123   const gdcm::DataSet& nestedds2 = item->GetNestedDataSet();
 124
 125   // Contour type [Contour Geometric Type]
 126   gdcm::Attribute<0x3006,0x0042> contgeotype;
 127   contgeotype.SetFromDataSet( nestedds2 );
 128
 129   if (contgeotype.GetValue() != "CLOSED_PLANAR " && contgeotype.GetValue() != "POINT ") { ///WARNING to the space after the name ...
 130     //std::cerr << "Skip this contour : type=" << mType << std::endl;
 131     return false;
 132   }
 133   if (contgeotype.GetValue() == "POINT ") {
 134     std::cerr << "Warning: POINT type not fully supported. (don't use GetMesh() with this!)"
 135       << std::endl;
 136   }
 137
 138   gdcm::Attribute<0x3006,0x0046> numcontpoints;
 139   numcontpoints.SetFromDataSet( nestedds2 );
 140   // Number of points [Number of Contour Points]
 141   mNbOfPoints = numcontpoints.GetValue();
 142   // DD(mNbOfPoints);
 143
 144   gdcm::Attribute<0x3006,0x0050> at;
 145   gdcm::Tag tcontourdata(0x3006,0x0050);
 146   const gdcm::DataElement & contourdata = nestedds2.GetDataElement( tcontourdata );
 147   at.SetFromDataElement( contourdata );
 148   const double* points = at.GetValues();
 149   //  unsigned int npts = at.GetNumberOfValues() / 3;
 150   assert(at.GetNumberOfValues() == static_cast<unsigned int>(mNbOfPoints)*3);
 151
 152   // Organize values
 153   mData = vtkSmartPointer<vtkPoints>::New();
 154   mData->SetDataTypeToDouble();
 155   mData->SetNumberOfPoints(mNbOfPoints);
 156   for(unsigned int i=0; i<mNbOfPoints; i++) {
 157     double p[3];
 158     p[0] = points[i*3];
 159     p[1] = points[i*3+1];
 160     p[2] = points[i*3+2]+0.5;
 161     mData->SetPoint(i, p);
 162     if (mZ == -1) mZ = p[2];
 163     if (std::fabs(p[2] - mZ) > mTolerance) {
 164       DD(i);
 165       DD(p[2]);
 166       DD(mZ);
 167       std::cout << "ERROR ! contour not in the same slice" << std::endl;
 168       assert(p[2] == mZ);
 169     }
 170   }
 171
 172   return true;
 173
 174 }
 175 #else
 176 bool clitk::DicomRT_Contour::Read(gdcm::SQItem * item)
 177 {
 178
 179   // Contour type [Contour Geometric Type]
 180   mType = item->GetEntryValue(0x3006,0x0042);
 181   // DD(mType);
 182   if (mType != "CLOSED_PLANAR " && mType != "POINT ") { ///WARNING to the space after the name ...
 183     //std::cerr << "Skip this contour : type=" << mType << std::endl;
 184     return false;
 185   }
 186   if (mType == "POINT ") {
 187     std::cerr << "Warning: POINT type not fully supported. (don't use GetMesh() with this!)"
 188       << std::endl;
 189   }
 190
 191   // Number of points [Number of Contour Points]
 192   mNbOfPoints = parse_value<int>(item->GetEntryValue(0x3006,0x0046));
 193   // DD(mNbOfPoints);
 194
 195   // Read values [Contour Data]
 196   std::vector<float> points = parse_string<float>(item->GetEntryValue(0x3006,0x0050),'\\');
 197   assert(points.size() == static_cast<unsigned int>(mNbOfPoints)*3);
 198
 199   // Organize values
 200   mData = vtkSmartPointer<vtkPoints>::New();
 201   mData->SetDataTypeToDouble();
 202   mData->SetNumberOfPoints(mNbOfPoints);
 203   for(unsigned int i=0; i<mNbOfPoints; i++) {
 204     double p[3];
 205     p[0] = points[i*3];
 206     p[1] = points[i*3+1];
 207     p[2] = points[i*3+2]+0.5;
 208     mData->SetPoint(i, p);
 209     if (mZ == -1) mZ = p[2];
 210     if (std::fabs(p[2] - mZ) > mTolerance) {
 211       DD(i);
 212       DD(p[2]);
 213       DD(mZ);
 214       std::cout << "ERROR ! contour not in the same slice" << std::endl;
 215       assert(p[2] == mZ);
 216     }
 217   }
 218
 219   return true;
 220 }
 221 #endif
 222 //--------------------------------------------------------------------
 223
 224
 225 //--------------------------------------------------------------------
 226 vtkPolyData * clitk::DicomRT_Contour::GetMesh()
 227 {
 228   if (!mMeshIsUpToDate) {
 229     ComputeMeshFromDataPoints();
 230   }
 231   return mMesh;
 232 }
 233 //--------------------------------------------------------------------
 234
 235
 236 //--------------------------------------------------------------------
 237 void clitk::DicomRT_Contour::SetMesh(vtkPolyData * mesh)
 238 {
 239   mMesh = mesh;
 240 }
 241 //--------------------------------------------------------------------
 242
 243
 244 //--------------------------------------------------------------------
 245 void clitk::DicomRT_Contour::SetTransformMatrix(vtkMatrix4x4* matrix)
 246 {
 247   mTransformMatrix = matrix;
 248 }
 249 //--------------------------------------------------------------------
 250 //--------------------------------------------------------------------
 251 double clitk::DicomRT_Contour::GetTolerance()
 252 {
 253   return mTolerance;
 254 }
 255 //--------------------------------------------------------------------
 256 void clitk::DicomRT_Contour::SetTolerance(double tol)
 257 {
 258   mTolerance = tol;
 259 }
 260 //--------------------------------------------------------------------
 261 //--------------------------------------------------------------------
 262 void clitk::DicomRT_Contour::ComputeMeshFromDataPoints()
 263 {
 264 //  DD("ComputeMesh Contour");
 265   mMesh = vtkSmartPointer<vtkPolyData>::New();
 266   mMesh->Allocate(); //for cell structures
 267   mPoints = vtkSmartPointer<vtkPoints>::New();
 268   mMesh->SetPoints(mPoints);
 269   vtkIdType ids[2];
 270   for (unsigned int idx=0 ; idx<mNbOfPoints ; idx++) {
 271     //double pointIn[4];
 272     //for (unsigned int j=0 ; j<3; ++j)
 273     //  pointIn[j] = mData->GetPoint(idx)[j];
 274     //pointIn[3] = 1.0;
 275     /*double pointOut[4];
 276     mTransformMatrix->MultiplyPoint(pointIn, pointOut);
 277     std::cout << pointOut[0] << " " << pointOut[1] << " " << pointOut[2] << " " << pointOut[3] << std::endl;
 278     mMesh->GetPoints()->InsertNextPoint(pointOut[0],
 279                                         pointOut[1],
 280                                         pointOut[2]);*/
 281     mMesh->GetPoints()->InsertNextPoint(mData->GetPoint(idx)[0],
 282                                         mData->GetPoint(idx)[1],
 283                                         mData->GetPoint(idx)[2]);
 284     //std::cout << mData->GetPoint(idx)[0] << " " << mData->GetPoint(idx)[1] << " " << mData->GetPoint(idx)[2] << std::endl;
 285     ids[0]=idx;
 286     ids[1]=(ids[0]+1) % mNbOfPoints; //0-1,1-2,...,n-1-0
 287     mMesh->GetLines()->InsertNextCell(2,ids);
 288   }
 289   //std::cout << std::endl;
 290   mMeshIsUpToDate = true;
 291 }
 292 //--------------------------------------------------------------------
 293
 294
 295 //--------------------------------------------------------------------
 296 void clitk::DicomRT_Contour::ComputeDataPointsFromMesh()
 297 {
 298   DD("ComputeDataPointsFromMesh");
 299
 300
 301   /*todo
 302
 303   mMesh = vtkSmartPointer<vtkPolyData>::New();
 304   mMesh->Allocate(); //for cell structures
 305   mPoints = vtkSmartPointer<vtkPoints>::New();
 306   mMesh->SetPoints(mPoints);
 307   vtkIdType ids[2];
 308   for (unsigned int idx=0 ; idx<mNbOfPoints ; idx++) {
 309     mMesh->GetPoints()->InsertNextPoint(mData->GetPoint(idx)[0],
 310                                         mData->GetPoint(idx)[1],
 311                                         mData->GetPoint(idx)[2]);
 312     ids[0]=idx;
 313     ids[1]=(ids[0]+1) % mNbOfPoints; //0-1,1-2,...,n-1-0
 314     mMesh->GetLines()->InsertNextCell(2,ids);
 315   }
 316   mMeshIsUpToDate = true;
 317 */
 318 }
 319 //--------------------------------------------------------------------