tools/clitkDicomRTPlan2Gate.cxx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef CLITKDICOMRTPLAN2GATE_CXX
  19 #define CLITKDICOMRTPLAN2GATE_CXX
  20
  21 #include "clitkDicomRTPlan2Gate_ggo.h"
  22
  23 #include "clitkIO.h"
  24 #include "clitkCommon.h"
  25
  26 #include "gdcmFile.h"
  27 #include "gdcmReader.h"
  28 #include "gdcmAttribute.h"
  29
  30 int main(int argc, char * argv[])
  31 {
  32
  33   // init command line
  34   GGO(clitkDicomRTPlan2Gate, args_info);
  35   CLITK_INIT;
  36   char* input_file = args_info.input_arg;
  37   char* output_file = args_info.output_arg;
  38
  39   gdcm::Reader reader;
  40   reader.SetFileName(input_file);
  41   if( !reader.Read() )
  42   {
  43     std::cerr << "Failed to read: " << input_file << std::endl;
  44     return 1;
  45   }
  46
  47   const gdcm::DataSet & ds = reader.GetFile().GetDataSet();
  48
  49   //Check file type
  50   gdcm::MediaStorage ms;
  51   ms.SetFromFile(reader.GetFile());
  52   if( ms != gdcm::MediaStorage::RTIonPlanStorage )
  53   {
  54     std::cerr << "The file <" << input_file << "> is not a RT Ion Plan." << std::endl;
  55     exit(0);
  56   }
  57
  58   std::ofstream output;
  59   output.open(output_file);
  60   output << "#TREATMENT-PLAN-DESCRIPTION" << "\n";
  61
  62   gdcm::Attribute<0x300a,0x2> label;
  63   label.SetFromDataSet(ds);
  64   gdcm::Attribute<0x300a,0x3> planName;
  65   planName.SetFromDataSet(ds);
  66
  67   output << "#PlanName\n"<< planName.GetValue() << "\n";
  68
  69   gdcm::Tag tFractionGroupSQ(0x300a,0x70);
  70   gdcm::Attribute<0x300a,0x71> fractionID;
  71   gdcm::Attribute<0x300a,0x78> numberOfFractions;
  72   gdcm::Attribute<0x300a,0x80> numberOfBeams;
  73
  74   gdcm::Tag tReferencedBeamSQ(0x300c,0x4);
  75   gdcm::Attribute<0x300a,0x86> meterSet;
  76   gdcm::Attribute<0x300c,0x6> beamID;
  77
  78   const gdcm::DataElement &FractionGroupSQ = ds.GetDataElement( tFractionGroupSQ );
  79   gdcm::SmartPointer<gdcm::SequenceOfItems> sqf = FractionGroupSQ.GetValueAsSQ();
  80   gdcm::Item &ReferencedBeamItem = sqf->GetItem(1);
  81   gdcm::DataSet &ReferencedBeamData = ReferencedBeamItem.GetNestedDataSet();
  82   numberOfFractions.SetFromDataSet(ReferencedBeamData);
  83
  84   output << "#NumberOfFractions\n"<< numberOfFractions.GetValue() << "\n";
  85
  86   for( unsigned int fidx = 1; fidx <= sqf->GetNumberOfItems(); ++fidx )
  87   {
  88     gdcm::Item &ReferencedBeamItem = sqf->GetItem(fidx);
  89     gdcm::DataSet &ReferencedBeamData = ReferencedBeamItem.GetNestedDataSet();
  90
  91     fractionID.SetFromDataSet(ReferencedBeamData);
  92     numberOfFractions.SetFromDataSet(ReferencedBeamData);
  93     numberOfBeams.SetFromDataSet(ReferencedBeamData);
  94     output << "##FractionID\n"<<  fractionID.GetValue() << "\n";
  95     output << "##NumberOfFields\n"<<  numberOfBeams.GetValue() << "\n";
  96
  97     const gdcm::DataElement &ReferencedBeamSQ = ReferencedBeamData.GetDataElement(tReferencedBeamSQ);
  98     gdcm::SmartPointer<gdcm::SequenceOfItems> sqb = ReferencedBeamSQ.GetValueAsSQ();
  99     for( unsigned int bidx = 1; bidx <= sqb->GetNumberOfItems(); ++bidx )
 100     {
 101       gdcm::Item &ReferencedBeamGroupItem = sqb->GetItem(bidx);
 102       gdcm::DataSet &ReferencedBeamGroupData = ReferencedBeamGroupItem.GetNestedDataSet();
 103
 104       meterSet.SetFromDataSet(ReferencedBeamGroupData);
 105       beamID.SetFromDataSet(ReferencedBeamGroupData);
 106       output << "###FieldsID\n"<<  beamID.GetValue() << "\n";
 107     }
 108   }
 109
 110   gdcm::Tag tPatientSetupSQ(0x300a,0x180);
 111   gdcm::Attribute<0x0018,0x5100> position;
 112   gdcm::Attribute<0x300a,0x1b0> setupTechnique;
 113   gdcm::Attribute<0x300a,0x182> setupID;
 114
 115   const gdcm::DataElement &PatientSetupSQ = ds.GetDataElement( tPatientSetupSQ );
 116   gdcm::SmartPointer<gdcm::SequenceOfItems> sqs = PatientSetupSQ.GetValueAsSQ();
 117   for( unsigned int sidx = 1; sidx <= sqs->GetNumberOfItems(); ++sidx )
 118   {
 119     gdcm::Item &PatientSetupItem = sqs->GetItem(sidx);
 120     gdcm::DataSet &PatientSetupData = PatientSetupItem.GetNestedDataSet();
 121
 122     position.SetFromDataSet(PatientSetupData);
 123     setupTechnique.SetFromDataSet(PatientSetupData);
 124     setupID.SetFromDataSet(PatientSetupData);
 125   }
 126
 127   gdcm::Tag tIonBeamSQ(0x300a,0x3a2);
 128   gdcm::Attribute<0x300a,0xb2> machine;
 129   gdcm::Attribute<0x300a,0xc6> radiationType;
 130   gdcm::Attribute<0x300a,0xce> treatmentType;
 131   gdcm::Attribute<0x300a,0xc0> ionBeamID;
 132   gdcm::Attribute<0x300a,0x10e> finalCumulativeMeterSetWeight;
 133   gdcm::Attribute<0x300a,0x110> numberOfControlPoints;
 134   gdcm::Attribute<0x300c,0x6a> patientSetupID;
 135   gdcm::Attribute<0x300a,0xd0> numberOfWedges;
 136   gdcm::Attribute<0x300a,0xe0> numberOfCompensators;
 137   gdcm::Attribute<0x300a,0xed> numberOfBoli;
 138   gdcm::Attribute<0x300a,0xf0> numberOfBlocks;
 139   gdcm::Attribute<0x300a,0x312> numberOfRangeShifters;
 140   gdcm::Attribute<0x300a,0x330> numberOfLateralSpreadingDevices;
 141   gdcm::Attribute<0x300a,0x340> numberOfRangeModulators;
 142
 143   gdcm::Tag tSnoutSQ(0x300a,0x30c);
 144   gdcm::Attribute<0x300a,0x30f> snoutID;
 145
 146   gdcm::Tag tControlPointSQ(0x300a,0x3a8);
 147   gdcm::Attribute<0x300a,0x15> energyUnit;
 148   gdcm::Attribute<0x300a,0x114> energyValue;
 149   gdcm::Attribute<0x300a,0x112> controlPointIndex;
 150   gdcm::Attribute<0x300a,0x134> cumulativeMetersetWeight;
 151   gdcm::Attribute<0x300a,0x392> numberOfScannedSpots;
 152   gdcm::Attribute<0x300a,0x390> spotTunnedID;
 153
 154   gdcm::Tag tSpotPositionsData(0x300a,0x394);
 155   gdcm::Attribute<0x300a,0x394> spotPositions;
 156   gdcm::Tag tSpotWeightsData(0x300a,0x396);
 157   gdcm::Attribute<0x300a,0x396> spotWeights;
 158
 159   gdcm::Attribute<0x300a,0x11e> gantryAngle;
 160   gdcm::Attribute<0x300a,0x122> patientSupportAngle;
 161   gdcm::Attribute<0x300a,0x12c> isocenterPosition;
 162
 163   int totalCumMeterSet = 0;
 164
 165   const gdcm::DataElement &IonBeamSQ = ds.GetDataElement( tIonBeamSQ );
 166   gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = IonBeamSQ.GetValueAsSQ();
 167
 168   for( unsigned int iidx = 1; iidx <= sqi->GetNumberOfItems(); ++iidx )
 169   {
 170     gdcm::Item &IonBeamItem = sqi->GetItem(iidx);
 171     gdcm::DataSet &IonBeamData = IonBeamItem.GetNestedDataSet();
 172     finalCumulativeMeterSetWeight.SetFromDataSet(IonBeamData);
 173     totalCumMeterSet += finalCumulativeMeterSetWeight.GetValue();
 174   }
 175
 176   output << "#TotalMetersetWeightOfAllFields\n" <<  totalCumMeterSet << "\n";
 177
 178   for( unsigned int iidx = 1; iidx <= sqi->GetNumberOfItems(); ++iidx )
 179   {
 180     gdcm::Item &IonBeamItem = sqi->GetItem(iidx);
 181     gdcm::DataSet &IonBeamData = IonBeamItem.GetNestedDataSet();
 182
 183     machine.SetFromDataSet(IonBeamData);
 184     radiationType.SetFromDataSet(IonBeamData);
 185     treatmentType.SetFromDataSet(IonBeamData);
 186     ionBeamID.SetFromDataSet(IonBeamData);
 187     finalCumulativeMeterSetWeight.SetFromDataSet(IonBeamData);
 188     numberOfControlPoints.SetFromDataSet(IonBeamData);
 189     patientSetupID.SetFromDataSet(IonBeamData);
 190     numberOfWedges.SetFromDataSet(IonBeamData);
 191     numberOfCompensators.SetFromDataSet(IonBeamData);
 192     numberOfBoli.SetFromDataSet(IonBeamData);
 193     numberOfBlocks.SetFromDataSet(IonBeamData);
 194     numberOfRangeShifters.SetFromDataSet(IonBeamData);
 195     numberOfLateralSpreadingDevices.SetFromDataSet(IonBeamData);
 196     numberOfRangeModulators.SetFromDataSet(IonBeamData);
 197
 198     const gdcm::DataElement &SnoutSQ = IonBeamData.GetDataElement(tSnoutSQ);
 199     gdcm::SmartPointer<gdcm::SequenceOfItems> sqsn = SnoutSQ.GetValueAsSQ();
 200     for( unsigned int snidx = 1; snidx <= sqsn->GetNumberOfItems(); ++snidx )
 201     {
 202       gdcm::Item &SnoutItem = sqsn->GetItem(snidx);
 203       gdcm::DataSet &SnoutData = SnoutItem.GetNestedDataSet();
 204
 205       snoutID.SetFromDataSet(SnoutData);
 206     }
 207
 208     const gdcm::DataElement &ControlPointSQ = IonBeamData.GetDataElement(tControlPointSQ);
 209     gdcm::SmartPointer<gdcm::SequenceOfItems> sqcp = ControlPointSQ.GetValueAsSQ();
 210     gdcm::Item &ControlPointItem = sqcp->GetItem(1);
 211     gdcm::DataSet &ControlPointData = ControlPointItem.GetNestedDataSet();
 212
 213     gantryAngle.SetFromDataSet(ControlPointData);
 214     patientSupportAngle.SetFromDataSet(ControlPointData);
 215     isocenterPosition.SetFromDataSet(ControlPointData);
 216
 217     output << "\n#FIELD-DESCRIPTION" << "\n";
 218     output << "###FieldID\n"<< ionBeamID.GetValue() << "\n";
 219     output << "###FinalCumulativeMeterSetWeight\n"<<  finalCumulativeMeterSetWeight.GetValue() << "\n";
 220     output << "###GantryAngle\n" << gantryAngle.GetValue() << "\n";
 221     output << "###PatientSupportAngle\n" << patientSupportAngle.GetValue() << "\n";
 222     output << "###IsocenterPosition\n" << isocenterPosition.GetValue(0);
 223     output << " " << isocenterPosition.GetValue(1);
 224     output << " " << isocenterPosition.GetValue(2) << "\n";
 225     output << "###NumberOfControlPoints\n" << numberOfControlPoints.GetValue() << "\n\n";
 226     output << "#SPOTS-DESCRIPTION" << "\n";
 227
 228     for( unsigned int cpidx = 1; cpidx <= sqcp->GetNumberOfItems(); ++cpidx )
 229     {
 230       gdcm::Item &ControlPointItem = sqcp->GetItem(cpidx);
 231       gdcm::DataSet &ControlPointData = ControlPointItem.GetNestedDataSet();
 232
 233       if ( cpidx == 1 )
 234       {
 235         gantryAngle.SetFromDataSet(ControlPointData);
 236         patientSupportAngle.SetFromDataSet(ControlPointData);
 237         isocenterPosition.SetFromDataSet(ControlPointData);
 238       }
 239
 240       if ( cpidx < sqcp->GetNumberOfItems() )
 241       {
 242         energyUnit.SetFromDataSet(ControlPointData);
 243         energyValue.SetFromDataSet(ControlPointData);
 244         controlPointIndex.SetFromDataSet(ControlPointData);
 245         cumulativeMetersetWeight.SetFromDataSet(ControlPointData);
 246         numberOfScannedSpots.SetFromDataSet(ControlPointData);
 247         spotTunnedID.SetFromDataSet(ControlPointData);
 248
 249         const gdcm::DataElement &spotPositionsData = ControlPointData.GetDataElement(tSpotPositionsData);
 250         spotPositions.SetFromDataElement(spotPositionsData);
 251         const gdcm::DataElement &spotWeightsData = ControlPointData.GetDataElement(tSpotWeightsData);
 252         spotWeights.SetFromDataElement(spotWeightsData);
 253       }
 254
 255       if ( cpidx == sqcp->GetNumberOfItems() )
 256       {
 257         energyUnit.SetValue( "0" );
 258         energyValue.SetValue( 0 );
 259         controlPointIndex.SetFromDataSet(ControlPointData);
 260         cumulativeMetersetWeight.SetFromDataSet(ControlPointData);
 261         numberOfScannedSpots.SetValue( 0 );
 262         spotTunnedID.SetValue( "0" );
 263       }
 264
 265       output << "####ControlPointIndex\n" << controlPointIndex.GetValue() << "\n";
 266       output << "####SpotTunnedID\n" << spotTunnedID.GetValue() << "\n";
 267       output << "####CumulativeMetersetWeight\n" << cumulativeMetersetWeight.GetValue() << "\n";
 268       output << "####Energy (MeV)\n"<< energyValue.GetValue() << "\n";
 269       output << "####NbOfScannedSpots\n" << numberOfScannedSpots.GetValue() << "\n";
 270       output << "####X Y Weight\n";
 271       for ( unsigned int sp=0; sp < numberOfScannedSpots.GetValue(); ++sp)
 272       {
 273         output << spotPositions.GetValue(2*sp) << " " << spotPositions.GetValue(2*sp+1) << " " << spotWeights.GetValue(sp) << "\n";
 274       }
 275
 276     }
 277
 278   }
 279
 280   return 0;
 281 }
 282
 283 #endif