# gdcmrtionplan.cxx

    /*=========================================================================

      Program: GDCM (Grassroots DICOM). A DICOM library
      Module:  $URL$

      Copyright (c) 2006-2010 Mathieu Malaterre
      All rights reserved.
      See Copyright.txt or http://gdcm.sourceforge.net/Copyright.html for details.

         This software is distributed WITHOUT ANY WARRANTY; without even
         the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
         PURPOSE.  See the above copyright notice for more information.

    =========================================================================*/
    #include "vtkImageData.h"
    #include "vtkPointData.h"
    #include "vtkPolyData.h"
    #include "vtkProperty.h"
    #include "vtkPolyDataMapper.h"
    #include "vtkActor.h"
    #include "vtkRenderer.h"
    #include "vtkCellArray.h"
    #include "vtkPoints.h"
    #include "vtkDoubleArray.h"
    #include <vtkXMLImageDataWriter.h>
    #include <vtkXMLPolyDataWriter.h>
    #include <vtkRenderWindowInteractor.h>
    #include <vtkImageColorViewer.h>

    #include "gdcmReader.h"
    #include "gdcmAttribute.h"

    /*
     This example is just for fun. We found a RT Ion Plan Storage and simply extracted the viz stuff for VTK

        RTIonPlanStorage, // 1.2.840.10008.5.1.4.1.1.481.8
    */
    int main(int argc, char *argv[])
    {
      if( argc < 3 )
        {
        std::cerr << argv[0] << " filename.dcm outfile.vti\n";
        return 1;
        }
      const char * filename = argv[1];
      const char * outfilename = argv[2];
      const char * outfilename2 = argv[3];

      gdcm::Reader reader;
      reader.SetFileName( filename );
      if( !reader.Read() )
        {
        return 1;
        }

      gdcm::MediaStorage ms;
      ms.SetFromFile( reader.GetFile() );
      if( ms != gdcm::MediaStorage::RTIonPlanStorage )
        {
        return 1;
        }

    /*
    (300a,03a2) SQ                                                    # u/l,1 Ion Beam Sequence
      (fffe,e000) na (Item with undefined length)
        (0008,1040) LO [Test]                                         # 4,1 Institutional Department Name
        (300a,00b2) SH (no value)                                     # 0,1 Treatment Machine Name
        (300a,00b3) CS [MU]                                           # 2,1 Primary Dosimeter Unit
        (300a,00c0) IS [1 ]                                           # 2,1 Beam Number
        (300a,00c2) LO [1 ]                                           # 2,1 Beam Name
        (300a,00c4) CS [STATIC]                                       # 6,1 Beam Type
        (300a,00c6) CS [PROTON]                                       # 6,1 Radiation Type
        (300a,00ce) CS [TREATMENT ]                                   # 10,1 Treatment Delivery Type
        (300a,00d0) IS [0 ]                                           # 2,1 Number of Wedges
        (300a,00e0) IS [1 ]                                           # 2,1 Number of Compensators
        (300a,00ed) IS [0 ]                                           # 2,1 Number of Boli
        (300a,00f0) IS [1 ]                                           # 2,1 Number of Blocks
        (300a,0110) IS [2 ]                                           # 2,1 Number of Control Points
        (300a,02ea) SQ                                                # u/l,1 Ion Range Compensator Sequence
          (fffe,e000) na (Item with undefined length)
            (300a,00e1) SH [lucite]                                   # 6,1 Material ID
            (300a,00e4) IS [1 ]                                       # 2,1 Compensator Number
            (300a,00e5) SH [75hdhe5 ]                                 # 8,1 Compensator ID
            (300a,00e7) IS [35]                                       # 2,1 Compensator Rows
            (300a,00e8) IS [37]                                       # 2,1 Compensator Columns
            (300a,00e9) DS [3.679991\4.249288 ]                       # 18,2 Compensator Pixel Spacing
            (300a,00ea) DS [-76.00\62.50]                             # 12,2 Compensator Position
            (300a,00ec) DS [52.13\52.13\52.13\53.18\54.04\54.04\47.11\40.06\40.06\38.79\34.87\33.28\33.28\33.28\33.28\35.43\35.43\34.54\34.54\34.71\36.10\38.62\44.88\44.88\44.88\45.00\45.00\45.00\45.66\45.66\46.42\39.77\39.77\39.77\39.77\39.77\43.52\52.13\52.13\52.13\53.18\53.52\54.0]         # 7618,1-n Compensator Thickness Data
            (300a,02e0) CS [ABSENT]                                   # 6,1 Compensator Divergence
            (300a,02e1) CS [SOURCE_SIDE ]                             # 12,1 Compensator Mounting Position
            (300a,02e4) FL 39.2                                       # 4,1 Isocenter to Compensator Tray Distance
            (300a,02e5) FL 2.12                                       # 4,1 Compensator Column Offset
            (300a,02e8) FL 4.76                                       # 4,1 Compensator Milling Tool Diameter
          (fffe,e00d)
    */
      const gdcm::DataSet& ds = reader.GetFile().GetDataSet();
      gdcm::Tag tbeamsq(0x300a,0x03a2);
      if( !ds.FindDataElement( tbeamsq ) )
        {
        return 1;
        }
      const gdcm::DataElement &beamsq = ds.GetDataElement( tbeamsq );
      //std::cout << beamsq << std::endl;
      gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = beamsq.GetValueAsSQ();
      if( !sqi || !sqi->GetNumberOfItems() )
        {
        return 1;
        }

      //for(unsigned int pd = 0; pd < sqi->GetNumberOfItems(); ++pd)
      //  {
        //const gdcm::Item & item = sqi->GetItem(1); // Item start at #1
        const gdcm::Item & item = sqi->GetItem(1); // Item start at #1
        const gdcm::DataSet& nestedds = item.GetNestedDataSet();
        //std::cout << nestedds << std::endl;
        gdcm::Tag tcompensatorsq(0x300a,0x02ea);
        if( !nestedds.FindDataElement( tcompensatorsq ) )
          {
          return 1;
          }
        const gdcm::DataElement &compensatorsq = nestedds.GetDataElement( tcompensatorsq );
        //std::cout << compensatorsq << std::endl;
        gdcm::SmartPointer<gdcm::SequenceOfItems> ssqi = compensatorsq.GetValueAsSQ();
        const gdcm::Item & item2 = ssqi->GetItem(1); // Item start at #1
        const gdcm::DataSet& nestedds2 = item2.GetNestedDataSet();
        //std::cout << nestedds2 << std::endl;
        gdcm::Tag tcompensatorthicknessdata(0x300a,0x00ec);
        if( !nestedds2.FindDataElement( tcompensatorthicknessdata ) )
          {
          return 1;
          }
        const gdcm::DataElement &compensatorthicknessdata = nestedds2.GetDataElement( tcompensatorthicknessdata );
        //  std::cout << compensatorthicknessdata << std::endl;
        gdcm::Attribute<0x300a,0x00ec> at;
        at.SetFromDataElement( compensatorthicknessdata );
        const double* pts = at.GetValues();
        //        (300a,00e7) IS [35]                                       # 2,1 Compensator Rows
        gdcm::Attribute<0x300a,0x00e7> at1;
        const gdcm::DataElement &compensatorrows = nestedds2.GetDataElement( at1.GetTag() );
        at1.SetFromDataElement( compensatorrows );
        std::cout << at1.GetValue() << std::endl;
        //        (300a,00e8) IS [37]                                       # 2,1 Compensator Columns
        gdcm::Attribute<0x300a,0x00e8> at2;
        const gdcm::DataElement &compensatorcols = nestedds2.GetDataElement( at2.GetTag() );
        at2.SetFromDataElement( compensatorcols );
        std::cout << at2.GetValue() << std::endl;

            // (300a,00e9) DS [3.679991\4.249288 ]                       # 18,2 Compensator Pixel Spacing
        gdcm::Attribute<0x300a,0x00e9> at3;
        const gdcm::DataElement &compensatorpixelspacing = nestedds2.GetDataElement( at3.GetTag() );
        at3.SetFromDataElement( compensatorpixelspacing );
        std::cout << at3.GetValue(0) << std::endl;
            // (300a,00ea) DS [-76.00\62.50]                             # 12,2 Compensator Position
        gdcm::Attribute<0x300a,0x00ea> at4;
        const gdcm::DataElement &compensatorposition = nestedds2.GetDataElement( at4.GetTag() );
        at4.SetFromDataElement( compensatorposition );
        std::cout << at4.GetValue(0) << std::endl;

        vtkDoubleArray *d = vtkDoubleArray::New();
        d->SetArray( (double*)pts , at1.GetValue() * at2.GetValue() , 0 );

        vtkImageData *img = vtkImageData::New();
        img->Initialize();
        img->SetDimensions( at2.GetValue(), at1.GetValue(), 1 );
        //imgb->SetExtent(1, xdim, 1, ydim, 1, zdim);
        img->SetScalarTypeToDouble();
        img->SetSpacing( at3.GetValue(1), at3.GetValue(0), 1); // FIXME image is upside down
        img->SetOrigin( at4.GetValue(0), at4.GetValue(1), 1);
        img->SetNumberOfScalarComponents(1);
        img->GetPointData()->SetScalars(d);

        img->Update();
        img->Print(std::cout);

        vtkXMLImageDataWriter *writeb= vtkXMLImageDataWriter::New();
        writeb->SetInput( img );
        writeb->SetFileName( outfilename );
        writeb->Write( );
    /*
        (300a,03a6) SQ                                        # u/l,1 Ion Block Sequence
          (fffe,e000) na (Item with undefined length)
            (300a,00e1) SH [brass ]                           # 6,1 Material ID
            (300a,00f7) FL 95.03                              # 4,1 Isocenter to Block Tray Distance
            (300a,00f8) CS [APERTURE]                         # 8,1 Block Type
            (300a,00fa) CS [ABSENT]                           # 6,1 Block Divergence
            (300a,00fb) CS [SOURCE_SIDE ]                     # 12,1 Block Mounting Position
            (300a,00fc) IS [1 ]                               # 2,1 Block Number
            (300a,0100) DS [50.00 ]                           # 6,1 Block Thickness
            (300a,0104) IS [179 ]                             # 4,1 Block Number of Points
            (300a,0106) DS [1.7\50.0\14.3\50.0\16.7\49.4\18.7\48.2\19.4\47.7\20.1\47.1\21.0\47.0\22.3\47.0\23.7\46.8\25.7\46.2\27.0\45.6\27.2\45.4\28.2\44.6\28.9\44.2\29.7\43.9\31.5\43.5\33.0\42.8\33.7\42.4\35.2\41.3\38.2\40.4\39.6\39.7\40.0\39.5\41.5\37.9\42
    2\37.4\43.0\37.1\44.7\36] # 1934,2-2n Block Data
          (fffe,e00d)
        (fffe,e0dd)

    */
        gdcm::Tag tblocksq(0x300a,0x03a6);
        if( !nestedds.FindDataElement( tblocksq ) )
          {
          return 1;
          }
        const gdcm::DataElement &blocksq = nestedds.GetDataElement( tblocksq );
        //std::cout << blocksq << std::endl;
        gdcm::SmartPointer<gdcm::SequenceOfItems> sssqi = blocksq.GetValueAsSQ();
        const gdcm::Item & item3 = sssqi->GetItem(1); // Item start at #1
        const gdcm::DataSet& nestedds3 = item3.GetNestedDataSet();

        gdcm::Tag tblockdata(0x300a,0x0106);
        if( !nestedds3.FindDataElement( tblockdata ) )
          {
          return 1;
          }
        const gdcm::DataElement &blockdata = nestedds3.GetDataElement( tblockdata );
        //  std::cout << blockdata << std::endl;
        gdcm::Attribute<0x300a,0x0106> at_;
        at_.SetFromDataElement( blockdata );

        vtkDoubleArray *scalars = vtkDoubleArray::New();
        scalars->SetNumberOfComponents(3);

        gdcm::Attribute<0x300a,0x0104> bnpts; //  IS [179 ]                                     # 4,1 Block Number of Points
        if( !nestedds3.FindDataElement( bnpts.GetTag() ) )
          {
          return 1;
          }
        const gdcm::DataElement &blocknpts = nestedds3.GetDataElement( bnpts.GetTag() );
        bnpts.SetFromDataElement(  blocknpts );
      std::cout << bnpts.GetValue() << std::endl;

        vtkPolyData *output = vtkPolyData::New();
        vtkPoints *newPts = vtkPoints::New();
        vtkCellArray *polys = vtkCellArray::New();
        const double *ptr = at_.GetValues();
          //unsigned int npts = bnpts.GetNumberOfValues() / 2;
          unsigned int npts = bnpts.GetValue();
        vtkIdType *ptIds =  new vtkIdType[npts];
        for(unsigned int i = 0; i < npts; ++i)
          {
            float x[3] = {};
            x[0] = ptr[2*i+0];
            x[1] = ptr[2*i+1];
            //x[2] = pts[i+2];
          vtkIdType ptId = newPts->InsertNextPoint( x );
          //std::cout << x[0] << "," << x[1] << "," << x[2] << std::endl;
            ptIds[i ] = ptId;
          }
          vtkIdType cellId = polys->InsertNextCell(npts , ptIds);
            delete[] ptIds;

        output->SetPoints(newPts);
        newPts->Delete();
        output->SetPolys(polys);
        polys->Delete();
        //output->GetCellData()->SetScalars(scalars);
        //scalars->Delete();
        output->Update();
        output->Print( std::cout );

      //  }

       vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();

       vtkImageColorViewer *viewer = vtkImageColorViewer::New();
       viewer->SetInput(img);
       viewer->SetupInteractor(iren);
       viewer->SetSize(600, 600);
       viewer->GetRenderer()->ResetCameraClippingRange();
       viewer->Render();
       viewer->GetRenderer()->ResetCameraClippingRange();

      vtkPolyDataMapper *cubeMapper = vtkPolyDataMapper::New();
      //vtkPolyDataMapper2D* cubeMapper = vtkPolyDataMapper2D::New();
          cubeMapper->SetInput( output );
          cubeMapper->SetScalarRange(0,7);
      vtkActor *cubeActor = vtkActor::New();
      //vtkActor2D* cubeActor = vtkActor2D::New();
          cubeActor->SetMapper(cubeMapper);
      vtkProperty * property = cubeActor->GetProperty();
      property->SetRepresentationToWireframe();

    viewer->GetRenderer()->AddActor( cubeActor );

        vtkXMLPolyDataWriter *writec= vtkXMLPolyDataWriter::New();
        writec->SetInput( output );
        writec->SetFileName( outfilename2 );
        writec->Write( );

       iren->Initialize();
       iren->Start();

      return 0;
    }