// 	$Id: dicomInfo.cc,v 1.36 2008-01-22 19:08:11 valerio Exp $

//============================================================================
// DCM is part of the CTN library, curtesy of
//
//		Mallinckrodt Institute of Radiology
//		Washington University in St. Louis MO
//		http://wuerlim.wustl.edu/
//
//============================================================================

//****************************************************************************
//
// Modification History (most recent first)
// mm/dd/yy  Who  What
//
// 10/17/11  VPL  Add verbose parameter
// 10/17/11  VPL  Add optional haveNumSlices to constructor
// 03/22/05  VPL  Allow processing of multiple series
//                Add NYU License agreement
// 03/15/05  VPL  Adjust some items after the fact if not mosaic
// 03/08/05  VPL  Make sure that we read shadow set only for mosaic images
//
//****************************************************************************

#include "dicomInfo.h"
#include <cstdlib>
#include <iostream>
#include <cstring>

using namespace std;

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Replace white spaces and slashes with underscores in string
///   
/// @param   strng - the string to be "despaced"
///
/// @return  strng - the cleaned string
///
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

string &DeSpace(string &strng)
{
   size_type locsp = 0;
   while ( locsp < strng.size() )
   {
	   if ( (strng[locsp] == ' ') || (strng[locsp] == '/') )
		   strng[locsp] = '_';
	   ++locsp;
   }

   return strng;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Xctor for DICOMImage class
///   
/// @param   path - path to the image
///          verbose - show all messages
///          haveNumSlices - don't need to read the number of slices from file ?
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

DICOMImage::DICOMImage(const char *path, bool verbose, bool haveNumSlices) :
	error_(false),
	handle_(NULL),
	imageNum_(0),
	numGroups_(0),
	mosaic_(false),
	d3_(false),
	sliceThickness_(0.0),
	sliceSpacing_(0.0),
	acqNum_(1),
	acqEchoNum_(0),
	columns_(-1),
	rows_(-1),
	xFOV_(-1),
	yFOV_(-1),
	repTime_(0),
	phaseEncoding_(0),
	intercept_(0.),
	slope_(0.),
	volumeData_(NULL),
	manufacturer_(UNDEFINED),
	shadowSet_(NULL),
	verbose_(verbose),
	numSlices_(0)
{
	if ( ! OpenFile(path) )
	{
		error_ = true;
		string msg = string("DICOMImage::DICOMImage: Cannot open file: ") + string(path);
		cerr << endl << "\t**** " << msg << endl;
	}
	else
	{
		imgPath_ = path;

		// Get the info we care about (if any fail, the error flag will be set)
	   
		GetMosaic()
			&& GetManufacturer()
			&& GetImageNumber()
			&& GetACQNumber()
			&& GetACQEchoNumber()
			&& GetACQMatrix()
			&& GetD3()
			&& GetSliceThickness()
			&& (haveNumSlices || GetNumSlices())
			&& GetImgOrientation()
			&& GetFOV()
			&& GetColumns()
			&& GetRows()
			&& MosaicAdjustments()
			&& GetPosNorm()
			&& GetSliceSpacing()
			&& GetComments()
			&& GetRepetitionTime()
			&& GetPhaseEncoding()
			&& GetImgRescaleIntercept()
			&& GetImgRescaleSlope()
			&& GetGroupsInfo()
			&& GetSeriesDescription()
			&& GetSeriesNumber()
			&& GetSeriesID()
			&& GetACQDescription()
			&& GetSubjectID()
			&& GetSubjectName()
			&& GetSeriesDate()
			&& GetSeriesTime()
			&& GetAccession();
	}

	if (handle_ != NULL)
		DCM_CloseObject(&handle_);
	handle_ = NULL;
}

//////////////////////////////////////////////////////////////////////////////
///
///
/// Purpose: Fix a couple of parameters for mosaics
///   
/// @param NONE
///
/// Returns: always true
/// 
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::MosaicAdjustments()
{
	if (mosaic_)
	{
		// The total number of entries in the mosaic will be the smallest perfect square
		// greater than the number of slices. Each side is the sqaure root of it

		mosaicGridSize_ = int(ceil(sqrt(float(NumSlices()))));

		columns_ /= mosaicGridSize_;
		rows_ /= mosaicGridSize_;
		xFOV_ /= mosaicGridSize_;
		yFOV_ /= mosaicGridSize_;
	}

	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
///
/// Purpose: Open image file
///   
/// @param path - path to image
///
/// Returns: false == error
/// 
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::OpenFile(const char *path)
{
	bool retValue = true;
	unsigned long options = DCM_ORDERLITTLEENDIAN | DCM_FORMATCONVERSION | DCM_VRMASK;
	
	// Try opening as PART10, if it fails it's might be bcause it does not have
	// a preable and the try it that way
	
	if ( DCM_OpenFile(path, options | DCM_PART10FILE, &handle_) != DCM_NORMAL )
	{
		DCM_CloseObject(&handle_);
		COND_PopCondition(TRUE);
		
		if ( DCM_OpenFile(path, options, &handle_) != DCM_NORMAL )
			retValue = false;
	}

	return retValue;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the name of the manufacturer (0008,0070)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: For now we know only of SIEMENS and GE
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetManufacturer()
{
   char *data = NULL;
   if (GetData(DCM_IDMANUFACTURER, &data) == 0)
   {
	  string msg("DICOMImage::GetManufacturer: Cannot get value of manufacturer in image ");
	  msg += imgPath_;
	  cerr << endl << "\t**** " << msg << endl;
	  error_ = true;
	  if (data != NULL) delete [] data;
	  return false;
   }
   
   if (strcmp(data, "SIEMENS") == 0)
	   manufacturer_ = SIEMENS;
   else if (strcmp(data, "GE MEDICAL SYSTEMS") == 0)
	   manufacturer_ = GE;

   delete [] data;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the acquisition repetiton time (0018,0080)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: DICOM stores the time in a string formated as hhmmss.fraction
//         we will store it in double representing the seconds.
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetRepetitionTime()
{
   // The number is really stored as a string
  
   char *data = NULL;
  
   if (GetData(DCM_ACQREPETITIONTIME, &data) == 0)
   {
	  string msg("DICOMImage::GetRepetitonTime: Cannot get value of acquisition repetition time in image ");
	  msg += imgPath_;
	  cerr << endl << "\t**** " << msg << endl;
	  error_ = true;
	  if (data != NULL) delete [] data;
	  return false;
   }

   repTime_ = atoi(data);
   delete [] data;   
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the phase encoding (0018,1312)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: DICOM stores the time in a string as either ROW, COL or OTHER. In
///        the case of OTHER, we will default to ROW.
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetPhaseEncoding()
{
   // The number is really stored as a string
  
   char *data = NULL;
  
   if (GetData(DCM_ACQPHASEENCODINGDIRECTION, &data) == 0)
   {
	  string msg("DICOMImage::GetPhaseEncoding: Cannot get value of phase encoding direction time in image ");
	  msg += imgPath_;
	  cerr << endl << "\t**** " << msg << endl;
	  error_ = true;
	  if (data != NULL) delete [] data;
	  return false;
   }

   if ( strcmp(data, "COL") == 0 )
	   phaseEncoding_ = 2;
   else
	   phaseEncoding_ = 1;

   delete [] data;   
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the image rescale intercept (0028,1052)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: This is an optional parameter. If not present, it will be set
///        to 0 and the return is always true.
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetImgRescaleIntercept()
{
   // The number is really stored as a string
  
   char *data = NULL;
  
   if (GetData(DCM_IMGRESCALEINTERCEPT, &data) == 0)
   {
          string msg("DICOMImage::GetImgRescaleIntercept: image ");
          msg += imgPath_;
		  msg += "does not contain image rescale intercept parameter, intercept set to 0";
          if ( verbose_ ) cerr << endl << "\t**** " << msg << endl;
          error_ = false;
		  intercept_ = 0.;
   }
   else
	   intercept_ = atof(data);

   if (data != NULL) delete [] data;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the image rescale slope (0028,1053)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: This is an optional parameter. If not present, it will be set
///        to 0 and the return is always true.
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetImgRescaleSlope()
{
   // The number is really stored as a string
  
   char *data = NULL;
  
   if (GetData(DCM_IMGRESCALESLOPE, &data) == 0)
   {
	   string msg("DICOMImage::GetImgRescaleSlope: image ");
          msg += imgPath_;
		  msg += "does not contain image rescale slope parameter, slope set to 0";
          if ( verbose_ ) cerr << endl << "\t**** " << msg << endl;
          error_ = false;
		  COND_PopCondition(TRUE);				// Clear error from stack
		  slope_ = 0.;
   }
   else
	   slope_ = atof(data);

   if (data != NULL) delete [] data;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the series instance id (0020, 000e)
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSeriesID()
{
   char *strn = NULL;
  
   if (GetData(DCM_RELSERIESINSTANCEUID, &strn) == 0)
   {
	  string msg("DICOMImage::GetSeriesID: Cannot get value of series id in image ");
	  msg += imgPath_;
	  cerr << endl << "\t**** " << msg << endl;
	  error_ = true;
	  if (strn != NULL) delete [] strn;
	  return false;
   }

   seriesID_ = strn;
   delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Build an acquisition description string like SPM5
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: SPM5 concatenates a number of different parameters.
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetACQDescription()
{
   char *strn = NULL;
  
   if (GetData(DCM_ACQTRANSMITTINGCOIL, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
	   acqDescription_ = strn;

   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }
   
   if (GetData(DCM_ACQMRACQUISITIONTYPE, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
	   acqDescription_ += string(" ") + strn;

   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }

   if (GetData(DCM_ACQSCANNINGSEQUENCE, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
	   acqDescription_ += string(" ") + strn;

   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }

   if (GetData(DCM_ACQREPETITIONTIME, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
	   acqDescription_ += string(" TR=") + strn;

   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }

   if (GetData(DCM_ACQECHOTIME, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
	   acqDescription_ += string("/TE=") + strn;

   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }

   if (GetData(DCM_ACQFLIPANGLE, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
	   acqDescription_ += string("/FA=") + strn;

   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }

   tm acqDatTime;
   memset(&acqDatTime, 0, sizeof(acqDatTime));
   
   if (GetData(DCM_IDACQUISITIONDATE, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
   {
	   acqDatTime.tm_mday = atoi(&(strn[6]));
	   strn[6] = '\0';
	   acqDatTime.tm_mon = atoi(&(strn[4])) - 1;
	   strn[4] = '\0';
	   acqDatTime.tm_year = atoi(strn) - 1900;
   }

   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }

   if (GetData(DCM_IDACQUISITIONTIME, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
   {
	   // In theory this could be just the time or the date + time
	   
	   char *dec = strchr(strn, '.');
	   if ( dec != NULL )
	   {
		   *dec = '\0';
		   dec -= 2;
		   acqDatTime.tm_sec = atoi(dec);
		   *dec = '\0';
		   dec -= 2;
		   acqDatTime.tm_min = atoi(dec);
		   *dec = '\0';
		   dec -= 2;
		   acqDatTime.tm_hour = atoi(dec);
	   }
   }
   
   if ( strn != NULL )
   {
	   delete [] strn;
	   strn = NULL;
   }

   if ( (acqDatTime.tm_mday != 0) || (acqDatTime.tm_mon != 0) || (acqDatTime.tm_year != 0)
		|| (acqDatTime.tm_sec != 0) || (acqDatTime.tm_min != 0) || (acqDatTime.tm_hour != 0) )
   {
	   strn = new char[32];
	   strftime(strn, 31, "%d-%b-%Y %H:%M:%S", &acqDatTime);
	   acqDescription_ += string(" ") + strn;
	   delete [] strn;
	   strn = NULL;
   }
   
   if ( mosaic_ ) acqDescription_ += string(" Mosaic");

   if ( acqDescription_.length() == 0 )
	   acqDescription_ = "DINIfTI";

   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the subject ID (0010, 0020)
///   
/// @param   NONE 
///   
/// @return  always true
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSubjectID()
{
   char *strn = NULL;
  
   if (GetData(DCM_PATID, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  subjID_ = "";
   }
   else
	   subjID_ = strn;
   delete [] strn;
   DeSpace(subjID_);
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the subject name (0010, 0010)
///   
/// @param   NONE 
///   
/// @return  always true
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSubjectName()
{
   char *strn = NULL;
  
   if (GetData(DCM_PATNAME, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  subjName_ = "";
   }
   else
	   subjName_ = strn;
   delete [] strn;
   DeSpace(subjName_);
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the series date (0008, 0021)
///   
/// @param   NONE 
///   
/// @return  always true
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSeriesDate()
{
   char *strn = NULL;
  
   if (GetData(DCM_IDSERIESDATE, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  seriesDate_ = "";
   }
   else
	   seriesDate_ = strn;
   delete [] strn;
   DeSpace(seriesDate_);
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the series time (0008, 0031)
///   
/// @param   NONE 
///   
/// @return  always true
///   
/// Notes: Cut off past the minutes
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSeriesTime()
{
   char *strn = NULL;
  
   if (GetData(DCM_IDSERIESTIME, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  seriesTime_ = "";
   }
   else
   {
	   strn[4] = '\0';
	   seriesTime_ = strn;
   }
   
   delete [] strn;
   DeSpace(seriesTime_);
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value accession number (0008, 0050)
///   
/// @param   NONE 
///   
/// @return  always true
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetAccession()
{
   char *strn = NULL;
  
   if (GetData(DCM_IDACCESSIONNUMBER, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  accession_ = "";
   }
   else
	   accession_ = strn;
   delete [] strn;
   DeSpace(accession_);
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the series description (0008, 103e)
///   
/// @param   NONE 
///   
/// @return  always true
///   
/// Notes: If we can't find the description, make something up
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSeriesDescription()
{
   char *strn = NULL;
  
   if (GetData(DCM_IDSERIESDESCR, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  seriesDescription_ = "";
   }
   else
	   seriesDescription_ = strn;
   delete [] strn;
   DeSpace(seriesDescription_);
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the image number (0020, 0013)
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetImageNumber()
{
   // The number is really stored as a string
  
   char *strn = NULL;
  
   if (GetData(DCM_RELIMAGENUMBER, &strn) == 0)
   {
	  string msg("DICOMImage::ImageNumber: Cannot get value of image number in image ");
	  msg += imgPath_;
	  cerr << endl << "\t**** " << msg << endl;
	  error_ = true;
	  if (strn != NULL) delete [] strn;
	  return false;
   }

   // Convert string to numerical value

   imageNum_ = atoi(strn);
   delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the series number (0020, 0011)
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSeriesNumber()
{
   // The number is really stored as a string
  
   char *strn = NULL;
  
   if (GetData(DCM_RELSERIESNUMBER, &strn) == 0)
   {
	  string msg("DICOMImage::SeriesNumber: Cannot get value of image number in image ");
	  msg += imgPath_;
	  cerr << endl << "\t**** " << msg << endl;
	  error_ = true;
	  if (strn != NULL) delete [] strn;
	  return false;
   }

   // Convert string to numerical value

   seriesNumber_ = strn;
   delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the relative acquisition number (0020, 0012)
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetACQNumber()
{
   // The number is really stored as a string
  
   char *strn = NULL;
  
   if (GetData(DCM_RELACQUISITIONNUMBER, &strn) == 0)
   {
	   // Defaults to 1

	   COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
   {
	   // Convert string to numerical value

	   acqNum_ = atoi(strn);
   }
   
   if (strn != NULL) delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the acquisition echo number (0018, 0086)
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: This applies to field maps
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetACQEchoNumber()
{
   // The number is really stored as a string
  
   char *strn = NULL;
  
   if (GetData(DCM_ACQECHONUMBER, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  acqEchoNum_ = 0;
   }
   else
   {
	   // Convert string to numerical value

	   acqEchoNum_ = atoi(strn);
   }
   
   if (strn != NULL) delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the acquisition matrix (0018, 1310)
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetACQMatrix()
{
   U16 values[4] = {0,0,0,0};	
   if (GetData(DCM_ACQACQUISITIONMATRIX, values) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }

   // Values are Phase Columns/Phase Rows/Frequency Column/Frequency Row
   // two will always be 0, the other two will give the accurate count

   if ( values[0] == 0 )
   {
	   acqRows_ = values[1];
	   acqColumns_ = values[2];
   }
   else
   {
	   acqColumns_ = values[0];
	   acqRows_ = values[3];
   }
   
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the image orientation (0020, 0037)
///   
/// @param   NONE 
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetImgOrientation()
{
	// The number is really stored as a string
  
	char *strn = NULL;
  
	if (GetData(DCM_RELIMAGEORIENTATIONPATIENT, &strn) == 0)
	{
		// Default to 1,1,1 for both rows and columns

		rowCosines_.assign(3,1);
		colCosines_.assign(3,1);
		COND_PopCondition(TRUE);				// Clear error from stack
	}
	else
	{
		// Convert string to numerical value
		// Two sets of direction cosines (x, y, z) for first row and first column
		// with respect to the patient, separated by back slashes.
		// First set is row values, second the colum values

		char *tok = strtok(strn, "\\");
		rowCosines_.push_back(atof(tok));
		tok = strtok(NULL, "\\");
		rowCosines_.push_back(atof(tok));
		tok = strtok(NULL, "\\");
		rowCosines_.push_back(atof(tok));

		tok = strtok(NULL, "\\");
		colCosines_.push_back(atof(tok));
		tok = strtok(NULL, "\\");
		colCosines_.push_back(atof(tok));
		tok = strtok(NULL, "\\");
		colCosines_.push_back(atof(tok));
	}
	
	if (strn != NULL) delete [] strn;
	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the mosaic flag (0008, 0008)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetMosaic()
{
   // We have to parse what the header contains
  
   char *strn = NULL;
  
   if (GetData(DCM_IDIMAGETYPE, &strn) == 0)
   {
	  string msg("DICOMImage::Mosaic: Cannot get value of image type in image ");
	  msg += imgPath_;
	  cerr << endl << "\t**** " << msg << endl;
	  error_ = true;
	  if (strn != NULL) delete [] strn;
	  return false;
   }

   // See if the word "MOSAIC" is there

   mosaic_ = (strstr(strn, "MOSAIC") != NULL);
   delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the acquisition type flag (0018, 0023)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetD3()
{
	// We have to parse what the header contains
  
	char *strn = NULL;
  
   if (GetData(DCM_ACQMRACQUISITIONTYPE, &strn) == 0)
   {
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  d3_ = false;
   }
   else
   {
	   // See if the string "3D" is there

	   d3_ = (strstr(strn, "3D") != NULL);
   }
   
   if (strn != NULL) delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the slice thickness (0018,0050)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSliceThickness()
{
   // The number is really stored as a string
  
   char *strn = NULL;
  
   if (GetData(DCM_ACQSLICETHICKNESS, &strn) == 0)
   {
	  string msg("DICOMImage::SliceThickness: Cannot get value of slice thickness in image ");
	  error_ = false;
	  COND_PopCondition(TRUE);				// Clear error from stack
	  sliceThickness_ = 1.;
   }
   else
   {
	   // Convert string to numerical value

	   sscanf(strn, "%f", &sliceThickness_);
   }
   
   if (strn != NULL) delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the slice spacing (0018,0088)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: If images are contigous, no spacing info will be available
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetSliceSpacing()
{
   // The number is really stored as a string
  
   char *strn = NULL;
  
   if (GetData(DCM_ACQSLICESPACING, &strn) == 0)
   {
	  sliceSpacing_ = 0;
	  error_ = false;
	  if (strn != NULL) delete [] strn;
	  COND_PopCondition(TRUE);				// Clear error from stack
   }
   else
   {
	  // Convert string to numerical value

	  sscanf(strn, "%f", &sliceSpacing_);
   }

   delete [] strn;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for the comments (0020,4000)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: This is used by raw images to store (that's why we don't care about
/// catching an error).
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetComments()
{
   char *data = NULL;
   if (GetData(DCM_RELIMAGECOMMENTS, &data) != 0)
	 comments_ = data;
   else
	  COND_PopCondition(TRUE);				// Clear error from stack
   
   error_ = false;
   if (data != NULL) delete [] data;
   return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for columns (0028, 0011)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetColumns()
{
	U16 data;
	
	if (GetData(DCM_IMGCOLUMNS, &data) == 0)
	{
		columns_ = 0;
		error_ = true;
		return false;
	}
	columns_ = data;
	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the value for rows (0028, 0010)
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: 
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetRows()
{
	U16 data;
	if (GetData(DCM_IMGROWS, &data) == 0)
	{
		rows_ = 0;
		error_ = true;
		return false;
	}

	rows_ = data;
	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Compute values of FOV's
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: Read pixel spacing values (0028, 0030) and use previously read
/// 	   columns and rows
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetFOV()
{
	// Values are stored in a string separated by backslashes
  
	char *strn = NULL;
  
	if (GetData(DCM_IMGPIXELSPACING, &strn) == 0)
	{
		xFOV_ = yFOV_ = 0;
		error_ = true;
		if (strn != NULL) delete [] strn;
		COND_PopCondition(TRUE);				// Clear error from stack
	}
	else
	{
		// Extract numerical values and use them to compute FOVs
		// Make sure we have read the values for rows and columns (which should
		// be positive

		if ( rows_ < 0 )
		{
			if ( GetRows() )
			{
				if ( columns_ < 0 )
					if ( ! GetColumns() )
					{
						delete [] strn;
						return false;
					}
			}
			else
			{
				delete [] strn;
				return false;
			}
		}
		
		char *tok = strtok(strn, "\\");
		yFOV_ = atof(tok) * float(rows_);
		tok = strtok(NULL, "\\");
		xFOV_ = atof(tok) * float(columns_);
		delete [] strn;
	}

	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get the shadow set data
///   
///   
/// @return  data - pointer to data
///   
/// Notes: Read it only once. We always return a copy which the user must
/// delete. 
///
//////////////////////////////////////////////////////////////////////////////

char *DICOMImage::GetShadowSet()
{
	static U32 len;
	
	if (shadowSet_ == NULL)
	{
		if (DCM_GetElementSize(&handle_, DCM_MAKETAG(0x0029,0x1020), &len) != DCM_NORMAL)
		{
			COND_PopCondition(TRUE);				// Clear error from stack
			return NULL;
		}
  
		char *storage = NULL;
		if ((len = GetData(DCM_MAKETAG(0x0029,0x1020), &storage)) == 0)
		{
			if (storage != NULL) delete [] storage;
			return NULL;
		}

		// Find the beginning of what we're interested in
		// Start from the back of the buffer and find matching
		// "ASCCONV END"/"ASCCONV BEGIN" block
		// Cannot use standard string search methods 'cause this is full of NULLs

		U32 startp = len - 1;
		char *strn = &(storage[startp]);
  
		while ( startp > 0 )
		{
			if ((*strn == '#') && (strncmp(strn, "### ASCCONV END", 15) == 0))
				break;
			--strn;
			--startp;
		}

		while ( startp > 0 )
		{
			if ((*strn == '#') && (strncmp(strn, "### ASCCONV BEGIN", 17) == 0))
				break;
			--strn;
			--startp;
		}

		len = len - startp + 1;
		
		if (len > 0)
		{
			shadowSet_ = new char [len];
			memcpy(shadowSet_, strn, sizeof(char)*(len));
		}

		delete [] storage;
	}
	else
		len = strlen(shadowSet_);
	
	char *retValue = new char [len];
	memcpy(retValue, shadowSet_, sizeof(char)*(len));
	return retValue;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Read the Siemens special data (0029, 1020) and get the number of
///          groups in this directory plus each group's size and starting
///          image number.
///   
/// @param   
///   
/// @return  false if error
///   
/// Notes: This looks and feels like a hack .... ;-(
///        We only know how to deal with Siemens shaodw set
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetGroupsInfo()
{
	char *storage = GetShadowSet();

	if ( (manufacturer_ != SIEMENS) || (storage == NULL) )
	{
		numGroups_ = 1;
		groupSize_.resize(1);
		groupSize_[0] = 1;
		groupStartImage_.resize(1);
		groupStartImage_[0] = 1;
		return true;
	}

	bool found = false;

	// Find all of the info of the form "sGroupArray.asGroup[nn].nSize" and
	// store the last "nn" value found.
  
	char *tok = strtok(storage, "\n");
	while ( tok )
	{
		// Make sure we don't go past the end
		if ((*tok == '#') && (strncmp(tok, "### ASCCONV END ###", 18) == 0))
			break;
		
		if ( ! strncmp(tok, "sGroupArray.asGroup[", 20) )
		{
			char *valuep = strstr(tok, "[");
			int groupNum;
			sscanf(valuep, "[%i].", &groupNum);
			if ( groupNum >= numGroups_ )
			{
				numGroups_ = groupNum + 1;                // The indexes start at 0, we want the total #
				if ( groupSize_.size() < numGroups_ )
				{
					groupSize_.resize(numGroups_);
					groupStartImage_.resize(numGroups_);
				}
			}
			
			valuep = strstr(tok, "= ");
			int parValue;
			sscanf(valuep, "= %d", &parValue);
			char *parp = strstr(tok, "]");
			if ( ! strncmp(parp, "].nSize", 7) )
			{
				groupSize_[groupNum] = parValue;
			}
			else if ( ! strncmp(parp, "].nLow", 6) )
			{
				// What is stored is the last image number of th previous group
				
				groupStartImage_[groupNum] = parValue + 1; 
			}
			
			found = true;
		}
		  
		tok = strtok(NULL, "\n");
	}

	if ( ! found )
	{
		string msg("DICOMImage::GetGroupsInfo: Cannot get value of Siemens info in image ");
		msg += imgPath_;
		cerr << endl << "\t**** " << msg << endl;
		error_ = true;
		delete [] storage;
		storage = NULL;
		return false;
	}

	// The first group has no start image number, set it to 1 for consistency

	groupStartImage_[0] = 1;
	
	delete [] storage;
	storage = NULL;
	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: If this is a Siemens image read the shadow set (0029, 1020) and
//           get the number of slices in this file.
///   
/// @param   NONE
///   
/// @return  Always true
///   
/// Notes: If this is not mosaic, we will have to fix it later.
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetNumSlices()
{
	if ( manufacturer_ == SIEMENS )
	{
		char *storage = GetShadowSet();
		if (storage == NULL)
		{
			numSlices_ = 1;
			error_ = true;
		}
		else
		{
			bool found = false;
	  
			// Believe it or not, we are looking only for one thing ....
	  
			char *tok = strtok(storage, "\n");
			while (tok)
			{
				if ( mosaic_ || !d3_ )
				{
					if (! strncmp(tok, "sSliceArray.lSize", 17))
					{
						char *valuep = strstr(tok, "=");
						sscanf(valuep, "=%i", &numSlices_);
						found = true;
						break;
					}
				}
				else
				{
					if (! strncmp(tok, "sSliceArray.asSlice[0].dThickness", 33))
					{
						char *valuep = strstr(tok, "=");
						float fval;
						sscanf(valuep, "=%f", &fval);
						numSlices_ = int(0.5 + fval / sliceThickness_);
						found = true;
						break;
					}
				}
		
				tok = strtok(NULL, "\n");
			}

			if (! found)
				error_ = true;

			delete [] storage;
			storage = NULL;
		}
	}
	else
	{
		char *strn = NULL;
		if (GetData(DCM_RELIMAGESINACQUISITION, &strn) == 0)
		{
			numSlices_ = 1;
			COND_PopCondition(TRUE);				// Clear error from stack
			error_ = true;
		}
		else
		{
			numSlices_ = atoi(strn);
		}
		
		delete [] strn;
	}

	if ( error_ )
	{
		error_ = false;
		string msg("DICOMImage::NumSlices: Cannot get number of slices from ");
		msg += imgPath_;
		cerr << endl << "\t**** " << msg << endl
			 << "\t**** Setting value to 1 (can be changed with '-s' option on command line.)" << endl;
	}

	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Store sag/cor/tra positons and normals.
///   
/// @param   NONE
///   
/// @return  false if error
///   
/// Notes: If this is SIEMENS, we'll have a list of values which we will get from
///        the Siemens shadow set, otherwise read it from the DICOM header.
///        The signs on the normal will have to be adjusted once we have the
/// 	   full set using AdjustSet().
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetPosNorm()
{
	sagNorm_ = rowCosines_.at(1) * colCosines_.at(2) - rowCosines_.at(2) * colCosines_.at(1);
	corNorm_ = rowCosines_.at(2) * colCosines_.at(0) - rowCosines_.at(0) * colCosines_.at(2);
	traNorm_ = rowCosines_.at(0) * colCosines_.at(1) - rowCosines_.at(1) * colCosines_.at(0);

	char *storage = NULL;
	
	if ( (mosaic_) && (manufacturer_ == SIEMENS) && ((storage = GetShadowSet()) != NULL) )
		return GetPosNormSiemens(storage);
	else
	{
		char *str = NULL;
		if (GetData(DCM_RELIMAGEPOSITIONPATIENT, &str) == 0)
		{
			error_ = true;
			if (str != NULL) delete [] str;
			COND_PopCondition(TRUE);				// Clear error from stack
			return false;
		}
		else
		{
			// Extract numerical values and push them on the lists

			char *tok = strtok(str, "\\");
			float posValue = atof(tok);
			sagPos_.push_back(posValue);

			tok = strtok(NULL, "\\");
			posValue = atof(tok);
			corPos_.push_back(posValue);
			
			tok = strtok(NULL, "\\");
			posValue = atof(tok);
			traPos_.push_back(posValue);
			
			delete [] str;
		}
	}
	
	return true;
}


//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Store sag/cor/tra positons and normals.
///   
/// @param   char *storage - the shadow set
///   
/// @return  false if error
///   
/// Notes: Get a list of values from the Siemens shadow set.
///        Siemens stores the value of the middle of the slice. The DICOM
///        standard specifies the middle of the upper-left-hand pixel, so we
///        have to adjust half FOV minus half pixel.
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetPosNormSiemens(char *storage)
{
	// Find the info of the form "sSliceArray.asSlice[NN].sNormal.dSag/dCor/dTra"
	// and "sSliceArray.asSlice[NN].sPosition.dSag/dCor/dTra" for each slice (if
	// this a Mosaic, there will be more than one slice)
   
	char *tok = strtok(storage, "\n");

	// Since the number of entries varies (sic), find the start of
	// sSliceArray.asSlice and stop parsing at the end of it
   
	bool parsed = false;
	float *values = new float [numSlices_ * 3];
	memset(values, 0, sizeof(float) * numSlices_ * 3);

	if (tok)
		do 
		{
			if (! strncmp(tok, "sSliceArray.asSlice[", 20))
			{
				parsed = true;

				// Find out what slice number we're talking about

				char *valuep = strstr(tok, "[");
				int sliceNum;
				sscanf(valuep, "[%i].", &sliceNum);

				// For mosaic images load all slice info,
				// for others we want the slice that corresponds to this image
				// Always store at least the first, just in case there is only one

				if ( (mosaic_ && (sliceNum < numSlices_))
					 || (!mosaic_ && (sliceNum == (imageNum_ - 1)))
					 || (sliceNum == 0) )
				{
					if (strstr(tok, "].sPosition.d") == NULL)
						continue;
				
					// Figure out which of the three axis it is

					int axis = 0;
					if (strstr(tok, ".dCor") != NULL)
						axis = 1;
					else if (strstr(tok, ".dTra") != NULL)
						axis = 2;

					// Time to read the value and store it
					
					valuep = strstr(tok, "=");

					// For each slice, we will order the entries by axis

					int index = axis;
					if (mosaic_)
						index += sliceNum * 3;
					sscanf(valuep, "=%f", &(values[index]));
				}
			}
			else
				if (parsed) break;
		
		} while (tok = strtok(NULL, "\n"));

	// Compute adjustment values
	
	// First make sure that we've read FOVs (which are always positive) 

	if ((xFOV_ < 0) || (yFOV_ < 0))
		if (!GetFOV())
		{
			delete [] storage;
			delete [] values;
			return false;
		}

	// Siemens draws the images starting at x=-FOV/2, y=-FOV/2
	// therefore, we need to move FOV/2 from the center of the slice
	float rowProp = ((float)rows_ / 2)    * ((float)xFOV_ / rows_);
	float colProp = ((float)columns_ / 2) * ((float)yFOV_ / columns_);
	
	float sagAdjust = - rowProp * rowCosines_.at(0) - colProp * colCosines_.at(0);
	float corAdjust = - rowProp * rowCosines_.at(1) - colProp * colCosines_.at(1);
	float traAdjust = - rowProp * rowCosines_.at(2) - colProp * colCosines_.at(2);
	
	// Now that we have our values, let's store them ...

	for (int slice = 0, index = 0; slice < numSlices_; ++slice)
	{
		sagPos_.push_back(values[index] + sagAdjust);
		++index;
		corPos_.push_back(values[index] + corAdjust);
		++index;
		traPos_.push_back(values[index] + traAdjust);
		++index;
	}

	delete [] storage;
	storage = NULL;
	delete [] values;
	values = NULL;
	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get data from image
///   
/// @param   tag    - which element
/// @param   storage - where to store info
///   
/// @return  U32 len - length of item or 0 if error
///   
/// Notes: Overloaded string version
///
//////////////////////////////////////////////////////////////////////////////

U32 DICOMImage::GetData(DCM_TAG tag, char **storage)
{
	DCM_ELEMENT element;
	U32 len;
  
	if (DCM_GetElement(&handle_, tag, &element) != DCM_NORMAL
		|| DCM_GetElementSize(&handle_, tag, &len) != DCM_NORMAL)
	{
		error_ = true;
		return 0;
	}

	void *ctx = NULL;
	if (*storage == NULL)
		*storage = new char [len+1];
   
	element.d.string = *storage;

	U32 newLen;
	if (DCM_GetElementValue(&handle_, &element, &newLen, &ctx) != DCM_NORMAL)
	{
		error_ = true;
		return 0;
	}

	if (newLen > len)
	{
		string msg("DICOMImage::GetData: Inconsistent data size");
		cerr << endl << "\t**** " << msg << endl;
		error_ = true;
		return 0;
	}

	// Terminate and strip trailing blanks

	(*storage)[len] = '\0';

	if (len > 0)
	{
		while ((len > 0) && ((*storage)[--len] == ' '))
			(*storage)[len] = '\0';
  
		return len+1;
	}
	else
		return 0;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get data from image
///   
/// @param   tag  - which element
/// @param   storage - where to store info
///   
/// @return  bool - false if error
///   
/// Notes: Overloaded U32 version
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetData(DCM_TAG tag, U32 *storage)
{
	DCM_ELEMENT element;
	if (DCM_GetElement(&handle_, tag, &element) != DCM_NORMAL)
	{
		error_ = true;
		return false;
	}

	U32 len;
	void *ctx = NULL;
	element.d.ul = storage;
	if (DCM_GetElementValue(&handle_, &element, &len, &ctx) != DCM_NORMAL)
	{
		error_ = true;
		return false;
	}
  
	return true;
}

//////////////////////////////////////////////////////////////////////////////
///
/// Purpose: Get data from image
///   
/// @param   tag  - which element
/// @param   storage - where to store info
///   
/// @return  bool - false if error
///   
/// Notes: Overloaded U16 version
///
//////////////////////////////////////////////////////////////////////////////

bool DICOMImage::GetData(DCM_TAG tag, U16 *storage)
{
	DCM_ELEMENT element;
	if (DCM_GetElement(&handle_, tag, &element) != DCM_NORMAL)
	{
		error_ = true;
		return false;
	}

	U32 len;
	void *ctx = NULL;
	element.d.us = storage;
	if (DCM_GetElementValue(&handle_, &element, &len, &ctx) != DCM_NORMAL)
	{
		error_ = true;
		return false;
	}
  
	return true;
}

//****************************************************************************
//
// Purpose: Read the pixel data for one slice
//   
// Parameters: U16 *data  - where we store it
//             int slice  - which slice in set (defaults to 0)
//   
// Returns: bool false == error
//   
// Notes: If _handle == NULL we need to open the file
//
//****************************************************************************

bool DICOMImage::ReadSlice(U16 *data, int slice)
{
	// Maybe we've read the volume, maybe not
		
	if ( volumeData_ == NULL )
		if ( ! ReadVolume() )
			return false;
	
	if (slice >= numSlices_)
	{
		error_ = true;
		string msg = string("DICOMImage::ReadSlice: Trying to read past end of file: ") + imgPath_;
		cerr << endl << "\t**** " << msg << endl;
		return false;
	}

	// Find the slice in the set
	// If this is not mosaic, there is only one slice

	if (mosaic_)
	{
		// Mosaic side in pixels

		int mosPix = mosaicGridSize_ * columns_;
		
		// Determine the mosaic row and column for the slice and the increment

		int mosRow = slice / mosaicGridSize_;
		int mosCol = slice % mosaicGridSize_;
		int mosInc = mosPix - columns_;

		// We can now copy the data from the correct spot
		
		int mosInd = mosCol * columns_ + mosPix * mosRow * rows_;
		int datInd = 0;
		for (int row = 0; row < rows_; ++row, mosInd += mosInc)
			for (int col = 0; col < columns_; ++col, ++mosInd, ++datInd)
				data[datInd] = volumeData_[mosInd];
	}
	else
	{
		memcpy((void *)data, (const void *)volumeData_, columns_ * rows_ * sizeof(U16));
	}
	
	return true;
}

//****************************************************************************
//
// Purpose: Read the whole volume and store it
//   
// Parameters: NONE
//   
// Returns: bool false = error
//   
// Notes: 
//
//****************************************************************************

bool DICOMImage::ReadVolume()
{
	if (handle_ == NULL)
	{
		if ( ! OpenFile(imgPath_.c_str()) )
		{
			error_ = true;
			string msg = string("DICOMImage::ReadVolume: Cannot open file: ") + imgPath_;
			cerr << endl << "\t**** " << msg << endl;
			return false;
		}
	}		

	DCM_ELEMENT element;
	if ( (DCM_GetElement(&handle_, DCM_PXLPIXELDATA, &element) != DCM_NORMAL ) || (element.length > 0xFFFFFFFE) )
	{
		error_ = true;
		string msg = string("DICOMImage::ReadVolumes: Error reading file: ") + imgPath_;
		cerr << endl << "\t**** " << msg << endl;
		return false;
	}

	// Element's length is in bytes

	volumeData_ = new U16 [element.length / 2];
	element.d.ow = volumeData_;

	U32 len;
	void *ctx = NULL;
	
	if (DCM_GetElementValue(&handle_, &element, &len, &ctx) != DCM_NORMAL)
	{
		error_ = true;
		string msg = string("DICOMImage::ReadVolume: Error reading file: ") + imgPath_;
		cerr << endl << "\t**** " << msg << endl;

		delete [] volumeData_;
		volumeData_ = NULL;
		return false;
	}

	DCM_CloseObject(&handle_);
	handle_ = NULL;

	return true;
}