# csa2img.cxx

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

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

      Copyright (c) 2006-2009 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.

    =========================================================================*/
    /*
     * I do not know what the format is, just guessing from info found on the net:
     *
     * http://atonal.ucdavis.edu/matlab/fmri/spm5/spm_dicom_convert.m
     *
     * This example is an attempt at understanding the format used by SIEMENS
     * their "SIEMENS CSA NON-IMAGE" DICOM file (1.3.12.2.1107.5.9.1)
     *
     * Everything done in this code is for the sole purpose of writing interoperable
     * software under Sect. 1201 (f) Reverse Engineering exception of the DMCA.
     * If you believe anything in this code violates any law or any of your rights,
     * please contact us (gdcm-developers@lists.sourceforge.net) so that we can
     * find a solution.
     *
     */
    #include "gdcmReader.h"
    #include "gdcmImageReader.h"
    #include "gdcmImageWriter.h"
    #include "gdcmCSAHeader.h"
    #include "gdcmAttribute.h"
    #include "gdcmPrivateTag.h"

    #include <math.h>

    /*
     *  gdcmSampleData/images_of_interest/MR-sonata-3D-as-Tile.dcm
     */
    bool reorganize_mosaic(const unsigned short *input, const unsigned int *inputdims, unsigned int square,
      const unsigned int *outputdims, unsigned short *output )
    {
      for(unsigned x = 0; x < outputdims[0]; ++x)
        {
        for(unsigned y = 0; y < outputdims[1]; ++y)
          {
          for(unsigned z = 0; z < outputdims[2]; ++z)
            {
            output[ x + y*outputdims[0] + z*outputdims[0]*outputdims[1] ] =
              input[ (x + z * outputdims[0]) + (y + (z/square)*outputdims[0])*inputdims[0] ];
            }
          }
        }
      return true;
    }

    int main(int argc, char *argv [])
    {
      if( argc < 2 ) return 1;
      // gdcmDataExtra/gdcmNonImageData/exCSA_Non-Image_Storage.dcm
      // PHANTOM.MR.CARDIO_COEUR_S_QUENCE_DE_REP_RAGE.9.257.2008.03.20.14.53.25.578125.43151705.IMA
      const char *filename = argv[1];
      gdcm::Reader reader;
      reader.SetFileName( filename );
      if( !reader.Read() )
        {
        std::cerr << "Failed to read: " << filename << std::endl;
        return 1;
        }

      gdcm::CSAHeader csa;
      const gdcm::DataSet& ds = reader.GetFile().GetDataSet();

      const gdcm::PrivateTag &t1 = csa.GetCSAImageHeaderInfoTag();
      //std::cout << t1 << std::endl;
      const gdcm::PrivateTag &t2 = csa.GetCSASeriesHeaderInfoTag();

      if( ds.FindDataElement( t1 ) )
        {
        csa.LoadFromDataElement( ds.GetDataElement( t1 ) );
        csa.Print( std::cout );
        }
      //const gdcm::Image &inputimage = reader.GetImage();
    #if 1
      int dims[2];
      const gdcm::CSAElement &csael = csa.GetCSAElementByName( "Columns" );
      std::cout << csael << std::endl;
      const gdcm::ByteValue *bv = csael.GetByteValue();
      gdcm::Element<gdcm::VR::IS, gdcm::VM::VM1> el;
      el.Set( csael.GetValue() );
      dims[0] = el.GetValue();
      std::cout << "Columns:" << el.GetValue() << std::endl;

      const gdcm::CSAElement &csael2 = csa.GetCSAElementByName( "Rows" );
      std::cout << csael2 << std::endl;
      gdcm::Element<gdcm::VR::IS, gdcm::VM::VM1> el2;
      el2.Set( csael2.GetValue() );
      dims[1] = el2.GetValue();
      std::cout << "Rows:" << el2.GetValue() << std::endl;

      const gdcm::CSAElement &csael3 = csa.GetCSAElementByName( "PixelSpacing" );
      std::cout << csael3 << std::endl;
      gdcm::Element<gdcm::VR::DS, gdcm::VM::VM2> el3;
      el3.Set( csael3.GetValue() );
      double spacing[2];
      spacing[0] = el3.GetValue(0);
      spacing[1] = el3.GetValue(1);
      std::cout << "PixelSpacing:" << el3.GetValue() << "," << el3.GetValue(1) << std::endl;

      gdcm::ImageWriter writer;

      gdcm::Image &image = writer.GetImage();
      image.SetNumberOfDimensions( 2 ); // good default
      image.SetDimension(0, dims[0] );
      image.SetDimension(1, dims[1] );
      image.SetSpacing(0, spacing[0] );
      image.SetSpacing(1, spacing[1] );
      gdcm::PixelFormat pixeltype = gdcm::PixelFormat::INT16; // bytepix = spm_type('int16','bits')/8;

      unsigned long l = image.GetBufferLength();
      const int p =  l / (dims[0] * dims[1]);

      //image.SetNumberOfDimensions( 3 );
      //image.SetDimension(2, p / pixeltype.GetPixelSize() );

      gdcm::PhotometricInterpretation pi;
      pi = gdcm::PhotometricInterpretation::MONOCHROME2;
      //pixeltype.SetSamplesPerPixel(  );
      image.SetPhotometricInterpretation( pi );
      image.SetPixelFormat( pixeltype );
      //image.SetIntercept( inputimage.GetIntercept() );
      //image.SetSlope( inputimage.GetSlope() );

      //gdcm::DataElement pixeldata( gdcm::Tag(0x7fe1,0x1010) );
      //pixeldata.SetByteValue( &outbuf[0], outbuf.size() );
      gdcm::PrivateTag csanonimaget(0x7fe1,0x10,"SIEMENS CSA NON-IMAGE");
      const gdcm::DataElement &pixeldata = ds.GetDataElement( csanonimaget );
      image.SetDataElement( pixeldata );

      std::string outfilename = "outcsa.dcm";
      //writer.SetFile( reader.GetFile() );
      writer.SetFileName( outfilename.c_str() );
      if( !writer.Write() )
        {
        std::cerr << "could not write: " << outfilename << std::endl;
        return 1;
        }

    #else
      // The following is deprecated.
      // See gdcmconv --mosaic instead

      //int dims[2] = { 448, 448 };
      //dims[0] /= 7;
      //dims[1] /= 7;
        unsigned int dims[3] = {};
      //double spacing[2] = {1,1};
        {
        const gdcm::DataElement& de = ds.GetDataElement( gdcm::Tag(0x0028, 0x0011) );
        gdcm::Attribute<0x0028,0x0011> at;
        at.SetFromDataElement( de );
        //pixeldata.SetDimension(0, at.GetValue() );
        dims[0] = at.GetValue();
        }

        {
        const gdcm::DataElement& de = ds.GetDataElement( gdcm::Tag(0x0028, 0x0010) );
        gdcm::Attribute<0x0028,0x0010> at;
        at.SetFromDataElement( de );
        //pixeldata.SetDimension(1, at.GetValue() );
        dims[1] = at.GetValue();
        }

      // SliceThickness ??
      const gdcm::CSAElement &csael4 = csa.GetCSAElementByName( "NumberOfImagesInMosaic" );
      std::cout << csael4 << std::endl;
      gdcm::Element<gdcm::VR::IS, gdcm::VM::VM1> el4;
      el4.Set( csael4.GetValue() );
      int numberOfImagesInMosaic = el4.GetValue();

      unsigned int div = (unsigned int)ceil(sqrt( (double)numberOfImagesInMosaic ) );
      dims[0] /= div;
      dims[1] /= div;
      dims[2] = numberOfImagesInMosaic;

      std::cout << "NumberOfImagesInMosaic:" << numberOfImagesInMosaic << std::endl;
      const double *spacing = inputimage.GetSpacing();

      unsigned long l = inputimage.GetBufferLength();
      std::vector<char> buf;
      buf.resize(l);
      inputimage.GetBuffer( &buf[0] );
      gdcm::DataElement pixeldata( gdcm::Tag(0x7fe1,0x1010) );

      std::vector<char> outbuf;
      outbuf.resize(l);
      //outbuf = buf;

      bool b = reorganize_mosaic((unsigned short*)&buf[0], inputimage.GetDimensions(), div, dims, (unsigned short*)&outbuf[0] );

      pixeldata.SetByteValue( &outbuf[0], outbuf.size() );
      //const gdcm::DataElement & pixeldata = ds.GetDataElement( gdcm::Tag(0x7fe1,0x1010) );
      //const gdcm::DataElement & pixeldata = ds.GetDataElement( gdcm::Tag(0x7fe0,0x0010) );
      //const gdcm::VL &l = pixeldata.GetVL();
      const int p =  l / (dims[0] * dims[1]);
      std::cout << "VL:" << l << std::endl;
      std::cout << "pixel:" << p << std::endl;

      gdcm::ImageWriter writer;

      gdcm::Image &image = writer.GetImage();
      image.SetNumberOfDimensions( 2 ); // good default
      image.SetDimension(0, dims[0] );
      image.SetDimension(1, dims[1] );
      image.SetSpacing(0, spacing[0] );
      image.SetSpacing(1, spacing[1] );
      gdcm::PixelFormat pixeltype = inputimage.GetPixelFormat(); //gdcm::PixelFormat::INT16; // bytepix = spm_type('int16','bits')/8;

      image.SetNumberOfDimensions( 3 );
      //image.SetDimension(2, p / pixeltype.GetPixelSize() );
      image.SetDimension(2, numberOfImagesInMosaic );

      gdcm::PhotometricInterpretation pi;
      pi = gdcm::PhotometricInterpretation::MONOCHROME2;
      //pixeltype.SetSamplesPerPixel(  );
      image.SetPhotometricInterpretation( pi );
      image.SetPixelFormat( pixeltype );
      image.SetIntercept( inputimage.GetIntercept() );
      image.SetSlope( inputimage.GetSlope() );

      image.SetDataElement( pixeldata );

      std::string outfilename = "outcsa.dcm";
      //writer.SetFile( reader.GetFile() );
      writer.SetFileName( outfilename.c_str() );
      if( !writer.Write() )
        {
        std::cerr << "could not write: " << outfilename << std::endl;
        return 1;
        }
    #endif

      return 0;
    }