# LargeVRDSExplicit.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 "gdcmReader.h"
    #include "gdcmWriter.h"
    #include "gdcmAttribute.h"
    #include "gdcmFileExplicitFilter.h"
    #include "gdcmSequenceOfItems.h"

    bool interpolate(const double * pts, unsigned int npts, std::vector<double> &out )
    {
      out.clear();
      for(unsigned int i = 0; i < 2*npts; ++i )
        {
        const unsigned int j = i / 2;
        if( i % 2 )
          {
          if( j != npts - 1 )
            {
            assert( j >= 0 && 3*j+5 < 3*npts );
            const double midpointx = (pts[3*j+0] + pts[3*j+3]) / 2;
            const double midpointy = (pts[3*j+1] + pts[3*j+4]) / 2;
            const double midpointz = (pts[3*j+2] + pts[3*j+5]) / 2;
            out.push_back( midpointx );
            out.push_back( midpointy );
            out.push_back( midpointz );
            }
          }
        else
          {
          assert( j >= 0 && j < npts );
          out.push_back( pts[3*j+0] );
          out.push_back( pts[3*j+1] );
          out.push_back( pts[3*j+2] );
          }
        }
      assert( out.size() == 2 * npts * 3 - 3 );
      return true;
    }

    int main(int argc, char *argv[])
    {
      if( argc < 3 )
        {
        std::cerr << argv[0] << " input.dcm output.dcm" << std::endl;
        return 1;
        }
      const char *filename = argv[1];
      const char *outfilename = argv[2];
      gdcm::Reader reader;
      reader.SetFileName( filename );
      if( !reader.Read() )
        {
        return 1;
        }

      gdcm::File &file = reader.GetFile();
      gdcm::DataSet &ds = file.GetDataSet();

      gdcm::FileExplicitFilter fef;
      //fef.SetChangePrivateTags( changeprivatetags );
      fef.SetFile( reader.GetFile() );
      if( !fef.Change() )
        {
        std::cerr << "Failed to change: " << filename << std::endl;
        return 1;
        }

      // (3006,0039) SQ (Sequence with undefined length #=4)     # u/l, 1 ROIContourSequence
      gdcm::Tag tag(0x3006,0x0039);

      const gdcm::DataElement &roicsq = ds.GetDataElement( tag );
      gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = roicsq.GetValueAsSQ();
      //sqi->SetNumberOfItems( 1 );
      const gdcm::Item & item = sqi->GetItem(1); // Item start at #1
      const gdcm::DataSet& nestedds = item.GetNestedDataSet();

      gdcm::Tag tcsq(0x3006,0x0040);
      if( !nestedds.FindDataElement( tcsq ) )
        {
        return 0;
        }
      const gdcm::DataElement& csq = nestedds.GetDataElement( tcsq );
      gdcm::SmartPointer<gdcm::SequenceOfItems> sqi2 = csq.GetValueAsSQ();
      if( !sqi2 || !sqi2->GetNumberOfItems() )
        {
        return 0;
        }
      unsigned int nitems = sqi2->GetNumberOfItems();
      gdcm::Item & item2 = sqi2->GetItem(1); // Item start at #1

      gdcm::DataSet& nestedds2 = item2.GetNestedDataSet();
      //item2.SetVLToUndefined();
      //std::cout << nestedds2 << std::endl;
      // (3006,0050) DS [43.57636\65.52504\-10.0\46.043102\62.564945\-10.0\49.126537\60.714... # 398,48 ContourData
      gdcm::Tag tcontourdata(0x3006,0x0050);
      const gdcm::DataElement & contourdata = nestedds2.GetDataElement( tcontourdata );
      //std::cout << contourdata << std::endl;

      //const gdcm::ByteValue *bv = contourdata.GetByteValue();
      gdcm::Attribute<0x3006,0x0046> ncontourpoints;
      ncontourpoints.Set( nestedds2 );

      gdcm::Attribute<0x3006,0x0050> at;
      at.SetFromDataElement( contourdata );
      const double* pts = at.GetValues();
      unsigned int npts = at.GetNumberOfValues() / 3;

      std::vector<double> out( pts, pts + npts * 3 );
      std::vector<double> out2;

      //const unsigned int niter = 7;
      const unsigned int niter = 8;
      for( unsigned int i = 0; i < niter; ++i)
        {
        bool b = interpolate(&out[0], out.size() / 3, out2);
        //const double *pout = &out[0];
        out = out2;
        out2.clear();
        }
      assert( out.size() % 3 == 0 );

      gdcm::Attribute<0x3006,0x0050> at_interpolate;
      at_interpolate.SetNumberOfValues( out.size() / 3 );
      at_interpolate.SetValues( &out[0], out.size() );

      ncontourpoints.SetValue( at_interpolate.GetNumberOfValues() / 3 );
      nestedds2.Replace( at_interpolate.GetAsDataElement() );
      nestedds2.Replace( ncontourpoints.GetAsDataElement() );

      //assert(0);

      // Let's take item one and subdivide it

      gdcm::TransferSyntax ts = gdcm::TransferSyntax::ImplicitVRLittleEndian;
      ts = gdcm::TransferSyntax::ExplicitVRLittleEndian;

      gdcm::FileMetaInformation &fmi = file.GetHeader();
      const char *tsuid = gdcm::TransferSyntax::GetTSString( ts );
      // const char * is ok since padding is \0 anyway...
      gdcm::DataElement de( gdcm::Tag(0x0002,0x0010) );
      de.SetByteValue( tsuid, strlen(tsuid) );
      de.SetVR( gdcm::Attribute<0x0002, 0x0010>::GetVR() );
      fmi.Replace( de );
      fmi.Remove( gdcm::Tag(0x0002,0x0012) ); // will be regenerated
      fmi.Remove( gdcm::Tag(0x0002,0x0013) ); //  '   '    '
      fmi.SetDataSetTransferSyntax(ts);

      gdcm::Writer w;
      w.SetFile( file );
      w.SetFileName( outfilename );
      if (!w.Write() )
        {
        return 1;
        }

      return 0;
    }