# VolumeSorter.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 "gdcmSorter.h"
    #include "gdcmIPPSorter.h"
    #include "gdcmScanner.h"
    #include "gdcmDataSet.h"
    #include "gdcmAttribute.h"
    #include "gdcmTesting.h"

    bool mysort1(gdcm::DataSet const & ds1, gdcm::DataSet const & ds2 )
    {
      gdcm::Attribute<0x0020,0x000d> at1;
      at1.Set( ds1 );
      gdcm::Attribute<0x0020,0x000d> at2;
      at2.Set( ds2 );
      return at1 < at2;
    }

    bool mysort2(gdcm::DataSet const & ds1, gdcm::DataSet const & ds2 )
    {
      gdcm::Attribute<0x0020,0x000e> at1;
      at1.Set( ds1 );
      gdcm::Attribute<0x0020,0x000e> at2;
      at2.Set( ds2 );
      return at1 < at2;
    }

    bool mysort3(gdcm::DataSet const & ds1, gdcm::DataSet const & ds2 )
    {
      // This is a floating point number is the comparison ok ?
      gdcm::Attribute<0x0020,0x0037> at1;
      at1.Set( ds1 );
      gdcm::Attribute<0x0020,0x0037> at2;
      at2.Set( ds2 );
      return at1 < at2;
    }

    bool mysort4(gdcm::DataSet const & ds1, gdcm::DataSet const & ds2 )
    {
      // Do the IPP sorting here
      gdcm::Attribute<0x0020,0x0032> ipp1;
      gdcm::Attribute<0x0020,0x0037> iop1;
      ipp1.Set( ds1 );
      iop1.Set( ds1 );
      gdcm::Attribute<0x0020,0x0032> ipp2;
      gdcm::Attribute<0x0020,0x0037> iop2;
      ipp2.Set( ds2 );
      iop2.Set( ds2 );
      if( iop1 != iop2 )
        {
        return false;
        }

      // else
      double normal[3];
      normal[0] = iop1[1]*iop1[5] - iop1[2]*iop1[4];
      normal[1] = iop1[2]*iop1[3] - iop1[0]*iop1[5];
      normal[2] = iop1[0]*iop1[4] - iop1[1]*iop1[3];
      double dist1 = 0;
      for (int i = 0; i < 3; ++i) dist1 += normal[i]*ipp1[i];
      double dist2 = 0;
      for (int i = 0; i < 3; ++i) dist2 += normal[i]*ipp2[i];

      std::cout << dist1 << "," << dist2 << std::endl;
      return dist1 < dist2;

    }

    int main(int argc, char *argv[])
    {
      const char *extradataroot = gdcm::Testing::GetDataExtraRoot();
      std::string dir1;
      if( argc < 2 )
        {
        if( !extradataroot )
          {
          return 1;
          }
        dir1 = extradataroot;
        dir1 += "/gdcmSampleData/ForSeriesTesting/VariousIncidences/ST1";
        }
      else
        {
        dir1 = argv[1];
        }

      gdcm::Directory d;
      d.Load( dir1.c_str(), true ); // recursive !
      const gdcm::Directory::FilenamesType &l1 = d.GetFilenames();
      const unsigned int nfiles = l1.size();
      std::cout << nfiles << std::endl;

      //if( nfiles != 280 )
      //  {
      //  return 1;
      //  }

      //d.Print( std::cout );

      gdcm::Scanner s;
      const gdcm::Tag t1(0x0020,0x000d); // Study Instance UID
      const gdcm::Tag t2(0x0020,0x000e); // Series Instance UID
      //const gdcm::Tag t3(0x0010,0x0010); // Patient's Name
      s.AddTag( t1 );
      s.AddTag( t2 );
      //s.AddTag( t3 );
      //s.AddTag( t4 );
      //s.AddTag( t5 );
      //s.AddTag( t6 );
      bool b = s.Scan( d.GetFilenames() );
      if( !b )
        {
        std::cerr << "Scanner failed" << std::endl;
        return 1;
        }

      //s.Print( std::cout );

      // Only get the DICOM files:
      gdcm::Directory::FilenamesType l2 = s.GetKeys();
      const unsigned int nfiles2 = l2.size();
      std::cout << nfiles2 << std::endl;

      if ( nfiles2 > nfiles )
        {
        return 1;
        }

      gdcm::Sorter sorter;
      sorter.SetSortFunction( mysort1 );
      sorter.StableSort( l2 );

      sorter.SetSortFunction( mysort2 );
      sorter.StableSort( sorter.GetFilenames() );

      sorter.SetSortFunction( mysort3 );
      sorter.StableSort( sorter.GetFilenames() );

      sorter.SetSortFunction( mysort4 );
      sorter.StableSort( sorter.GetFilenames() );

      //sorter.Print( std::cout );

      // Let's try to check our result:
      // assume that IPP is precise enough so that we can test floating point equality:
      unsigned int nvalues = 0;
    {
      gdcm::Scanner s;
      s.AddTag( gdcm::Tag(0x20,0x32) ); // Image Position (Patient)
      //s.AddTag( gdcm::Tag(0x20,0x37) ); // Image Orientation (Patient)
      s.Scan( d.GetFilenames() );

      //s.Print( std::cout );

      const gdcm::Scanner::ValuesType &values = s.GetValues();
      nvalues = values.size();
      std::cout << "There are " << nvalues << " different type of values" << std::endl;
      assert( nfiles2 % nvalues == 0 );
      std::cout << "Series is composed of " << (nfiles/nvalues) << " different 3D volumes" << std::endl;
    }

      gdcm::Directory::FilenamesType sorted_files = sorter.GetFilenames();

      // Which means we can take nvalues files at a time and execute gdcm::IPPSorter on it:
      gdcm::IPPSorter ippsorter;
      gdcm::Directory::FilenamesType sub( sorted_files.begin(), sorted_files.begin() + nvalues);
      std::cout << sub.size() << std::endl;
      std::cout << sub[0] << std::endl;
      std::cout << sub[nvalues-1] << std::endl;
      ippsorter.SetComputeZSpacing( false );
      if( !ippsorter.Sort( sub ) )
        {
        std::cerr << "Could not sort" << std::endl;
        return 1;
        }

      std::cout << "IPPSorter:" << std::endl;
      ippsorter.Print( std::cout );

      return 0;
    }