# DiscriminateVolume.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 "gdcmScanner.h"
    #include "gdcmTesting.h"
    #include "gdcmIPPSorter.h"
    #include "gdcmDirectionCosines.h"
    #include <cmath>

    /*
     * The following example is a basic sorted which should work in generic cases.
     * It sort files based on:
     * Study Instance UID
     *   Series Instance UID
     *     Frame of Reference UID
     *       Image Orientation (Patient)
     *         Image Position (Patient) (Sorting based on IPP + IOP)
     */

    namespace gdcm {
      const Tag t1(0x0020,0x000d); // Study Instance UID
      const Tag t2(0x0020,0x000e); // Series Instance UID
      const Tag t3(0x0020,0x0052); // Frame of Reference UID
      const Tag t4(0x0020,0x0037); // Image Orientation (Patient)

    class DiscriminateVolume
    {
    private:
      std::vector< Directory::FilenamesType > SortedFiles;
      std::vector< Directory::FilenamesType > UnsortedFiles;

      Directory::FilenamesType GetAllFilenamesFromTagToValue(
        Scanner const & s, Directory::FilenamesType const &filesubset, Tag const &t, const char *valueref)
    {
      Directory::FilenamesType theReturn;
      if( valueref )
        {
        size_t len = strlen( valueref );
        Directory::FilenamesType::const_iterator file = filesubset.begin();
        for(; file != filesubset.end(); ++file)
          {
          const char *filename = file->c_str();
          const char * value = s.GetValue(filename, t);
          if( value && strncmp(value, valueref, len ) == 0 )
            {
            theReturn.push_back( filename );
            }
          }
        }
      return theReturn;
    }

    void ProcessAIOP(Scanner const & s, Directory::FilenamesType const & subset, const char *iopval)
    {
      std::cout << "IOP: " << iopval << std::endl;
      IPPSorter ipp;
      ipp.SetComputeZSpacing( true );
      ipp.SetZSpacingTolerance( 1e-3 ); // ??
      bool b = ipp.Sort( subset );
      if( !b )
        {
        // If you reach here this means you need one more parameter to discriminiat this
        // series. Eg. T1 / T2 intertwinted. Multiple Echo (0018,0081)
        std::cerr << "Failed to sort: " << subset.begin()->c_str() << std::endl;
        for(
          Directory::FilenamesType::const_iterator file = subset.begin();
          file != subset.end(); ++file)
          {
          std::cerr << *file << std::endl;
          }
        UnsortedFiles.push_back( subset );
        return ;
        }
      ipp.Print( std::cout );
      SortedFiles.push_back( ipp.GetFilenames() );
    }

    void ProcessAFrameOfRef(Scanner const & s, Directory::FilenamesType const & subset, const char * frameuid)
    {
      // In this subset of files (belonging to same series), let's find those
      // belonging to the same Frame ref UID:
      Directory::FilenamesType files = GetAllFilenamesFromTagToValue(
        s, subset, t3, frameuid);

      std::set< std::string > iopset;

      for(
        Directory::FilenamesType::const_iterator file = files.begin();
        file != files.end(); ++file)
        {
        //std::cout << *file << std::endl;
        const char * value = s.GetValue(file->c_str(), gdcm::t4 );
        assert( value );
        iopset.insert( value );
        }
      size_t n = iopset.size();
      if ( n == 0 )
        {
        assert( files.empty() );
        return;
        }

      std::cout << "Frame of Ref: " << frameuid << std::endl;
      if ( n == 1 )
        {
        ProcessAIOP(s, files, iopset.begin()->c_str() );
        }
      else
        {
        const char *f = files.begin()->c_str();
        std::cerr << "More than one IOP: " << f << std::endl;
        // Make sure that there is actually 'n' different IOP
        gdcm::DirectionCosines ref;
        gdcm::DirectionCosines dc;
        for(
          std::set< std::string >::const_iterator it = iopset.begin();
          it != iopset.end(); ++it )
          {
          ref.SetFromString( it->c_str() );
          for(
            Directory::FilenamesType::const_iterator file = files.begin();
            file != files.end(); ++file)
            {
            std::string value = s.GetValue(file->c_str(), gdcm::t4 );
            if( value != it->c_str() )
              {
              dc.SetFromString( value.c_str() );
              const double crossdot = ref.CrossDot(dc);
              const double eps = std::fabs( 1. - crossdot );
              if( eps < 1e-6 )
                {
                std::cerr << "Problem with IOP discrimination: " << file->c_str()
                  << " " << it->c_str() << std::endl;
                return;
                }
              }
            }
          }
          // If we reach here this means there is actually 'n' different IOP
        for(
          std::set< std::string >::const_iterator it = iopset.begin();
          it != iopset.end(); ++it )
          {
          const char *iopvalue = it->c_str();
          Directory::FilenamesType iopfiles = GetAllFilenamesFromTagToValue(
            s, files, t4, iopvalue );
          ProcessAIOP(s, iopfiles, iopvalue );
          }
        }
    }

    void ProcessASeries(Scanner const & s, const char * seriesuid)
    {
      std::cout << "Series: " << seriesuid << std::endl;
      // let's find all files belonging to this series:
      Directory::FilenamesType seriesfiles = GetAllFilenamesFromTagToValue(
        s, s.GetFilenames(), t2, seriesuid);

      gdcm::Scanner::ValuesType vt3 = s.GetValues(t3);
      for(
        gdcm::Scanner::ValuesType::const_iterator it = vt3.begin()
        ; it != vt3.end(); ++it )
        {
        ProcessAFrameOfRef(s, seriesfiles, it->c_str());
        }
    }

    void ProcessAStudy(Scanner const & s, const char * studyuid)
    {
      std::cout << "Study: " << studyuid << std::endl;
      gdcm::Scanner::ValuesType vt2 = s.GetValues(t2);
      for(
        gdcm::Scanner::ValuesType::const_iterator it = vt2.begin()
        ; it != vt2.end(); ++it )
        {
        ProcessASeries(s, it->c_str());
        }
    }
    public:

    void Print( std::ostream & os )
    {
      os << "Sorted Files: " << std::endl;
      for(
        std::vector< Directory::FilenamesType >::const_iterator it = SortedFiles.begin();
        it != SortedFiles.end(); ++it )
        {
        os << "Group: " << std::endl;
        for(
          Directory::FilenamesType::const_iterator file = it->begin();
          file != it->end(); ++file)
          {
          os << *file << std::endl;
          }
        }
      os << "Unsorted Files: " << std::endl;
      for(
        std::vector< Directory::FilenamesType >::const_iterator it = UnsortedFiles.begin();
        it != UnsortedFiles.end(); ++it )
        {
        os << "Group: " << std::endl;
        for(
          Directory::FilenamesType::const_iterator file = it->begin();
          file != it->end(); ++file)
          {
          os << *file << std::endl;
          }
        }

    }

      std::vector< Directory::FilenamesType > const & GetSortedFiles() const { return SortedFiles; }
      std::vector< Directory::FilenamesType > const & GetUnsortedFiles() const { return UnsortedFiles; }

    void ProcessIntoVolume( Scanner const & s )
    {
      gdcm::Scanner::ValuesType vt1 = s.GetValues( gdcm::t1 );
      for(
        gdcm::Scanner::ValuesType::const_iterator it = vt1.begin()
        ; it != vt1.end(); ++it )
        {
        ProcessAStudy( s, it->c_str() );
        }

    }

    };

    } // namespace gdcm

    int main(int argc, char *argv[])
    {
      std::string dir1;
      if( argc < 2 )
        {
        const char *extradataroot = NULL;
    #ifdef GDCM_BUILD_TESTING
        extradataroot = gdcm::Testing::GetDataExtraRoot();
    #endif
        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 !

      gdcm::Scanner s;
      s.AddTag( gdcm::t1 );
      s.AddTag( gdcm::t2 );
      s.AddTag( gdcm::t3 );
      s.AddTag( gdcm::t4 );
      bool b = s.Scan( d.GetFilenames() );
      if( !b )
        {
        std::cerr << "Scanner failed" << std::endl;
        return 1;
        }

      gdcm::DiscriminateVolume dv;
      dv.ProcessIntoVolume( s );
      dv.Print( std::cout );

      return 0;
    }