# threadgdcm.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 "gdcmImageReader.h"
    #include "gdcmDirectory.h"
    #include "gdcmSystem.h"

    #include "vtkImageData.h"
    #include "vtkStructuredPointsWriter.h"

    #include <pthread.h>

    struct threadparams
    {
      const char **filenames;
      unsigned int nfiles;
      char *scalarpointer;
    // TODO I should also pass in the dim of the reference image just in case
    };

    void *ReadFilesThread(void *voidparams)
    {
      const threadparams *params = static_cast<const threadparams *> (voidparams);

      const unsigned int nfiles = params->nfiles;
      for(unsigned int file = 0; file < nfiles; ++file)
        {
        /*
        // TODO: update progress
        pthread_mutex_lock(&params->lock);
        //section critique
        ReadingProgress+=params->stepProgress;
        pthread_mutex_unlock(&params->lock);
        */
        const char *filename = params->filenames[file];
        //std::cerr << filename << std::endl;

        gdcm::ImageReader reader;
        reader.SetFileName( filename );
        try
          {
          if( !reader.Read() )
            {
            std::cerr << "Failed to read: " << filename << std::endl;
            break;
            }
          }
        catch( ... )
          {
          std::cerr << "Failed to read: " << filename << std::endl;
          break;
          }

        const gdcm::Image &image = reader.GetImage();
        unsigned long len = image.GetBufferLength();
        char * pointer = params->scalarpointer;
    #if 0
        char *tempimage = new char[len];
        image.GetBuffer(tempimage);

        memcpy(pointer + file*len, tempimage, len);
        delete[] tempimage;
    #else
        char *tempimage = pointer + file * len;
        image.GetBuffer(tempimage);
    #endif
        }

      return voidparams;
    }

    void ShowFilenames(const threadparams &params)
    {
      std::cout << "start" << std::endl;
      for(unsigned int i = 0; i < params.nfiles; ++i)
        {
        const char *filename = params.filenames[i];
        std::cout << filename << std::endl;
        }
      std::cout << "end" << std::endl;
    }

    void ReadFiles(unsigned int nfiles, const char *filenames[])
    {
      // \precondition: nfiles > 0
      assert( nfiles > 0 );
      const char *reference= filenames[0]; // take the first image as reference

      gdcm::ImageReader reader;
      reader.SetFileName( reference );
      if( !reader.Read() )
        {
        // That would be very bad...
        assert(0);
        }

      const gdcm::Image &image = reader.GetImage();
      gdcm::PixelFormat pixeltype = image.GetPixelFormat();
      unsigned long len = image.GetBufferLength();
      const unsigned int *dims = image.GetDimensions();
      unsigned short pixelsize = pixeltype.GetPixelSize();
      (void)pixeltype;
      assert( image.GetNumberOfDimensions() == 2 );

      vtkImageData *output = vtkImageData::New();
      output->SetDimensions(dims[0], dims[1], nfiles);

      switch( pixeltype )
        {
      case gdcm::PixelFormat::INT8:
    #if (VTK_MAJOR_VERSION >= 5) || ( VTK_MAJOR_VERSION == 4 && VTK_MINOR_VERSION > 5 )
        output->SetScalarType ( VTK_SIGNED_CHAR );
    #else
        output->SetScalarType ( VTK_CHAR );
    #endif
        break;
      case gdcm::PixelFormat::UINT8:
        output->SetScalarType ( VTK_UNSIGNED_CHAR );
        break;
      case gdcm::PixelFormat::INT16:
        output->SetScalarType ( VTK_SHORT );
        break;
      case gdcm::PixelFormat::UINT16:
        output->SetScalarType ( VTK_UNSIGNED_SHORT );
        break;
      case gdcm::PixelFormat::INT32:
        output->SetScalarType ( VTK_INT );
        break;
      case gdcm::PixelFormat::UINT32:
        output->SetScalarType ( VTK_UNSIGNED_INT );
        break;
      default:
        assert(0);
        }

      output->SetNumberOfScalarComponents ( pixeltype.GetSamplesPerPixel() );

      output->AllocateScalars();
      char * scalarpointer = static_cast<char*>(output->GetScalarPointer());

      const unsigned int nthreads = 4;
      threadparams params[nthreads];

      //pthread_mutex_t lock;
      //pthread_mutex_init(&lock, NULL);

      pthread_t *pthread = new pthread_t[nthreads];

      // There is nfiles, and nThreads
      assert( nfiles > nthreads );
      const unsigned int partition = nfiles / nthreads;
      for (unsigned int thread=0; thread < nthreads; ++thread)
        {
        params[thread].filenames = filenames + thread * partition;
        params[thread].nfiles = partition;
        if( thread == nthreads - 1 )
          {
          // There is slightly more files to process in this thread:
          params[thread].nfiles += nfiles % nthreads;
          }
        assert( thread * partition < nfiles );
        params[thread].scalarpointer = scalarpointer + thread * partition * len;
        //assert( params[thread].scalarpointer < scalarpointer + 2 * dims[0] * dims[1] * dims[2] );
        // start thread:
        int res = pthread_create( &pthread[thread], NULL, ReadFilesThread, &params[thread]);
        if( res )
          {
          std::cerr << "Unable to start a new thread, pthread returned: " << res << std::endl;
          assert(0);
          }
        //ShowFilenames(params[thread]);
        }
    // DEBUG
      unsigned int total = 0;
      for (unsigned int thread=0; thread < nthreads; ++thread)
        {
        total += params[thread].nfiles;
        }
      assert( total == nfiles );
    // END DEBUG

      for (unsigned int thread=0;thread<nthreads;thread++)
        {
        pthread_join( pthread[thread], NULL);
        }
      delete[] pthread;

      //pthread_mutex_destroy(&lock);

      // For some reason writing down the file is painfully slow...
      vtkStructuredPointsWriter *writer = vtkStructuredPointsWriter::New();
      writer->SetInput( output );
      writer->SetFileName( "/tmp/threadgdcm.vtk" );
      writer->SetFileTypeToBinary();
      //writer->Write();
      writer->Delete();

      //output->Print( std::cout );
      output->Delete();
    }

    int main(int argc, char *argv[])
    {
      if( argc < 2 )
        {
        std::cerr << argv[0] << " [directory|list of filenames]\n";
        return 1;
        }

      // Check if user pass in a single directory
      if( argc == 2 && gdcm::System::FileIsDirectory( argv[1] ) )
        {
        gdcm::Directory d;
        d.Load( argv[1] );
        gdcm::Directory::FilenamesType l = d.GetFilenames();
        const unsigned int nfiles = l.size();
        const char **filenames = new const char* [ nfiles ];
        for(unsigned int i = 0; i < nfiles; ++i)
          {
          filenames[i] = l[i].c_str();
          }
        ReadFiles(nfiles, filenames);
        delete[] filenames;
        }
      else
        {
        // Simply copy all filenames into the vector:
        const char **filenames = const_cast<const char**>(argv+1);
        const unsigned int nfiles = argc - 1;
        ReadFiles(nfiles, filenames);
        }

      return 0;
    }