/***************************************************************************
 *   Copyright (C) 2008 by Yang Li,,,   *
 *   yangli@yangli-desktop   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/


#include <iostream>
#include "DTIRegistrationCLP.h"
#include <string.h>
#include <math.h>
#include <algorithm>
#include <stdio.h>
#include "itkDiffusionTensor3D.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkNrrdImageIO.h"
#include "itksys/Directory.hxx"
#include "itksys/SystemTools.hxx"

int main(int argc, char *argv[]) 
{
  PARSE_ARGS;
  char command[500];
  int PaddedSize[3];
  int MdlMargin[3]; //left margin
  int SubjMargin[3]; //left margin

  	//get the input image information
  	typedef itk::DiffusionTensor3D<float> PixelType;
	typedef itk::Image< PixelType, 3 > ImageType;
	typedef itk::NrrdImageIO ImageIOTypeNrrd;
	typedef itk::ImageFileReader< ImageType > ReaderType;
	
	ReaderType::Pointer MdlReader = ReaderType::New();
        ReaderType::Pointer SubjReader = ReaderType::New();

	ImageIOTypeNrrd::Pointer NrrdIO = ImageIOTypeNrrd::New();
	ImageType::Pointer DTImageMdl = ImageType::New();
	ImageType::Pointer DTImageSubj = ImageType::New();
        
   	MdlReader->SetImageIO( NrrdIO );
	SubjReader->SetImageIO( NrrdIO );

        MdlReader->SetFileName( DTIRegistrationInputModel.c_str() );
	SubjReader->SetFileName( DTIRegistrationInputSubject.c_str() );

	MdlReader->Update();
	SubjReader->Update();
	DTImageMdl = MdlReader->GetOutput();
	DTImageSubj = SubjReader->GetOutput();

	ImageType::RegionType MdlRegion = DTImageMdl->GetLargestPossibleRegion();
	ImageType::RegionType SubjRegion = DTImageSubj->GetLargestPossibleRegion();
	
	ImageType::SizeType MdlSize = MdlRegion.GetSize();
  	ImageType::SpacingType MdlSpacing = DTImageMdl->GetSpacing();

	ImageType::SizeType SubjSize = SubjRegion.GetSize();
  	ImageType::SpacingType SubjSpacing = DTImageSubj->GetSpacing();

        //decide the padding size
   	PaddedSize[0] = std::max(SubjSize[0],MdlSize[0])*2;
	PaddedSize[1] = std::max(SubjSize[1],MdlSize[1])*2;
	PaddedSize[2] = std::max(SubjSize[2],MdlSize[2])*2;
	
	MdlMargin[0] = std::floor((PaddedSize[0]-MdlSize[0])/2);
	MdlMargin[1] = std::floor((PaddedSize[1]-MdlSize[1])/2);
	MdlMargin[2] = std::floor((PaddedSize[2]-MdlSize[2])/2);

	SubjMargin[0] = std::floor((PaddedSize[0]-SubjSize[0])/2);
	SubjMargin[1] = std::floor((PaddedSize[1]-SubjSize[1])/2);
	SubjMargin[2] = std::floor((PaddedSize[2]-SubjSize[2])/2);


   //pad the input images
   std::sprintf(command,"DTIPad %s %d %d %d %d %d %d model.padded.dt", DTIRegistrationInputModel.c_str(), PaddedSize[0], PaddedSize[1], PaddedSize[2], MdlMargin[0]+1, MdlMargin[1]+1, MdlMargin[2]+1);
   std::cout<<command<<"\n";
   system(command);
   std::sprintf(command,"DTIPad %s %d %d %d %d %d %d subject.padded.dt", DTIRegistrationInputSubject.c_str(), PaddedSize[0], PaddedSize[1], PaddedSize[2], SubjMargin[0]+1, SubjMargin[1]+1, SubjMargin[2]+1);
   std::cout<<command<<"\n";
   system(command);
   
   //begin registration
   std::sprintf(command,"tensorgeometryregistration_v2.pl -M model.padded.dt -S subject.padded.dt -X%d,%d,%d -V%f,%f,%f -v%f,%f,%f -F",PaddedSize[0], PaddedSize[1], PaddedSize[2],MdlSpacing[0],MdlSpacing[1],MdlSpacing[2],SubjSpacing[0],SubjSpacing[1],SubjSpacing[2]);
   std::cout<<command<<"\n";
   system(command); 

   //crop the registered
   //system("pwd\n");
   //system("ls");
   std::sprintf(command,"DTConverter Warped.Subj.dt %d %d %d %f %f %f warped.nrrd", PaddedSize[0], PaddedSize[1], PaddedSize[2], MdlSpacing[0], MdlSpacing[1], MdlSpacing[2]);
   std::cout<<command<<"\n";
   system(command);
   std::sprintf(command,"DTICrop warped.nrrd %d %d %d %d %d %d %s",MdlMargin[0]+1,MdlMargin[0]+MdlSize[0],MdlMargin[1]+1,MdlMargin[1]+MdlSize[1],MdlMargin[2]+1,MdlMargin[2]+MdlSize[2], DTIRegistrationOutputRegistered.c_str());
   std::cout<<command<<"\n";
   system(command);
     
   //remove the garbage generated
   std::sprintf(command,"rm model.padded.dt subject.padded.dt Field.1.img.Mdl2Subj Warped.Subj.dt warped.nrrd");
   std::cout<<command<<"\n";
   system(command);
  //std::cout<<DTIRegistrationInputModel.c_str()<<"\n";
  //std::cout<<DTIRegistrationInputSubject.c_str()<<"\n";
  //std::cout<<DTIRegistrationOutputRegistered.c_str()<<"\n";
  //std::cout<<""<<"\n";
  return EXIT_SUCCESS;
}
