/*=========================================================================
                                                                                
  Program:   gdcm
  Module:    $RCSfile: Dense2007ToDicom.cxx,v $
  Language:  C++
  Date:      $Date: 2007/09/18 10:54:23 $
  Version:   $Revision: 1.2 $
                                                                                
  Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  l'Image). All rights reserved. See Doc/License.txt or
  http://www.creatis.insa-lyon.fr/Public/Gdcm/License.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 notices for more information.
                                                                                
=========================================================================*/

#include <fstream>
#include <iostream>
//#include <values.h>

#include "gdcmFile.h"
#include "gdcmFileHelper.h"
#include "gdcmDebug.h"
#include "gdcmDirList.h"

#include "gdcmArgMgr.h"

/**
  * \brief   
  *  Converts the "Dense" ".txt" (2007 version)  files into 16 bits Dicom Files,
  * Hope they don't change soon!
  */  


void LoadPeakStrain(std::ifstream &from, std::string imageName);
void LoadStrain(std::ifstream &from, std::string imageName);
void MakeDicomImage(float *tabVal, float *X, float *Y, float *Z, int NP, std::string dcmImageName);


int main(int argc, char *argv[])
{
   START_USAGE(usage)
   " \n Dense2007ToDicom :\n                                                  ",
   "        Converts the '.txt' files into 16 bits Dicom Files,               ",
   " usage:                                                                   ",
   " Dense2007ToDicom strain=...strain.txt  peak_strain=...peak_strain.txt    ",
   "                 [verbose] [debug]                                        ",
   "                                                                          ",
   " verbose  : user wants to run the program in 'verbose mode'               ",
   " debug    : *developer*  wants to run the program in 'debug mode'         ",
   FINISH_USAGE

   // ----- Initialize Arguments Manager ------
      
   GDCM_NAME_SPACE::ArgMgr *am = new GDCM_NAME_SPACE::ArgMgr(argc, argv);
  
   if (argc == 1 || am->ArgMgrDefined("usage")) 
   {
      am->ArgMgrUsage(usage); // Display 'usage'
      delete am;
      return 0;
   }

   const char *strain      = am->ArgMgrWantString("strain",usage);
   const char *peak_strain = am->ArgMgrWantString("peak_strain",usage);
   
   if (am->ArgMgrDefined("debug"))
      GDCM_NAME_SPACE::Debug::DebugOn();

   int verbose  = am->ArgMgrDefined("verbose");      

   // if unused Param we give up
   if ( am->ArgMgrPrintUnusedLabels() )
   { 
      am->ArgMgrUsage(usage);
      delete am;
      return 0;
   }
   delete am;  // we don't need Argument Manager any longer

   // ----- Begin Processing -----


   std::ifstream fromPeakStrain( peak_strain );             
   if ( !fromPeakStrain )
   {
      std::cout << "Can't open file" << peak_strain << std::endl;
      exit(0);
   }

   std::ifstream fromStrain( strain );      
   if ( !fromStrain )
   {
      std::cout << "Can't open file" << strain << std::endl;
      exit(0);
   }  
       
   std::cout << "Success in open file" << peak_strain << std::endl;
   LoadPeakStrain(fromPeakStrain, peak_strain);
   fromPeakStrain.close();  

   std::cout << "Success in open file" << strain << std::endl;
   LoadStrain(fromStrain, strain);
   fromStrain.close();      
   return 1;            
}

// =====================================================================================================================

void LoadPeakStrain(std::ifstream &from, std::string textFileName)
{
// in sax_base_slice0_peak_strain.txt :

/*
Number of material points (NP) = 181
Origin of (readout, phase enc, slice sel) coordinates in 3D =  87.3243 3.19392
88.2381
Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046
Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
followed by their peak Ecc strain, an array of NP elements,
followed by their peak Err strain, an array of NP elements,
followed by their peak E11 strain, an array of NP elements,
followed by their Peak E22 strain, an array of NP elements,
      42.0000      10.0000     0.000000
      ...
    -0.154905   -0.0840482    -0.157350    -0.221403    -0.168118    -0.131331
    -0.153781    -0.148481    -0.166602    -0.232858    -0.222650    -0.213712
    ...
*/  

   if (!from)
      return;

   std::string str1;
   int NP;

   //Number of material points (NP) = 181   
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> NP;
    
    std::cout << "NP : " << NP << std::endl; 

   //Origin of (readout, phase enc, slice sel) coordinates in 3D =  87.3243 3.19392 88.2381
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     
     float readout,  phase_enc, slice_sel;
     from >> readout;
     from >> phase_enc;
     from >> slice_sel;
     std::cout << " readout " << readout << " phase_enc " << phase_enc << " slice_sel " << slice_sel << std::endl;
     
    // Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
    
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    
    float readoutX, readoutY, readoutZ;
    from >> readoutX;
    from >> readoutY;       
    from >> readoutZ;
    std::cout << " readoutX " << readoutX <<  " readoutY " << readoutY <<  " readoutZ " << readoutZ << std::endl;
     
// Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046

     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     
    float phase_encX, phase_encY, phase_encZ;
    from >> phase_encX;
    from >> phase_encY;       
    from >> phase_encZ;
    std::cout << " phase_encX " << phase_encX <<  " phase_encY " << phase_encY <<  " phase_encZ " << phase_encZ << std::endl; 

// Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     
    float slice_selX, slice_selY, slice_selZ;
    from >> slice_selX;
    from >> slice_selY;       
    from >> slice_selZ;
    std::cout << " slice_selX " << slice_selX <<  " slice_selY " << slice_selY <<  " slice_selZ " << slice_selZ << std::endl; 



// Skip 5 lines :
/*
The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
followed by their peak Ecc strain, an array of NP elements,
followed by their peak Err strain, an array of NP elements,
followed by their peak E11 strain, an array of NP elements,
followed by their Peak E22 strain, an array of NP elements,
*/

std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------stop skipping ---------------- " << std::endl;
  
   float *X = new float[NP];
   float *Y = new float[NP];   
   float *Z = new float[NP];
      
   char c;
   int i;   
   for (i=0; i<NP; i++) {   
      from >> X[i];
      for (;;) {
        if (!from.get(c))
          break;
        if (!isspace(c)) {
          from.putback(c);
          break;
        }
     }  
      from >> Y[i];    
      for (;;) {
        if (!from.get(c))
          break;
        if (!isspace(c)) {
          from.putback(c);
          break;
        }
     }        
      from >> Z[i];              
   }

   std::cout << "--------------- Ecc_strain ------------------" << std::endl;
   float *ecc_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> ecc_strain[i]; 
       std::cout <<  ecc_strain[i] <<  std::endl;
   }

   std::cout << "--------------- Err_strain ------------------" << std::endl;
   float *err_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> err_strain[i]; 
       std::cout <<  err_strain[i] <<  std::endl;
   }
   
   std::cout << "--------------- E11_strain ------------------" << std::endl;
   float *e11_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> e11_strain[i]; 
       std::cout <<  e11_strain[i] <<  std::endl;
   }  
   
   std::cout << "--------------- E22_strain ------------------" << std::endl;
   float *e22_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> e22_strain[i]; 
       std::cout <<  e22_strain[i] <<  std::endl;
   }    
 
   std::string dcmImageName;    

//followed by their peak Ecc strain, an array of NP elements,
   dcmImageName = textFileName + "_peak_Ecc_strain.dcm";
   MakeDicomImage(ecc_strain, X, Y, Z, NP, dcmImageName);

//followed by their peak Err strain, an array of NP elements,
   dcmImageName = textFileName + "_peak_Err_strain.dcm";
   MakeDicomImage(err_strain, X, Y, Z, NP, dcmImageName);
   
//followed by their peak E11 strain, an array of NP elements,
   dcmImageName = textFileName + "_peak_E11_strain.dcm";
   MakeDicomImage(e11_strain, X, Y, Z, NP, dcmImageName);
   
//followed by their Peak E22 strain, an array of NP elements,
   dcmImageName = textFileName + "_peak_E22_strain.dcm";
   MakeDicomImage(e22_strain, X, Y, Z, NP, dcmImageName);         
}

// =====================================================================================================================

void LoadStrain(std::ifstream &from, std::string textFileName)
{

// in sax_base_slice0_strain.txt :
/*
Number of cine frames = 18
Temporal resolution = 32.0000 ms
First frame starts at 48.0000 ms
Number of material points (NP) = 181
Origin of (readout, phase enc, slice sel) coordinates in 3D =  87.324341 3.193918 88.238113 
Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046
Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
followed by their Ecc strain, an array of dimensions(NP, number of cine frames),
followed by their Err strain, an array of dimensions(NP, number of cine frames),
followed by their E11 strain, an array of dimensions(NP, number of cine frames),
followed by their E22 strain, an array of dimensions(NP, number of cine frames),
Note that RV Err, E11 and E22 strains are not calculated due to the small thickness.
      42.0000      10.0000     0.000000
      44.0000      10.0000     0.000000
      ...
   -0.0622793   -0.0840482    -0.157350    -0.196722    -0.105844    -0.131331
    -0.153781  -0.00940573   -0.0542236    -0.100403   -0.0369671   -0.0696840      
*/ 

   if (!from)
      return;

   std::string str1;
   int NP;    // Number of Pints
   int NCF;   // Number of cine frames
   float TR;  // Temporal resolution
   float FFS; // First frame starts
   
   // Number of cine frames = 18
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> NCF;

   // Temporal resolution = 32.0000 ms
    from >> str1;
    from >> str1;
    from >> str1;
    from >> TR;
    from >> str1;

   // First frame starts at 48.0000 ms
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> FFS;
    from >> str1;        

   //Number of material points (NP) = 181   
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> NP;
    
    std::cout << "NP : " << NP << std::endl; 

   //Origin of (readout, phase enc, slice sel) coordinates in 3D = 87.324341 3.193918 88.238113 
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     
     float readout,  phase_enc, slice_sel;
     from >> readout;
     from >> phase_enc;
     from >> slice_sel;
     std::cout << " readout " << readout << " phase_enc " << phase_enc << " slice_sel " << slice_sel << std::endl;
     
    // Readout direction in 3D =  -0.162314 -0.0771294 -0.983720
    
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    from >> str1;
    
    float readoutX, readoutY, readoutZ;
    from >> readoutX;
    from >> readoutY;       
    from >> readoutZ;
    std::cout << " readoutX " << readoutX <<  " readoutY " << readoutY <<  " readoutZ " << readoutZ << std::endl;
     
// Phase Enc. direction in 3D =  -0.540606 -0.827052 0.154046

     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     
    float phase_encX, phase_encY, phase_encZ;
    from >> phase_encX;
    from >> phase_encY;       
    from >> phase_encZ;
    std::cout << " phase_encX " << phase_encX <<  " phase_encY " << phase_encY <<  " phase_encZ " << phase_encZ << std::endl; 

// Slice select direction in 3D =  0.825469 -0.556809 -0.0925458
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     from >> str1;
     
    float slice_selX, slice_selY, slice_selZ;
    from >> slice_selX;
    from >> slice_selY;       
    from >> slice_selZ;
    std::cout << " slice_selX " << slice_selX <<  " slice_selY " << slice_selY <<  " slice_selZ " << slice_selZ << std::endl; 



// Skip 6 lines :
/*
The following are the (readout, phase enc, slice sel) coordinates (mm) of the grid points for which strains are calculated,
followed by their Ecc strain, an array of dimensions(NP, number of cine frames),
followed by their Err strain, an array of dimensions(NP, number of cine frames),
followed by their E11 strain, an array of dimensions(NP, number of cine frames),
followed by their E22 strain, an array of dimensions(NP, number of cine frames),
Note that RV Err, E11 and E22 strains are not calculated due to the small thickness.
*/
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------start skipping 1 line---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------stop skipping ---------------- " << std::endl;
   std::getline(from, str1);
   std::cout << "[" << str1 << "]" << std::endl;
std::cout << "------------stop skipping ---------------- " << std::endl;
  
   float *X = new float[NP];
   float *Y = new float[NP];   
   float *Z = new float[NP];
      
   char c;
   int i;   
   for (i=0; i<NP; i++) {
   
      from >> X[i];
      for (;;) {
        if (!from.get(c))
          break;
        if (!isspace(c)) {
          from.putback(c);
          break;
        }
     }  
      from >> Y[i];  
      for (;;) {
        if (!from.get(c))
          break;
        if (!isspace(c)) {
          from.putback(c);
          break;
        }
     }        
      from >> Z[i];              
   }

   std::cout << "--------------- Ecc_strain ------------------" << std::endl;
   float *ecc_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> ecc_strain[i]; 
       std::cout <<  ecc_strain[i] <<  std::endl;
   }

   std::cout << "--------------- Err_strain ------------------" << std::endl;
   float *err_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> err_strain[i]; 
       std::cout <<  err_strain[i] <<  std::endl;
   }
   
   std::cout << "--------------- E11_strain ------------------" << std::endl;
   float *e11_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> e11_strain[i]; 
       std::cout <<  e11_strain[i] <<  std::endl;
   }  
   
   std::cout << "--------------- E22_strain ------------------" << std::endl;
   float *e22_strain = new float[NP];
   for (i=0; i<NP; i++) {
       from >> e22_strain[i]; 
       std::cout <<  e22_strain[i] <<  std::endl;
   }    

   std::string dcmImageName;    

//followed by their Ecc strain, an array of NP elements,
   dcmImageName = textFileName + "_Ecc_strain.dcm";
   MakeDicomImage(ecc_strain, X, Y, Z, NP, dcmImageName);

//followed by their  Err strain, an array of NP elements,
   dcmImageName = textFileName + "_Err_strain.dcm";
   MakeDicomImage(err_strain, X, Y, Z, NP, dcmImageName);
   
//followed by their E11 strain, an array of NP elements,
   dcmImageName = textFileName + "_E11_strain.dcm";
   MakeDicomImage(e11_strain, X, Y, Z, NP, dcmImageName);
   
//followed by their E22 strain, an array of NP elements,
   dcmImageName = textFileName + "_E22_strain.dcm";
   MakeDicomImage(e22_strain, X, Y, Z, NP, dcmImageName);   
      
}


// =====================================================================================================================
    

void MakeDicomImage(float *tabVal, float *X, float *Y, float *Z, int NP, std::string dcmImageName)
{
   float minX = 999., minY = 999., minZ = 999.;
   float maxX = 0., maxY = 0., maxZ = 0.;
   int i;
   for (i=0; i<NP; i++) {
      // std::cout << X[i] << " " << Y[i] << " " << Z[i] <<  std::endl;
      if(maxX < X[i])
         maxX = X[i];
      if(maxY < Y[i])
         maxY = Y[i];
      if(maxZ < Z[i])
         maxZ = Z[i];
 
      if(minX > X[i])
         minX = X[i];
      if(minY > Y[i])
         minY = Y[i];
      if(minZ > Z[i])
         minZ = Z[i];
 
   }   
      std::cout << "Min X,Y,Z " << minX << " " << minY << " " << minZ <<  std::endl;
      std::cout << "Max X,Y,Z " << maxX << " " << maxY << " " << maxZ <<  std::endl;
      std::cout << "Size X,Y,Z " << maxX-minX << " " << maxY-minY << " " << maxZ-minZ <<  std::endl;      


    uint16_t *img = new uint16_t[int(maxX*maxY)];
    
    for(int i3=0;i3<int(maxX*maxY);i3++)
       img[i3] = 0; 
     
    for(int i2=0; i2<NP; i2++) {
    img[ int(maxX -(X[i2]-1))+  int((maxY -Y[i2])* maxX) ] = int(tabVal[i2]*100);
    }       
  
 // GDCM_NAME_SPACE::Debug::DebugOn();
  
       std::ostringstream str; 
       GDCM_NAME_SPACE::File *file;
       file = GDCM_NAME_SPACE::File::New();       
              
  // Set the image size
        str.str("");
        str << (int)maxX;
        file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns
        str.str("");
        str << (int)maxY;
        file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows
  // Set the pixel type
  //      16; //8, 16, 32
        file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated
        str.str("");
        str << 16; // may be 12 or 16 if componentSize =16
        file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored
        file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit

  // Set the pixel representation // 0/1 , 0=unsigned
        file->InsertEntryString("1",0x0028,0x0103, "US"); // Pixel Representation
  // Set the samples per pixel // 1:Grey level, 3:RGB
        file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel

/*
  // Set Rescale Intercept
        str.str("");
        str << div;  
        file->InsertEntryString(str.str(),0x0028,0x1052,"DS");

  // Set Rescale Slope
        str.str("");
        str << mini;  
        file->InsertEntryString(str.str(),0x0028,0x1053,"DS");
*/
    
    GDCM_NAME_SPACE::FileHelper *fileH;
    fileH = GDCM_NAME_SPACE::FileHelper::New(file);
    // cast is just to avoid warnings (*no* conversion)
    fileH->SetImageData((uint8_t *)img,int(maxX*maxY)*sizeof(uint16_t));
    fileH->SetWriteModeToRaw(); 
    fileH->SetWriteTypeToDcmExplVR();
        
    if( !fileH->Write(dcmImageName))
       std::cout << "Failed for [" << dcmImageName << "]\n"
                 << "           File is unwrittable" << std::endl;

    file->Print();
           
   delete img;
   file->Delete();
   fileH->Delete();  
}