//  $Id: niftiout.cc,v 1.17 2007-11-15 16:19:42 valerio Exp $	

//****************************************************************************
//
// Modification History (most recent first)
// mm/dd/yy  Who  What
//
// 10/17/11  VPL  Intercept and slope are now computed
// 04/25/07  VPL  ndim (a.k.a dim[0]) should be 3 for anatomical, 4 for time
//                series with 1 coil, 5 otherwise.
// 01/31/07  VPL  Protect against missing images at the beginning of series
// 10/24/06  VPL  Set hdr.ndim to 4 when only one coil.
//                Set hdr.nv, hdr.nw, hdr.du, hdr.dv, hdr.dw to 1.
//                Set hdr.scl_slope to 0.
// 03/22/05  VPL  Allow processing of multiple series
//                Add NYU License agreement
// 02/01/05  VPL  
//
//****************************************************************************

#include "niftiout.h"

//****************************************************************************
//
// Purpose: Compute slice distance
//   
// Parameters: nifti_image &hdr            - the header for NIfTI
//             list<DICOMImage> *imageList - pointer to the list of images
//   
// Returns: bool false - error
//   
// Notes: 
//
//****************************************************************************

static bool SliceDistance(nifti_image &hdr, IMAGELIST *imageList)
{
	bool retValue = true;
	IMAGELIST::iterator image = imageList->begin();

	if ( !image->Mosaic() && (imageList->size() != hdr.nz) )
	{
		cerr << endl << "\t**** Number of slices provided, " << imageList->size()
			 << ", does not match value in header, " << hdr.nz << "." << endl
			 << "\t     Value in header can be overridden with \"-s\" flag." << endl;
		return false;
	}
	
	if ( hdr.nz < 2 )                     // Only one slice, so use slice thinkness
	{
		cerr << endl << "\t**** SliceDistance: WARNING: Number of slices is less than two" << endl
			<< "\t**** Using slice thinkness as slice distance" << endl;
		hdr.dz = image->SliceThickness();
    }
	else
	{
		float x = image->SagPos();
		float y = image->CorPos();
		float z = image->TraPos();

		// If this is Mosaic, compute from slices within image
		// Else from first slice from contigous images

		if (image->Mosaic())
		{
			x -= image->SagPos(1);
			y -= image->CorPos(1);
			z -= image->TraPos(1);
		}
		else
		{
			image++;
			x -= image->SagPos();
			y -= image->CorPos();
			z -= image->TraPos();
		}
		
		hdr.dz = (float) sqrt(x*x + y*y + z*z);

		if( hdr.dz < 0.00001)				// if distance between slices is zero, use slice thickness
			hdr.dz = image->SliceThickness();
	}

	return retValue;
}

//****************************************************************************
//
// Purpose: Compute quaternion
//   
// Parameters: nifti_image &hdr            - the header for NIfTI
//             list<DICOMImage> *imageList - pointer to the list of images
//   
// Returns: NONE
//   
// Notes: Copied from n2n2.c by L&R Fleisher
//
//****************************************************************************

static void FindQuaternion(nifti_image &hdr, IMAGELIST *imageList)
{
	// Get the lowest slice and call it k = 0

	IMAGELIST::iterator image = imageList->begin();
	hdr.qto_xyz.m[0][0] = -image->RowSagCos();
	hdr.qto_xyz.m[0][1] = -image->ColSagCos();
	hdr.qto_xyz.m[0][2] = -image->SagNorm();

	hdr.qto_xyz.m[1][0] = -image->RowCorCos();
	hdr.qto_xyz.m[1][1] = -image->ColCorCos();
	hdr.qto_xyz.m[1][2] = -image->CorNorm();

	hdr.qto_xyz.m[2][0] = image->RowTraCos();
	hdr.qto_xyz.m[2][1] = image->ColTraCos();
	hdr.qto_xyz.m[2][2] = image->TraNorm();

	// offset = slicePosition(center) - M*(i,j,k of center)
	// since we are working with the first slice k = 0.

	hdr.qoffset_x = -image->SagPos();
	hdr.qoffset_y = -image->CorPos();
	hdr.qoffset_z =  image->TraPos();

	nifti_mat44_to_quatern(hdr.qto_xyz, &(hdr.quatern_b), &(hdr.quatern_c), &(hdr.quatern_d),
						   NULL, NULL, NULL, NULL, NULL, NULL, &(hdr.qfac));
}

//****************************************************************************
//
// Purpose: Fill in NIfTI header with default values
//   
// Parameters: nifti_image &hdr   - the header to fill
//             IMAGEMAP &imageMap - the image set ordered by timepoint
//   
// Returns: bool false == error
//   
// Notes: Copied from n2n2.c by L&R Fleisher
//
//    Look at nifti_image_write in nifiti1_io.c, and at nifti_1_header 
//    struct in nifti1.h for description of the fields of the data
//    header. The nifti_image struct items initialized below are (mostly)
//    in the order given in nifti1_io.h 
//
//    THE FUNCTION FILLS THE DEFAULT VALUES INTO nifti_image AS FOUND IN 
//    THE INPUT.HERE IS THE LIST OF WHAT ONE MAY WANT TO CHANGE ELSEWHERE:
//        - nifti_type file format should be set according to user specifications
//          notice that only float storage type is supported, see a note below
//
//        - hdr.ndim is set to 5 when there is more than one coil, else it is
//          set to 4.
//
//   FIXME WILL NEVER BE FIXED
//        - hdr.dz is set to distance between slice centers contrary
//          to specification of ANALYZE/NIFTI standard which requires dz
//          be slice thickness! In other words, 
//          slicethicknes != z-grid spacing = dz.
//          This is probably correct and is a mistake in the standard.
//          Also, for 3D acquisitions what is slice thickness?
//          inside find_slice_thinkness, use get_slicethickness()
//        - hdr.dt is set to time between slices
//          If it is zero, it is set to time between volumes.
//          How to define the time of the slice????
//        - datatype is set for int storage type.
//          If other storage types are desired the field should be set 
//          accordingly
//        - slice_code
//        - slice_duration
//        - iname and fname if used, do not belong here
//        - toffset do not understand what it is.
//
//****************************************************************************

bool NIfTICreateHeader(nifti_image &hdr, IMAGEMAP &imageMap)
{
	bool retValue = true;

	IMAGELIST *imageList = imageMap.begin()->second;
	IMAGELIST::iterator image = imageList->begin();

	memset(&hdr, 0, sizeof(nifti_image));
	
	hdr.nx = image->Columns();
	hdr.ny = image->Rows();
	hdr.nz = image->NumSlices();				// number of slices per volume
	hdr.nt = imageMap.size(); 					// number of timepoints
	hdr.nu = 1;                                 // either just 1, or combined into 1

	if ( hdr.nt == 1 )
	{
		hdr.ndim = 3;
	}
	else
	{
		hdr.ndim = 4;
	}
	
	hdr.nv = 1;
	hdr.nw = 1;

	hdr.nvox = (size_t)(hdr.nx * hdr.ny * hdr.nz * hdr.nt);

	hdr.datatype = DT_SIGNED_SHORT;

	// hdr.nbyper should be sizeof(int) for this datatype, but NIFTI does
	// something else
	
	nifti_datatype_sizes(hdr.datatype, &(hdr.nbyper), &(hdr.swapsize));

	hdr.dx = image->XFOV() / ((float)hdr.nx);
	hdr.dy = image->YFOV() / ((float)hdr.ny);

	if (! SliceDistance(hdr, imageList))
    {
		cerr << endl << "\t**** CreateNIfTIHeader: Cannot find slice distance for "
			 << image->ImagePath() << "." << endl << endl;
		return false;
    }

	// hdr.dt is specified as time between slices (what does it mean??)
	// use time between volumes instead, as it is clearer

	if( hdr.nt > 1 )
    {
		hdr.dt = image->RepetitionTime();
    }
	else
    {
		cerr << endl << "\t**** CreateNIfTIHeader: WARNING: number of repetitions is less than two" << endl
			 << "\t**** time between volumes is set to zero" << endl;
		hdr.dt = 0.0;
    }

	hdr.du = 1.0;
	hdr.dv = 1.0;
	hdr.dw = 1.0;

	hdr.scl_slope = image->ImgRescaleSlope();
	hdr.scl_inter = image->ImgRescaleIntercept();

	hdr.cal_min = 0;
	hdr.cal_max = 0;

	// quaternions in scanner-based system
	
	hdr.qform_code = NIFTI_XFORM_SCANNER_ANAT;
	hdr.sform_code = 0;

	hdr.freq_dim  = image->FrequencyDim();
	hdr.phase_dim = image->PhaseDim();
	hdr.slice_dim = 3;
	hdr.slice_code = 0;
	hdr.slice_start = 0;
	hdr.slice_end = hdr.nz - 1;
	hdr.slice_duration = 0.0;

	FindQuaternion(hdr, imageList);

	hdr.toffset = 0.0;    						// Time is since 0:00 of the day
	hdr.xyz_units = NIFTI_UNITS_MM; 			// these are probably correct
	hdr.time_units = NIFTI_UNITS_MSEC;

	// no coil intent in the format specification 
	// should we make it NIFTI_INTENT_VECTOR ?
	
	hdr.intent_code = NIFTI_INTENT_NONE;
	hdr.intent_p1 = 0.0;
	hdr.intent_p2 = 0.0;
	hdr.intent_p3 = 0.0;
	hdr.intent_name[0] = '\0' ;
	hdr.aux_file[0] = '\0';

	hdr.byteorder = nifti_short_order();  /* use host byte order */

	hdr.data = NULL; 

	strncpy(hdr.descrip, image->ACQDescription().c_str(), 79);
	hdr.descrip[79] = '\0';

	return retValue;
}

//****************************************************************************
//
// Purpose: Write the data to the NIfTI file
//   
// Parameters: znzFile fp		  - the (opened) output file
//             nifti_image &hdr   - NIfTI header with all the good info
//             IMAGEMAP &imageMap - map with sorted list of images per timepoint
//   
// Returns: bool false == error
//   
// Notes: 
//
//****************************************************************************

bool NIfTIWriteData(znzFile fp, nifti_image &hdr, IMAGEMAP &imageMap)
{
	// For now start with just one coil mode .... COMB aka norm
	
	int numPix = hdr.nx * hdr.ny;

	hdr.data = (void *)(new unsigned short [numPix]);
	if (hdr.data == NULL)
    {
		cerr << "\t**** NIfTIWriteData: Error: memory allocation failure." << endl << endl;
		return false;
    }
	
	// Assume, order is the same for all repetitions
	
	for (IMAGEMAP::iterator imageList = imageMap.begin(); imageList != imageMap.end(); imageList++)
	{
		// If this a mosaic, load the volume at once

		IMAGELIST::iterator img = imageList->second->begin();
		const bool mosaic = img->Mosaic();
		
		// Loop through slices and read the slice values

		const int numSlices = img->NumSlices();
		
		for (int slice = 0; slice < numSlices; ++slice)
		{
			// This should never happen
			
			if (img == imageList->second->end())
			{
				cerr << "\t**** NIfTIWriteData: Error: "
					 << "Image list larger than number off slices." << endl << endl;

				// clean write buffer
				delete [] static_cast<unsigned short*>(hdr.data);
				
				return false;
			}
			
			// Clear out the array we use for reading
					
			memset(hdr.data, 0, numPix * sizeof(unsigned short));
	
			// All elements are assumed to be the same size or NIFTI/ANALYZE
			// will break
				
			if ((img->Columns() != hdr.nx) || (img->Rows() != hdr.ny))
			{
				cerr << "\t**** NIfTIWriteData: Error: "
					 << "NIFTI does not support acquisition of slices with different FOVs. " << endl
					 << "\t     Aborting conversion" << endl << endl;
				img->FreeVolume();

				// clean write buffer
				delete [] static_cast<unsigned short*>(hdr.data);
				
				return false;
			}
			if (!img->ReadSlice((unsigned short *)hdr.data, slice))
			{
				cerr << "\t**** NIfTIWriteData: Error: "
					 << "Error reading data from DICOM file " << img->ImagePath()
					 << "." << endl << "\t     Aborting conversion" << endl << endl;
				img->FreeVolume();

				// clean write buffer
				delete [] static_cast<unsigned short*>(hdr.data);
				
				return false;
			}

			if ( nifti_write_buffer(fp, hdr.data, numPix * hdr.nbyper) != (numPix * hdr.nbyper))
			{
				cerr << "\t**** NIfTIWriteData: Error: while processing image " << img->ImagePath()
					 << "." << endl << "\t     Data write failed. Aborting conversion" << endl << endl;
				img->FreeVolume();

				// clean write buffer
				delete [] static_cast<unsigned short*>(hdr.data);
				
				return false;
			}

			if (!mosaic)
			{
				img->FreeVolume();
				img++;
			}
		}
		if (img != imageList->second->end()) img->FreeVolume();
	}

	// clean write buffer
	delete [] static_cast<unsigned short*>(hdr.data);

	return true;
}