// MetroProMap class implementation file
//
// purpose: A class to decode surface profile information 
//          from a Zygo white light interferometer written
//          in native binary format by Zygo's MetroPro software.
//
// author: richard.t.jones at uconn.edu
// version: November 12, 2011


#include <iostream>
#include "MetroProMap.h"
#include <math.h>

#include <TH2D.h>


int verbose = 1;

void *bigendian(char *buf);

double MetroProMap::inches2mm = 25.4;
//double MetroProMap::pixels2mm = 2.25338e-3; // nominal value, x10 objective
//double MetroProMap::pixels2mm = 1.108e-3; // calibrated for x2.0 zoom, x10 objective
//double MetroProMap::pixels2mm = 1.705e-3; // calibrated for x1.3 zoom, x10 objective
//double MetroProMap::pixels2mm = 2.217e-3; // calibrated for x1.0 zoom, x10 objective
double MetroProMap::pixels2mm = 5.490e-3; // calibrated for x0.4 zoom, x10 objective
//double MetroProMap::pixels2mm = 4.392e-3; // calibrated for x0.5 zoom, x10 objective
//double MetroProMap::pixels2mm = 0.9e-3; // calibrated for x1.0 zoom, x2.5 objective
//double MetroProMap::pixels2mm = 9.014e-3; // calibrated for x0.5 zoom, x2.5 objective

MetroProMap::MetroProMap()
{
   fHeader = 0;
   fMagicNumber = 0;
   fHeaderFormat = 0;
   fHeaderSize = 0;
   fIntens_ac_org_x = 0;
   fIntens_ac_org_y = 0;
   fIntens_ac_width = 0;
   fIntens_ac_height = 0;
   fIntens_ac_n_buckets = 0;
   fIntens_ac_range = 0;
   fIntens_ac_n_bytes = 0;
   fPhase_cn_org_x = 0;
   fPhase_cn_org_y = 0;
   fPhase_cn_width = 0;
   fPhase_cn_height = 0;
   fPhase_cn_n_bytes = 0;
   fPhase_intf_scale_factor = 0;
   fPhase_wavelength_in = 0;
   fPhase_obliquity_factor = 0;
   fPhase_pixel_pitch = 0;
   fPhase_phase_res = 0;
   fIntensData = 0;
   fPhaseData = 0;
}

MetroProMap::MetroProMap(MetroProMap &map)
{
   fMagicNumber = map.fMagicNumber;
   fHeaderFormat = map.fHeaderFormat;
   fHeaderSize = map.fHeaderSize;
   fIntens_ac_org_x = map.fIntens_ac_org_x;
   fIntens_ac_org_y = map.fIntens_ac_org_y;
   fIntens_ac_width = map.fIntens_ac_width;
   fIntens_ac_height = map.fIntens_ac_height;
   fIntens_ac_n_buckets = map.fIntens_ac_n_buckets;
   fIntens_ac_range = map.fIntens_ac_range;
   fIntens_ac_n_bytes = map.fIntens_ac_n_bytes;
   fPhase_cn_org_x = map.fPhase_cn_org_x;
   fPhase_cn_org_y = map.fPhase_cn_org_y;
   fPhase_cn_width = map.fPhase_cn_width;
   fPhase_cn_height = map.fPhase_cn_height;
   fPhase_cn_n_bytes = map.fPhase_cn_n_bytes;
   fPhase_intf_scale_factor = map.fPhase_intf_scale_factor;
   fPhase_wavelength_in = map.fPhase_wavelength_in;
   fPhase_obliquity_factor = map.fPhase_obliquity_factor;
   fPhase_pixel_pitch = map.fPhase_pixel_pitch;
   fPhase_phase_res = map.fPhase_phase_res;
   fHeader = new char[fHeaderSize];
   for (int i=0; i<fHeaderSize; ++i) {
      fHeader[i] = map.fHeader[i];
   }
   int pcount = fIntens_ac_width*fIntens_ac_height;
   fIntensData = new short int[pcount];
   for (int i=0; i<pcount; ++i) {
      fIntensData[i] = map.fIntensData[i];
   }
   pcount = fPhase_cn_width*fPhase_cn_height;
   fPhaseData = new int[pcount];
   for (int i=0; i<pcount; ++i) {
      fPhaseData[i] = map.fPhaseData[i];
   }
}

MetroProMap::~MetroProMap()
{
   if (fHeader) {
      delete [] fHeader;
      fHeader = 0;
   }
   if (fIntensData) {
      delete [] fIntensData;
      fIntensData = 0;
   }
   if (fPhaseData) {
      delete [] fPhaseData;
      fPhaseData = 0;
   }
}

int MetroProMap::getRfactor() const
{
   if (fHeaderFormat == 1 && fPhase_phase_res == 0) {
      return 4096;
   }
   else if (fHeaderFormat == 1 && fPhase_phase_res == 1) {
      return 32768;
   }
   else if (fHeaderFormat > 1 && fPhase_phase_res == 0) {
      return 4096;
   }
   else if (fHeaderFormat > 1 && fPhase_phase_res == 1) {
      return 32768;
   }
   else if (fHeaderFormat > 1 && fPhase_phase_res == 2) {
      return 131072;
   }
   else {
      std::cerr << "Error in MetroProMap::getRfactor()"
                << " - unrecognized value for phase_res parameter"
                << std::endl;
      return -1;
   }
}

double MetroProMap::getAverageHeight(int inmm) const
{
   if (fPhaseData == 0) {
      return 0;
   }
   int sum0=0;
   double sum1=0;
   int npixels=fPhase_cn_width*fPhase_cn_height;
   for (int i=0; i < npixels; ++i) {
std::cout << i << " : " << fPhaseData[i] << std::endl;
      if (fPhaseData[i] < 2147483640) {
         sum1+= fPhaseData[i];
         sum0++;
      }
   }
   return sum1/sum0*getZygos2height(inmm);
}

double MetroProMap::getMaximumHeight(int inmm) const
{
   if (fPhaseData == 0) {
      return 0;
   }
   int zmax=0;
   int npixels=fPhase_cn_width*fPhase_cn_height;
   for (int i=0; i < npixels; ++i) {
      if (fPhaseData[i] < 2147483640) {
         zmax = (zmax > fPhaseData[i])? zmax : fPhaseData[i];
      }
   }
   return zmax*getZygos2height(inmm);
}

double MetroProMap::getMinimumHeight(int inmm) const
{
   if (fPhaseData == 0) {
      return 0;
   }
   int zmin=0x4fffffff;
   int npixels=fPhase_cn_width*fPhase_cn_height;
   for (int i=0; i < npixels; ++i) {
      if (fPhaseData[i] < 2147483640) {
         zmin = (zmin < fPhaseData[i])? zmin : fPhaseData[i];
      }
   }
   return zmin*getZygos2height(inmm);
}

double MetroProMap::getZygos2height(int inmm) const
{
   return fPhase_wavelength_in
         *fPhase_intf_scale_factor
         *fPhase_obliquity_factor
         /getRfactor()*1000
         /((inmm)? 1 : inches2mm);
}

int MetroProMap::getIntensity(int ix, int iy) const
{
   ix -= fIntens_ac_org_x;
   iy -= fIntens_ac_org_y;
   if (ix < 0 || ix >= fIntens_ac_width ||
       iy < 0 || iy >= fIntens_ac_height) {
      return 0;
   }
   else {
      return fIntensData[ix+iy*fIntens_ac_width];
   }
}

int MetroProMap::getPhase(int ix, int iy) const
{
   ix -= fPhase_cn_org_x;
   iy -= fPhase_cn_org_y;
   if (ix < 0 || ix >= fPhase_cn_width ||
       iy < 0 || iy >= fPhase_cn_height) {
      return 0;
   }
   else {
      return fPhaseData[ix+iy*fPhase_cn_width];
   }
}

double MetroProMap::getElevation(int ix, int iy, int inmm) const
{
   double phase = getPhase(ix,iy);
   return (phase < 2147483640)?  phase*getZygos2height(inmm): 0;
}

std::istream &operator>>(std::istream &istr, MetroProMap &map)
{
   char buf[256];
   union packed_t {
      char Char[4];
      short int Short[2];
      int Int;
      float Float;
   } *packed;

   istr.read(buf,256);

   packed = (union packed_t*)bigendian(&buf[0]);
   map.fMagicNumber = packed->Int;
   packed = (union packed_t*)bigendian(&buf[4]);
   map.fHeaderFormat = packed->Short[1];
   packed = (union packed_t*)bigendian(&buf[6]);
   map.fHeaderSize = packed->Int;
   if (map.fMagicNumber == (int)0x881b036f &&
       map.fHeaderFormat == 1 && map.fHeaderSize == 834) {
      if (map.fHeader) {
         delete [] map.fHeader;
      }
      map.fHeader = new char[834];
   }
   else if (map.fMagicNumber == (int)0x881b0370 &&
            map.fHeaderFormat == 2 && map.fHeaderSize == 834) {
      if (map.fHeader) {
         delete [] map.fHeader;
      }
      map.fHeader = new char[834];
   }
   else if (map.fMagicNumber == (int)0x881b0371 &&
            map.fHeaderFormat == 3 && map.fHeaderSize == 4096) {
      if (map.fHeader) {
         delete [] map.fHeader;
      }
      map.fHeader = new char[4096];
   }
   else {
      std::cerr << "Error in std::istream &operator>>(istream&,MetroProMap&"
                << " - invalid header read from input stream." << std::endl
                << "   magic number = " << std::hex
                << map.fMagicNumber << std::endl
                << "   header_format = " << std::dec
                << map.fHeaderFormat << std::endl
                << "   header_size = " << std::dec
                << map.fHeaderSize << std::endl;
      return istr;
   }

   for (int i=0; i < 256; ++i) {
      map.fHeader[i] = buf[i];
   }
   istr.read(&map.fHeader[256],map.fHeaderSize-256);
   if (istr.fail()) {
      std::cerr << "Error in std::istream &operator>>(istream&,MetroProMap&)"
                << " - i/o error reading header from input stream."
                << std::endl;
      return istr;
   }
   else if (istr.eof()) {
      std::cerr << "Error in std::istream &operator>>(istream&,MetroProMap&)"
                << " - end of file while reading header from input stream."
                << std::endl;
      return istr;
   }

#ifdef DEBUG_LISTING_OF_HEADER
   for (int i=0; i<map.fHeaderSize; i+=4) {
      packed = (union packed_t*)bigendian(&buf[i]);
      std::cerr << i << " : " << packed->Float << std::endl;
   }
#endif

   packed = (union packed_t*)bigendian(&buf[48]);
   map.fIntens_ac_org_x = packed->Short[1];
   map.fIntens_ac_org_y = packed->Short[0];
   packed = (union packed_t*)bigendian(&buf[52]);
   map.fIntens_ac_width = packed->Short[1];
   map.fIntens_ac_height = packed->Short[0];
   packed = (union packed_t*)bigendian(&buf[56]);
   map.fIntens_ac_n_buckets = packed->Short[1];
   map.fIntens_ac_range = packed->Short[0];
   packed = (union packed_t*)bigendian(&buf[60]);
   map.fIntens_ac_n_bytes = packed->Int;
   packed = (union packed_t*)bigendian(&buf[64]);
   map.fPhase_cn_org_x = packed->Short[1];
   map.fPhase_cn_org_y = packed->Short[0];
   packed = (union packed_t*)bigendian(&buf[68]);
   map.fPhase_cn_width = packed->Short[1];
   map.fPhase_cn_height = packed->Short[0];
   packed = (union packed_t*)bigendian(&buf[72]);
   map.fPhase_cn_n_bytes = packed->Int;
   packed = (union packed_t*)bigendian(&buf[164]);
   map.fPhase_intf_scale_factor = packed->Float;
   packed = (union packed_t*)bigendian(&buf[168]);
   map.fPhase_wavelength_in = packed->Float;
   packed = (union packed_t*)bigendian(&buf[176]);
   map.fPhase_obliquity_factor = packed->Float;
   packed = (union packed_t*)bigendian(&buf[184]);
   map.fPhase_pixel_pitch = packed->Float;
   packed = (union packed_t*)bigendian(&buf[216]);
   map.fPhase_phase_res = packed->Short[0];

   if (map.fIntensData) {
      delete [] map.fIntensData;
      map.fIntensData = 0;
   }
   if (map.fIntens_ac_n_bytes) {
      map.fIntensData = new short int[map.fIntens_ac_n_bytes/sizeof(short int)];
      istr.read((char*)map.fIntensData,map.fIntens_ac_n_bytes);
      int npixels=map.fIntens_ac_width*map.fIntens_ac_height;
      for (int i=0; i < npixels; i+=2) {
         packed = (union packed_t*)bigendian((char*)&map.fIntensData[i]);
         map.fIntensData[i] = packed->Short[1];
         map.fIntensData[i+1] = packed->Short[0];
      }
   }

   if (map.fPhaseData) {
      delete [] map.fPhaseData;
      map.fPhaseData = 0;
   }
   if (map.fPhase_cn_n_bytes) {
      map.fPhaseData = new int[map.fPhase_cn_n_bytes/sizeof(int)];
      istr.read((char*)map.fPhaseData,map.fPhase_cn_n_bytes);
      int npixels=map.fPhase_cn_width*map.fPhase_cn_height;
      for (int i=0; i < npixels; ++i) {
         packed = (union packed_t*)bigendian((char*)&map.fPhaseData[i]);
         map.fPhaseData[i] = packed->Int;
      }
   }
   if (istr.fail()) {
      std::cerr << "Error in std::istream &operator>>(istream&,MetroProMap&)"
                << " - i/o error reading header from input stream."
                << std::endl;
      return istr;
   }
   else if (istr.eof()) {
      std::cerr << "Error in std::istream &operator>>(istream&,MetroProMap&)"
                << " - end of file while reading header from input stream."
                << std::endl;
      return istr;
   }

   if (verbose) {
      std::cout << "istream &operator>>(istream&,MetroProMap) saw some data:"
                << std::endl
                << "   magic_number: " << std::hex
                << map.fMagicNumber << std::endl
                << "   header_format: " << std::dec
                << map.fHeaderFormat << std::endl
                << "   header_size: " << std::dec
                << map.fHeaderSize << std::endl
                << "   intensity section: " << std::endl
                << "      ac_org_x: " << map.fIntens_ac_org_x << std::endl
                << "      ac_org_y: " << map.fIntens_ac_org_y << std::endl
                << "      ac_width: " << map.fIntens_ac_width << std::endl
                << "      ac_height: " << map.fIntens_ac_height << std::endl
                << "      ac_n_buckets: " << map.fIntens_ac_n_buckets << std::endl
                << "      ac_range: " << map.fIntens_ac_range << std::endl
                << "      ac_n_bytes: " << map.fIntens_ac_n_bytes << std::endl
                << "   phase section: " << std::endl
                << "      cn_org_x: " << map.fPhase_cn_org_x << std::endl
                << "      cn_org_y: " << map.fPhase_cn_org_y << std::endl
                << "      cn_width: " << map.fPhase_cn_width << std::endl
                << "      cn_height: " << map.fPhase_cn_height << std::endl
                << "      cn_n_bytes: " << map.fPhase_cn_n_bytes << std::endl
                << "      intf_scale_factor: " 
                << map.fPhase_intf_scale_factor << std::endl
                << "      wavelength_in: " 
                << map.fPhase_wavelength_in << std::endl
                << "      obliquity_factor: " 
                << map.fPhase_obliquity_factor << std::endl
                << "      pixel_pitch: " 
                << map.fPhase_pixel_pitch << std::endl
                << "      phase_res: " << map.fPhase_phase_res << std::endl;
   }
   return istr;
}

char *bigendian_buf = 0;

void *bigendian(char *buf)
{
   if (bigendian_buf == 0) {
      bigendian_buf = new char[4];
   }
   bigendian_buf[3] = buf[0];  
   bigendian_buf[2] = buf[1];  
   bigendian_buf[1] = buf[2];  
   bigendian_buf[0] = buf[3];  
   return (void*)bigendian_buf;
}

TH2D *MetroProMap::makeHist2D(std::string name, std::string title)
{
   int nxbins = getWidth();
   int nybins = getHeight();
   double pitch = getPixelPitch(1);
   if (pitch > 0) {
      pixels2mm = pitch;
   }
   double xlo = 0;
   double xhi = nxbins*pixels2mm;
   double ylo = 0;
   double yhi = nybins*pixels2mm;
   TH2D *h2 = new TH2D(name.c_str(),title.c_str(),
                       nxbins,xlo,xhi,nybins,ylo,yhi);
   for (int iy=0; iy < nybins; ++iy) {
      for (int ix=0; ix < nxbins; ++ix) {
         h2->SetBinContent(ix+1,iy+1,getElevation(ix,iy));
      }
   }
   return h2;
}