//
// ablator: ROOT c++ class implementing the capability to analyze a surface
//          profile image of a diamond sample (eg. from a Zygo microscope)
//          and create a laser ablation sequence file that can be used to
//          mill that sample down to match an arbitrary target surface.
//
// author: richard.t.jones at uconn.edu
// version: february 10, 2015
//
// usage example:
//   $ root -l
//   root [0] .L Map2D.cc+O
//   root [1] .L ablator.C+O
//   root [2] TFile fin("ablator.root");
//   root [3] ablator abl();
//   root [4] Map2D *virm = abl.mill(virgin, target_smoothed);
//   root [5] virm->Draw("colz");
//

#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <string>
#include <list>

#include <TROOT.h>
#include <TRandom.h>
#include <TFFTComplexReal.h>
#include <TDecompSVD.h>

#include "ablator.h"

#define SMOOTH_REVERSIBILITY_CHECK 1
#define NAIVE_GAUSSIAN_FOCAL_SPOT 0

double ablator::get_cut_depth() const
{
   // Imagine an infinite plane populated with single laser pulses laid out 
   // in a regular rectangular array with pitch fRaster_dx, fRaster_dy. Compute
   // the average depth of the cut surface below the original surface once this
   // array of pulses has been deposited.

   int nx = fSpotHsize / fSpotHresol;
   int ny = fSpotVsize / fSpotVresol;
   Map2D templ(TH2D("pvt:spot_map", "private spot map",
                    nx, 0, fSpotHsize, ny, 0, fSpotVsize));
   Map2D *spot = make_beam_spot(&templ);
   double integral = spot->Integral();
   delete spot;
   return integral * (fSpotHresol * fSpotVresol)
                   / (fRaster_dx * fRaster_dy);
}

double ablator::get_max_pit_depth() const
{
   // Pulse-model-independent method for requesting the maximum depth (microns)
   // of the pit created in a virgin surface by a single laser pulse.

#if NAIVE_GAUSSIAN_FOCAL_SPOT
   return getSpotDepth();
#else
   return getSpotPitDepth();
#endif
}

Map2D *ablator::make_beam_spot(const Map2D *templ) const
{
   // Build a map in two dimensions of the pit depth profile generated
   // by a single laser pulse centered at the origin. The origin of the
   // spot is in the four corners of the map, so that the four quadrants
   // of the pit profile are located in the four corners. The output map
   // the same shape as input argument templ, with the x,y coordinates
   // assumed to be in mm and the z coordinate assumed to be in microns.

   // Allow the user to inject a custom beam spot
   // simply by having a map called "custom_beam_spot" 
   // defined within the current ROOT environment.

   Map2D *custom = (Map2D*)gROOT->FindObject("custom_beam_spot");
   if (custom != 0) {
      Map2D *spot = (Map2D*)custom->Clone("beam_spot");
      spot->SetTitle("single-pulse spot depth profile");
      return spot;
   }

   Map2D *spot = (Map2D*)templ->Clone("beam_spot");
   spot->SetTitle("single-pulse spot depth profile");

   int ndx = spot->GetNbinsX();
   int ndy = spot->GetNbinsY();
   for (int jx = 0; jx < ndx / 2; ++jx) {
      for (int jy = 0; jy < ndy / 2; ++jy) {

#if NAIVE_GAUSSIAN_FOCAL_SPOT
         double du = spot->GetXaxis()->GetBinCenter(jx+1) / fSpotHrms;
         double dv = spot->GetYaxis()->GetBinCenter(jy+1) / fSpotVrms;
         double depth = fSpotDepth * exp(-0.5 * (du*du + dv*dv));
         spot->SetBinContent(jx+1, jy+1, depth);
         spot->SetBinContent(ndx-jx, jy+1, depth);
         spot->SetBinContent(jx+1, ndy-jy, depth);
         spot->SetBinContent(ndx-jx, ndy-jy, depth);
#else
         double u = spot->GetXaxis()->GetBinCenter(jx) / fSpotHsigma;
         double xfactor = (fSpotPitDepth + fSpotHbase) * exp(-0.5 * u*u)
                                         - fSpotHbase;
         xfactor = (xfactor > 0)? xfactor : 0;
         double y = spot->GetYaxis()->GetBinCenter(jy);
         double v0 = y / fSpotVsigma0;
         double v1 = y / fSpotVsigma1;
         double yfactor = fSpotVampli0 * exp(-0.5 * v0*v0)
                        + fSpotVampli1 * exp(-0.5 * v1*v1);
         yfactor = 1.01 * yfactor - 0.01 * fSpotPitDepth;
         yfactor = (yfactor > 0)? yfactor : 0;
         double v2 = (y - fSpotVoffset2) / fSpotVsigma2;
         double yfactor2 = (fSpotVampli2 + fSpotVbase2) * exp(-0.5 * v2*v2)
                                         - fSpotVbase2;
         yfactor2 = (yfactor2 > 0)? yfactor2 : 0;
         yfactor += yfactor2;
         double depth = xfactor * yfactor / fSpotPitDepth;
         spot->SetBinContent(jx+1, jy+1, depth);
         spot->SetBinContent(ndx-jx, jy+1, depth);
         spot->SetBinContent(jx+1, ndy-jy, depth);
         spot->SetBinContent(ndx-jx, ndy-jy, depth);
#endif

      }
   }
   return spot;
}

Map2D *ablator::make_gaussian_spot(const Map2D *templ,
                                   double sigx, double sigy) const
{
   // Build a 2D gaussian smearing function with sigmas sigx,sigy.
   // The origin of the spot is in the four corners of the map, so that
   // the four quadrants of the pit profile are located in the four
   // corners. The output map has the same shape as input argument
   // templ, with the x,y coordinates assumed to be in mm and the z
   // coordinate adjusted so that the sum over all bins is unity.

   Map2D *spot = (Map2D*)templ->Clone("gaussian_spot");
   int ndx = templ->GetNbinsX();
   int ndy = templ->GetNbinsY();
   for (int jx = 0; jx < ndx / 2; ++jx) {
      for (int jy = 0; jy < ndy / 2; ++jy) {
         double du = spot->GetXaxis()->GetBinCenter(jx+1) / sigx;
         double dv = spot->GetYaxis()->GetBinCenter(jy+1) / sigy;
         double depth = fSpotDepth * exp(-0.5 * (du*du + dv*dv));
         spot->SetBinContent(jx+1, jy+1, depth);
         spot->SetBinContent(ndx-jx, jy+1, depth);
         spot->SetBinContent(jx+1, ndy-jy, depth);
         spot->SetBinContent(ndx-jx, ndy-jy, depth);
      }
   }
   spot->Rescale(1, 1, 1/spot->Integral());
   return spot;
}

Map2D *ablator::mill(const Map2D *part_map, Map2D *target_map)
{
   // Generate an ablation pulse sequence that mills a surface with an
   // initial profile given by part_map down as close as possible to the
   // profile specified in target_map.  Only regions of target_map which
   // are lower (smaller height values) than the corresponding region in
   // part_map will be touched during the ablation sequence, so if there
   // are areas of the surface that should not be touched, they must be
   // assigned a large value in target_map to protect them. A new Map2D
   // with the exterior dimensions of part_map and a bin size set by the
   // pitch of the laser raster is created to hold the number of laser
   // pulses for each point in the raster that is needed to bring the
   // part down to the level of the target. A new map is generated to
   // contain an estimate of the final shape expected for the ablated
   // surface, and a pointer to this map is returned at exit.  It is up
   // to the user to destroy this object when it is no longer needed.

   std::cout << "ablator::mill - per-pass cut depth is " 
             << get_cut_depth()
             << std::endl;

   // form the shape of the beam spot and take its DFT

   Map2D *spot = make_beam_spot(part_map);
   spot->Rescale(1, 1, 1/spot->Integral());
   Map2D *spotFFTr = (Map2D*)spot->Clone("spotFFTr");
   spot->FFT(spotFFTr,"RE R2C P");
   Map2D *spotFFTi = (Map2D*)spot->Clone("spotFFTi");
   spot->FFT(spotFFTi,"IM R2C P");

   // find the shape of the target surface and take its DFT

   int nxbins = part_map->GetNbinsX();
   int nybins = part_map->GetNbinsY();
   Map2D *shape = (Map2D*)part_map->Clone("shape");
   for (int ix = 1; ix <= nxbins; ++ix) { 
      for (int iy = 1; iy <= nybins; ++iy) { 
         double x = shape->GetXaxis()->GetBinCenter(ix);
         double y = shape->GetYaxis()->GetBinCenter(iy);
         double tar = target_map->Interpolate(x,y);
         shape->SetBinContent(ix, iy, tar);
      }
   }
   Map2D *shapeFFTr = (Map2D*)shape->Clone("shapeFFTr");
   shape->FFT(shapeFFTr,"RE R2C P");
   Map2D *shapeFFTi = (Map2D*)shape->Clone("shapeFFTi");
   shape->FFT(shapeFFTi,"IM R2C P");

   // take the smoothing out of the target map
 
   int dim[2] = {nxbins, nybins};
   TVirtualFFT *invFFT = TVirtualFFT::FFT(2, dim, "C2R");
   for (int ix = 1; ix <= nxbins; ++ix) {
      for (int iy = 1; iy <= nybins/2 + 1; ++iy) {
         double spot_r = spotFFTr->GetBinContent(ix,iy);
         double spot_i = spotFFTi->GetBinContent(ix,iy);
         double spot_mag2 = spot_r*spot_r + spot_i*spot_i;
         double shape_r = shapeFFTr->GetBinContent(ix,iy);
         double shape_i = shapeFFTi->GetBinContent(ix,iy);
         double ratio_r = (shape_r*spot_r + shape_i*spot_i)/spot_mag2;
         double ratio_i = (shape_i*spot_r - shape_r*spot_i)/spot_mag2;
         int index[2] = {ix - 1, iy - 1};
         if (spot_mag2 > 1e-20) {
            invFFT->SetPoint(index, ratio_r, ratio_i);
         }
         else {
            invFFT->SetPoint(index, 0, 0);
         }
      }
   }
   invFFT->Transform();
   Map2D *sharp = (Map2D*)shape->Clone("sharp");
   sharp = (Map2D*)TH2::TransformHisto(invFFT, sharp, "Re");
   sharp->Rescale(1, 1, 1.0/(nxbins*nybins));
   sharp->SetTitle("sharpened target image");

   // take the smoothing out of the part map

   Map2D *part_sharp = unsmooth(part_map);

   // now take the difference between the part and sharpened target maps
 
   for (int ix = 1; ix <= nxbins; ++ix) { 
      for (int iy = 1; iy <= nybins; ++iy) { 
         double x = sharp->GetXaxis()->GetBinCenter(ix);
         double y = sharp->GetYaxis()->GetBinCenter(iy);
         double par = part_sharp->Interpolate(x,y);
         double tar = sharp->Interpolate(x,y);
         shape->SetBinContent(ix, iy, (tar <= par+1e100)? par - tar : 0);
      }
   }

   // reconvolve the beam spot with the solved-for program shape

   Map2D *programFFTr = (Map2D*)shape->Clone("programFFTr");
   shape->FFT(programFFTr,"RE R2C P");
   Map2D *programFFTi = (Map2D*)shape->Clone("programFFTi");
   shape->FFT(programFFTi,"IM R2C P");
   invFFT = TVirtualFFT::FFT(2, dim, "C2R");
   for (int ix = 1; ix <= nxbins; ++ix) {
      for (int iy = 1; iy <= nybins/2 + 1; ++iy) {
         double spot_r = spotFFTr->GetBinContent(ix,iy);
         double spot_i = spotFFTi->GetBinContent(ix,iy);
         double prog_r = programFFTr->GetBinContent(ix,iy);
         double prog_i = programFFTi->GetBinContent(ix,iy);
         double prod_r = (prog_r*spot_r - prog_i*spot_i);
         double prod_i = (prog_i*spot_r + prog_r*spot_i);
         int index[2] = {ix - 1, iy - 1};
         invFFT->SetPoint(index, prod_r, prod_i);
      }
   }
   invFFT->Transform();
   Map2D *removed = (Map2D*)shape->Clone("removed");
   removed = (Map2D*)TH2::TransformHisto(invFFT, removed, "Re");
   removed->Rescale(1, 1, 1.0/(nxbins*nybins));
   removed->SetTitle("removed material");

   delete programFFTr;
   delete programFFTi;
   delete shapeFFTr;
   delete shapeFFTi;

   // generate a map of the milled part

   std::string name(part_map->GetName());
   name += "_milled";
   Map2D *result = (Map2D*)shape->Clone(name.c_str());
   result->ZeroSuppress(1e100);
   result->Fill(part_sharp);
   result->Add(shape, -1);
 
   // reconvolve the beam shape with the milled part result

   Map2D *resultFFTr = (Map2D*)result->Clone("resultFFTr");
   result->FFT(resultFFTr,"RE R2C P");
   Map2D *resultFFTi = (Map2D*)result->Clone("resultFFTi");
   result->FFT(resultFFTi,"IM R2C P");
   invFFT = TVirtualFFT::FFT(2, dim, "C2R");
   for (int ix = 1; ix <= nxbins; ++ix) {
      for (int iy = 1; iy <= nybins/2 + 1; ++iy) {
         double spot_r = spotFFTr->GetBinContent(ix,iy);
         double spot_i = spotFFTi->GetBinContent(ix,iy);
         double resu_r = resultFFTr->GetBinContent(ix,iy);
         double resu_i = resultFFTi->GetBinContent(ix,iy);
         double prod_r = (resu_r*spot_r - resu_i*spot_i);
         double prod_i = (resu_i*spot_r + resu_r*spot_i);
         int index[2] = {ix - 1, iy - 1};
         invFFT->SetPoint(index, prod_r, prod_i);
      }
   }
   invFFT->Transform();
   result = (Map2D*)TH2::TransformHisto(invFFT, result, "Re");
   result->Rescale(1, 1, 1.0/(nxbins*nybins));
   delete resultFFTr;
   delete resultFFTi;

   // rescale program to units of laser pulses on the raster grid

   double xmin = part_map->GetXaxis()->GetXmin();
   double xmax = part_map->GetXaxis()->GetXmax();
   double ymin = part_map->GetYaxis()->GetXmin();
   double ymax = part_map->GetYaxis()->GetXmax();
   int nxraster = (int)floor((xmax - xmin)/fRaster_dx + 0.5);
   int nyraster = (int)floor((ymax - ymin)/fRaster_dy + 0.5);
   xmax = xmin + fRaster_dx * nxraster;
   ymax = ymin + fRaster_dy * nyraster;
   Map2D *program = new Map2D(TH2D("program", "milling program",
                                   nxraster, xmin, xmax, nyraster, ymin, ymax));
   program->Fill(shape);
   program->Rescale(1, 1, 1/get_cut_depth());

   delete spot;
   delete spotFFTr;
   delete spotFFTi;
   delete invFFT;

   return result;
}

Map2D *ablator::smooth(Map2D *target, double smear_factor)
{
   // Tool to create a realistic target surface profile, starting from an
   // idealized profile with perfectly sharp edges, specified by the target
   // argument. This is accomplished by convoluting the target map with a
   // 2D Gaussian smearing function. A new map is generated containing
   // smoothed surface profile, and passed back by pointer in the return
   // value. It is up to the user code to destroy this object when it is
   // no longer needed. The smear_factor argument is used to rescale the
   // naive gaussian beam spot sigmas before taking the convolution with
   // the target map. If smear_factor==0 (default) then the actual beam
   // spot profile is used instead of a naive gaussian.

   // form the shape of the beam spot and take its DFT

   Map2D *fatspot;
   if (smear_factor > 0) {
      fatspot = make_gaussian_spot(target, smear_factor * fXedgeParam, 
                                           smear_factor * fYedgeParam);
   }
   else {
      fatspot = make_beam_spot(target);
      fatspot->Rescale(1, 1, 1/fatspot->Integral());
   }
   fatspot->SetName("fatspot");
   fatspot->SetTitle("smoothing spot");
   Map2D *fatspotFFTr = (Map2D*)fatspot->Clone("fatspotFFTr");
   fatspot->FFT(fatspotFFTr,"RE R2C P");
   Map2D *fatspotFFTi = (Map2D*)fatspot->Clone("fatspotFFTi");
   fatspot->FFT(fatspotFFTi,"IM R2C P");

   // find the shape of the target surface and take its DFT

   Map2D *targetFFTr = (Map2D*)target->Clone("targetFFTr");
   target->FFT(targetFFTr,"RE R2C P");
   Map2D *targetFFTi = (Map2D*)target->Clone("targetFFTi");
   target->FFT(targetFFTi,"IM R2C P");

   // use the convolution theorem to fold the target and beam spot shapes

   int nxbins = target->GetNbinsX();
   int nybins = target->GetNbinsY();
   int dim[2] = {nxbins, nybins};
   TVirtualFFT *invFFT = TVirtualFFT::FFT(2, dim, "C2R");
   for (int ix = 1; ix <= nxbins; ++ix) {
      for (int iy = 1; iy <= nybins / 2 + 1; ++iy) {
         double targ_r = targetFFTr->GetBinContent(ix,iy);
         double targ_i = targetFFTi->GetBinContent(ix,iy);
         double spot_r = fatspotFFTr->GetBinContent(ix,iy);
         double spot_i = fatspotFFTi->GetBinContent(ix,iy);
         double prod_r = targ_r*spot_r - targ_i*spot_i;
         double prod_i = targ_r*spot_i + targ_i*spot_r;
         int index[2] = {ix - 1, iy - 1};
         invFFT->SetPoint(index, prod_r, prod_i);
      }
   }
   invFFT->Transform();
   std::string name(target->GetName());
   name += "_smoothed";
   Map2D *result = (Map2D*)target->Clone(name.c_str());
   result = (Map2D*)TH2::TransformHisto(invFFT, result, "Re");
   result->Rescale(1, 1, 1.0/(nxbins*nybins));

#if SMOOTH_REVERSIBILITY_CHECK

   // run the convolution backwards to check reversibility
 
   Map2D *resultFFTr = (Map2D*)result->Clone("resultFFTr");
   result->FFT(resultFFTr,"RE R2C P");
   Map2D *resultFFTi = (Map2D*)result->Clone("resultFFTi");
   result->FFT(resultFFTi,"IM R2C P");

   invFFT = TVirtualFFT::FFT(2, dim, "C2R");
   for (int ix = 1; ix <= nxbins; ++ix) {
      for (int iy = 1; iy <= nybins / 2 + 1; ++iy) {
         double targ_r = resultFFTr->GetBinContent(ix,iy);
         double targ_i = resultFFTi->GetBinContent(ix,iy);
         double targ_mag2 = targ_r*targ_r + targ_i*targ_i;
         double spot_r = fatspotFFTr->GetBinContent(ix,iy);
         double spot_i = fatspotFFTi->GetBinContent(ix,iy);
         double spot_mag2 = spot_r*spot_r + spot_i*spot_i;
         double prod_r = (targ_r*spot_r + targ_i*spot_i)/spot_mag2;
         double prod_i = (targ_i*spot_r - targ_r*spot_i)/spot_mag2;
         int index[2] = {ix - 1, iy - 1};
         if (spot_mag2 > 1e-24) {
            invFFT->SetPoint(index, prod_r, prod_i);
         }
         else {
            invFFT->SetPoint(index, 0, 0);
         }
      }
   }
   invFFT->Transform();
   name = target->GetName();
   name += "_unsmoothed";
   Map2D *unsmoothed = (Map2D*)target->Clone(name.c_str());
   unsmoothed = (Map2D*)TH2::TransformHisto(invFFT, unsmoothed, "Re");
   unsmoothed->Rescale(1, 1, 1.0/(nxbins*nybins));

   delete resultFFTr;
   delete resultFFTi;

#endif

   delete fatspot;
   delete fatspotFFTr;
   delete fatspotFFTi;
   delete targetFFTr;
   delete targetFFTi;
   delete invFFT;

   return result;
}

Map2D *ablator::unsmooth(const Map2D *part, double noise_cutoff)
{
   // Tool to take a part image and deconvolve it with the beam spot profile.
   // This function is the numerical inverse of smooth() so it should give
   // well-behaved results for maps that were created by smooth() acting on
   // a source map, giving back the source map with minimal residual errors.
   // The errors are controlled by the noise_cutoff argument, which allows
   // the user to restrict the noise coming from high-frequency modes with
   // very small amplitude that get amplified in the deconvolution. 

   // form the shape of the beam spot and take its DFT

   Map2D *fatspot;
   fatspot = make_beam_spot(part);
   fatspot->Rescale(1, 1, 1/fatspot->Integral());
   fatspot->SetName("fatspot");
   fatspot->SetTitle("smoothing spot");
   Map2D *fatspotFFTr = (Map2D*)fatspot->Clone("fatspotFFTr");
   fatspot->FFT(fatspotFFTr,"RE R2C P");
   Map2D *fatspotFFTi = (Map2D*)fatspot->Clone("fatspotFFTi");
   fatspot->FFT(fatspotFFTi,"IM R2C P");

   // find the shape of the part surface and take its DFT

   Map2D *result = (Map2D*)part->Clone("result");
   Map2D *partFFTr = (Map2D*)part->Clone("partFFTr");
   result->FFT(partFFTr,"RE R2C P");
   Map2D *partFFTi = (Map2D*)part->Clone("partFFTi");
   result->FFT(partFFTi,"IM R2C P");

   // run the deconvolution
 
   double amp2_cutoff = 1e-24;
   int nxbins = part->GetNbinsX();
   int nybins = part->GetNbinsY();
   int dim[2] = {nxbins, nybins};
   TVirtualFFT *invFFT = TVirtualFFT::FFT(2, dim, "C2R");
   for (int ix = 1; ix <= nxbins; ++ix) {
      for (int iy = 1; iy <= nybins / 2 + 1; ++iy) {
         double part_r = partFFTr->GetBinContent(ix,iy);
         double part_i = partFFTi->GetBinContent(ix,iy);
         double part_mag2 = part_r*part_r + part_i*part_i;
         double spot_r = fatspotFFTr->GetBinContent(ix,iy);
         double spot_i = fatspotFFTi->GetBinContent(ix,iy);
         double spot_mag2 = spot_r*spot_r + spot_i*spot_i;
         double ratio_r = (part_r*spot_r + part_i*spot_i)/spot_mag2;
         double ratio_i = (part_i*spot_r - part_r*spot_i)/spot_mag2;
         int index[2] = {ix - 1, iy - 1};
         if (spot_mag2 > amp2_cutoff &&
             part_mag2 * noise_cutoff > nxbins * nybins)
         {
            invFFT->SetPoint(index, ratio_r, ratio_i);
         }
         else {
            invFFT->SetPoint(index, 0, 0);
         }
      }
   }
   invFFT->Transform();
   std::string name(part->GetName());
   name += "_unsmoothed";
   Map2D *unsmoothed = (Map2D*)part->Clone(name.c_str());
   unsmoothed = (Map2D*)TH2::TransformHisto(invFFT, unsmoothed, "Re");
   unsmoothed->Rescale(1, 1, 1.0/(nxbins*nybins));

   delete fatspot;
   delete fatspotFFTr;
   delete fatspotFFTi;
   delete partFFTr;
   delete partFFTi;
   delete invFFT;

   return unsmoothed;
}

Map2D *ablator::unsmooth2(const Map2D *part, double eigenvalue_cutoff)
{
   // Tool to take a part image and deconvolve it with the beam spot profile.
   // The exact algorithm would greatly amplify any noise present in the part
   // image, so a refined procedure based on singular value decomposition is
   // used to invert the smearing matrix and only eigenvalues above a certain
   // minimum cutoff are included when computing the inverse convolution.

   Map2D *spot = make_beam_spot(part);
   spot->Rescale(1, 1, 1 / spot->Integral());
   int nxbins = part->GetNbinsX();
   int nybins = part->GetNbinsY();
   TMatrixD P(nxbins,nybins);
   for (int ix=0; ix < nxbins; ++ix)
      for (int iy=0; iy < nybins; ++iy)
         P(ix,iy) = part->GetBinContent(ix+1,iy+1);

   // Assume that the convolution function is factorized
   // like F(x,y) = f(x) g(y), for example a Gaussian with
   // major, minor axes aligned with the x,y axes.

   TMatrixD xspot(nxbins,nxbins);
   for (int row = 0; row < nxbins; ++row) {
      for (int col = 0; col < nxbins; ++col) {
         int colrel = (nxbins + row - col) % nxbins;
         xspot(row,col) = spot->GetBinContent(1+colrel,1);
      }
   }
   TDecompSVD xsvd(xspot);
   if (! xsvd.Decompose()) {
      std::cerr << "ablator::unsmooth2 - Error in TDecompSVD, matrix"
                << " decomposition failed, cannot continue!"
                << std::endl;
      return 0;
   }
   const TVectorD &xsig = xsvd.GetSig();
#if 0
   std::cout << "xsvd eigenvalues are listed below: " << std::endl;
   for (int eig=0; eig < xsig.GetNrows(); ++eig) {
      std::cout << eig << " : " << xsig[eig] << std::endl;
   }
#endif
   const TMatrixD &xU = xsvd.GetU();
   check_orthogonality(xU);
   const TMatrixD &xV = xsvd.GetV();
   check_orthogonality(xV);
   TMatrixD xfilter(nxbins,nxbins);
   xfilter *= 0;
   for (int eig = 0; eig < xsvd.GetNrows(); ++eig) {
      if (xsig[eig] / xsig[0] < eigenvalue_cutoff)
         break;
      double xsiginv = 1 / xsig[eig];
      for (int row = 0; row < nxbins; ++row)
         for (int col = 0; col < nxbins; ++col)
            xfilter(row,col) += xU(row,eig) * xsiginv * xV(col,eig);
   }

   TMatrixD yspot(nybins,nybins);
   for (int row = 0; row < nybins; ++row) {
      for (int col = 0; col < nybins; ++col) {
         int colrel = (nybins + row - col) % nybins;
         yspot(row,col) = spot->GetBinContent(1,1+colrel);
      }
   }
   TDecompSVD ysvd(yspot);
   if (! ysvd.Decompose()) {
      std::cerr << "ablator::unsmooth2 - Error in TDecompSVD, matrix"
                << " decomposition failed, cannot continue!"
                << std::endl;
      return 0;
   }
   const TVectorD &ysig = ysvd.GetSig();
#if 0
   std::cout << "ysvd eigenvalues are listed below: " << std::endl;
   for (int eig=0; eig < ysig.GetNrows(); ++eig) {
      std::cout << eig << " : " << ysig[eig] << std::endl;
   }
#endif
   const TMatrixD &yU = ysvd.GetU();
   check_orthogonality(yU);
   const TMatrixD &yV = ysvd.GetV();
   check_orthogonality(yV);
   TMatrixD yfilter(nybins,nybins);
   yfilter *= 0;
   for (int eig = 0; eig < ysvd.GetNrows(); ++eig) {
      if (ysig[eig] / ysig[0] < eigenvalue_cutoff)
         break;
      double ysiginv = 1 / ysig[eig];
      for (int row = 0; row < nybins; ++row) {
         for (int col = 0; col < nybins; ++col) {
            yfilter(col,row) += yU(row,eig) * ysiginv * yV(col,eig);
         }
      }
   }

   TMatrixD Psharp(nxbins,nybins);
   Psharp.Mult(xfilter,P *= yfilter);
   std::string partname(part->GetName());
   partname += "_sharp";
   Map2D *partsharp = (Map2D*)part->Clone(partname.c_str());
   std::string title(part->GetTitle());
   title += " sharpened";
   partsharp->SetTitle(title.c_str());
   for (int ix=0; ix < nxbins; ++ix)
      for (int iy=0; iy < nybins; ++iy)
         partsharp->SetBinContent(ix+1,iy+1,Psharp(ix,iy));

   double norm = part->Integral() / partsharp->Integral();
   partsharp->Rescale(1,1,norm);
   return partsharp;
}

void ablator::check_orthogonality(const TMatrixD &matr) 
{
   // Checks if the input matrix is orthogonal, used only for debugging.

   int Nrows = matr.GetNrows();
   int Ncols = matr.GetNcols();
   for (int row1 = 0; row1 < Nrows; ++row1) {
      for (int row2 = row1; row2 < Nrows; ++row2) {
         double sum = 0;
         for (int col = 0; col < Ncols; ++col) {
            sum += matr(row1,col) * matr(row2,col);
         }
         if (row1 == row2 && fabs(sum - 1) > 1e-6) {
            std::cout << "check_orthogonality warning: diagonal element "
                      << row1 <<  " of Mtranspose M is " << sum
                      << " should be 1!" << std::endl;
         }
         else if (row1 != row2 && fabs(sum) > 1e-12) {
            std::cout << "check_orthogonality warning: off-diagonal element "
                      << row1 << "," << row2 << " of Mtranspose M is " << sum
                      << " should be 0!" << std::endl;
         }
      }
   }
}

Map2D *ablator::fillvoids(const Map2D *part, double max_void_sigmas)
{
   // This is a frontend for the Map2D::FillVoids method, with special
   // code added to deal with a peculiar weakness of the algorithm used
   // by FillVoids. Zygo images of an ablated diamond surface sometimes
   // contain narrow features whose walls are so steep that they do not
   // reflect enough light back to the interferometer for a height to be
   // measured. In such a case one typically sees a narrow spot or stripe
   // corresponding to the bottom or top of the feature, surrounded by
   // pixels without height data. Presented with such data, FillVoids tends
   // to extrapolate the smooth surrounding region into the voids right up
   // to the extremum pixels and leave behind a sharp hole or a ridge with
   // nearly vertical sides at the location of the extremum. This is not
   // physical. The algorithm below searches for these features and widens
   // them out to a transverse profile that is more consistent with their
   // having been formed by an ablation process. Both negative and positive
   // extrema are treated in a symmetric fashion.

   std::string partname(part->GetName());
   //partname += "_filled";
   Map2D *pfilled = (Map2D*)part->Clone(partname.c_str());
   //pfilled->FillVoids();

   int max_spike_width = 3;
   int max_step_height = 500;
   int max_void_span = 500;
   int min_void_span = 3;

   int nxbins = part->GetNbinsX();
   int nybins = part->GetNbinsY();
   int spike_count = 0;

   // sweep left to right over rows in x

   partname += "x";
   pfilled = (Map2D*)pfilled->Clone(partname.c_str());
   for (int iy = 1; iy <= nybins; ++iy) {
      double h1 = part->GetBinContent(1,iy);
      for (int ix = 2; ix < nxbins; ++ix) {
         double h0 = h1;
         h1 = part->GetBinContent(ix,iy);
         if (h0 == 0 && h1 != 0) {
            double hsum = h1;
            int iwidth = 1;
            while (iwidth <= max_spike_width && ix < nxbins) {
               hsum += h1 = part->GetBinContent(++ix,iy);
               if (h1 == 0)
                  break;
               ++iwidth;
            }
            if (h1 != 0)
               continue;
            double hleft = 0;
            double hright = 0;
            int ixleft = ix-iwidth-1;
            while (ixleft > 1) {
               hleft = part->GetBinContent(--ixleft,iy);
               if (hleft != 0)
                  break;
            }
            if (hleft == 0)
               continue;
            if (ix-iwidth-1-ixleft < min_void_span)
               continue;
            int ixright = ix;
            while (ixright < nxbins) {
               hright = part->GetBinContent(++ixright,iy);
               if (hright != 0)
                  break;
            }
            if (hright == 0)
               continue;
            if (ixright-ix < min_void_span)
               continue;
            if (fabs(hright - hleft) > max_step_height)
               continue;

            double hmax = hsum/iwidth;
            double xmax = ix - iwidth/2.;
            double sleft = (xmax - ixleft) * fRaster_dx / fXedgeParam;
            double zleft = exp(-0.5*sleft*sleft);
            double sright = (xmax - ixright) * fRaster_dx / fXedgeParam;
            double zright = exp(-0.5*sright*sright);

            // Solve for hicept, hslope, hgauss where
            //    h_i = hicept + hbase x_i + hgauss z_i
            //  where
            //    z_i = exp(-0.5 * [(x_i - x_0)/sigma]^2 )
            //  and i=1 marks the left, i=2 marks the right, and
            //  i=0 marks the center of the Gaussian interpolator.

            double S[3][3];
            S[0][0] = ixleft*zright - ixright*zleft;
            S[0][1] = ixright - xmax*zright;
            S[0][2] = xmax*zleft - ixleft;
            S[1][0] = zleft - zright;
            S[1][1] = zright - 1;
            S[1][2] = 1 - zleft;
            S[2][0] = ixright - ixleft;
            S[2][1] = xmax - ixright;
            S[2][2] = ixleft - xmax;
            double hicept = S[0][0]*hmax + S[0][1]*hleft + S[0][2]*hright;
            double hslope = S[1][0]*hmax + S[1][1]*hleft + S[1][2]*hright;
            double hgauss = S[2][0]*hmax + S[2][1]*hleft + S[2][2]*hright;
            double determ = S[0][0] + S[0][1] + S[0][2];
            hicept /= determ;
            hslope /= determ;
            hgauss /= determ;

            ++spike_count;

            if (ix-iwidth-1-ixleft < max_void_span) {
               for (int ivoid = ix-iwidth-1; ivoid > ixleft; --ivoid) {
                  double svoid = (xmax - ivoid) * fRaster_dx / fXedgeParam;
                  if (svoid > max_void_sigmas)
                     break;
                  double zvoid = exp(-0.5*svoid*svoid);
                  double hvoid = hicept + hslope*ivoid + hgauss*zvoid;
                  pfilled->SetBinContent(ivoid,iy,hvoid);
               }
            }
            if (ixright-ix < max_void_span) {
               for (int ivoid = ix; ivoid < ixright; ++ivoid) {
                  double svoid = (ivoid - xmax) * fRaster_dx / fXedgeParam;
                  if (svoid > max_void_sigmas)
                     break;
                  double zvoid = exp(-0.5*svoid*svoid);
                  double hvoid = hicept + hslope*ivoid + hgauss*zvoid;
                  pfilled->SetBinContent(ivoid,iy,hvoid);
               }
            }
         }
      }
   }
 
   // sweep bottom to top over column in y

   partname += "y";
   pfilled = (Map2D*)pfilled->Clone(partname.c_str());
   for (int ix = 1; ix <= nxbins; ++ix) {
      double h1 = part->GetBinContent(ix,1);
      for (int iy = 2; iy < nybins; ++iy) {
         double h0 = h1;
         h1 = part->GetBinContent(ix,iy);
         if (h0 == 0 && h1 != 0) {
            double hsum = h1;
            int iwidth = 1;
            while (iwidth <= max_spike_width && iy <= nybins) {
               hsum += h1 = part->GetBinContent(ix,++iy);
               if (h1 == 0)
                  break;
               ++iwidth;
            }
            if (h1 != 0)
               continue;
            double hlow = 0;
            double hhigh = 0;
            int iylow = iy-iwidth-1;
            while (iylow > 1) {
               hlow = part->GetBinContent(ix,--iylow);
               if (hlow != 0)
                  break;
            }
            if (hlow == 0)
               continue;
            if (iy-iwidth-1-iylow < min_void_span)
               continue;
            int iyhigh = iy;
            while (iyhigh < nybins) {
               hhigh = part->GetBinContent(ix,++iyhigh);
               if (hhigh != 0)
                  break;
            }
            if (hhigh == 0)
               continue;
            if (iyhigh-iy < min_void_span)
               continue;
            if (fabs(hhigh - hlow) > max_step_height)
               continue;

            double hmax = hsum/iwidth;
            double ymax = iy - iwidth/2.;
            double slow = (ymax - iylow) * fRaster_dy / fYedgeParam;
            double zlow = exp(-0.5*slow*slow);
            double shigh = (ymax - iyhigh) * fRaster_dy / fYedgeParam;
            double zhigh = exp(-0.5*shigh*shigh);

            // Solve for hicept, hslope, hgauss where
            //    h_i = hicept + hbase y_i + hgauss z_i
            //  where
            //    z_i = exp(-0.5 * [(y_i - y_0)/sigma]^2 )
            //  and i=1 marks the low, i=2 marks the high, and
            //  i=0 marks the center of the Gaussian interpolator.

            double S[3][3];
            S[0][0] = iylow*zhigh - iyhigh*zlow;
            S[0][1] = iyhigh - ymax*zhigh;
            S[0][2] = ymax*zlow - iylow;
            S[1][0] = zlow - zhigh;
            S[1][1] = zhigh - 1;
            S[1][2] = 1 - zlow;
            S[2][0] = iyhigh - iylow;
            S[2][1] = ymax - iyhigh;
            S[2][2] = iylow - ymax;
            double hicept = S[0][0]*hmax + S[0][1]*hlow + S[0][2]*hhigh;
            double hslope = S[1][0]*hmax + S[1][1]*hlow + S[1][2]*hhigh;
            double hgauss = S[2][0]*hmax + S[2][1]*hlow + S[2][2]*hhigh;
            double determ = S[0][0] + S[0][1] + S[0][2];
            hicept /= determ;
            hslope /= determ;
            hgauss /= determ;

            ++spike_count;

            if (iy-iwidth-1-iylow < max_void_span) {
               for (int ivoid = iy-iwidth-1; ivoid > iylow; --ivoid) {
                  double svoid = (ymax - ivoid) * fRaster_dy / fYedgeParam;
                  if (svoid > max_void_sigmas)
                     break;
                  double zvoid = exp(-0.5*svoid*svoid);
                  double hvoid = hicept + hslope*ivoid + hgauss*zvoid;
                  pfilled->SetBinContent(ix,ivoid,hvoid);
               }
            }
            if (iyhigh-iy < max_void_span) {
               for (int ivoid = iy; ivoid < iyhigh; ++ivoid) {
                  double svoid = (ivoid - ymax) * fRaster_dy / fYedgeParam;
                  if (svoid > max_void_sigmas)
                     break;
                  double zvoid = exp(-0.5*svoid*svoid);
                  double hvoid = hicept + hslope*ivoid + hgauss*zvoid;
                  pfilled->SetBinContent(ix,ivoid,hvoid);
               }
            }
         }
      }
   }
   std::cout << "Total pits drilled: " << spike_count << std::endl;
   return pfilled;
}

Map2D *ablator::discretize(const Map2D *pulse_map)
{
   // The milling methods generate pulse maps which tell how many pulses
   // should be placed on each pixel over the surface of the diamond to
   // achieve a desired target surface profile. These maps are continuous
   // but pulses must be laid down in integer multiples of pulses. This
   // function goes over a continuum map and replaces the floating point
   // pulse count in every cell with an integer value that results in a
   // minimal distortion of the ablated surface.

   int nxbins = pulse_map->GetNbinsX();
   int nybins = pulse_map->GetNbinsY();
   std::string name(pulse_map->GetName());
   name += "_int";
   Map2D *result = (Map2D*)pulse_map->Clone(name.c_str());

   double column_res[nxbins];
   for (int ix = 1; ix <= nxbins; ++ix) {
      column_res[ix] = 0;
   }
   for (int iy = 1; iy <= nybins; ++iy) {
      double row_res = 0;
      for (int ix = 1; ix <= nxbins; ++ix) {
         double z = pulse_map->GetBinContent(ix, iy);
         double zint = floor(z + column_res[ix] + row_res + 0.3);
         zint = (zint > 0)? zint : 0;
         result->SetBinContent(ix, iy, zint);
         column_res[ix] += (z - zint) * 0.9;
         row_res += (z - zint) * 0.1;
      }
   }
   return result;
}

int ablator::rasterize(const Map2D *program_map, const char *seqfile)
{
   // Takes an ablation program map and generates a pulse sequence
   // file suitable for input by the Labview laser ablation control
   // program.

   int nxbins = program_map->GetNbinsX();
   int nybins = program_map->GetNbinsY();
   double xmin = +1e9;
   double xmax = -1e9;
   double ymin = +1e9;
   double ymax = -1e9;
   int pmax = 0;
   int psum = 0;
   for (int i = 1; i <= nxbins; ++i) {
      for (int j = 1; j <= nybins; ++j) {
         int p = program_map->GetBinContent(i, j);
         if (p > 0) {
            double x = program_map->GetXaxis()->GetBinCenter(i);
            double y = program_map->GetYaxis()->GetBinCenter(j);
            xmin = (x < xmin)? x : xmin;
            ymin = (y < ymin)? y : ymin;
            xmax = (x > xmax)? x : xmax;
            ymax = (y > ymax)? y : ymax;
            pmax = (p > pmax)? p : pmax;
            psum += p;
         }
      }
   }
   std::cout << psum << " pulses generated in " << pmax << " passes."
             << std::endl;

   int width = ((xmax - xmin) / fRaster_dx) + 1;
   int width32 = (int)ceil(width / 32.);
   for (int pass = 1; pass <= pmax; ++pass) {
      std::string fname(seqfile);
      std::size_t p_seq = fname.rfind(".seq");
      if (p_seq == fname.npos) {
         std::stringstream str;
         str << seqfile << "_" << pass << ".seq";
         fname = str.str();
      }
      else {
         std::stringstream str;
         str << fname.substr(0, p_seq) << "_" << pass << ".seq";
         fname = str.str();
      }
      std::ofstream ofs(fname.c_str());
      if (!ofs.is_open()) {
         std::cerr << "Error in ablator::mill - unable to open output file "
                   << fname << ", cannot continue." << std::endl;
         return 0;
      }

      int pulses_this_pass = 0;
      int row = 0;
      for (double y = ymin; y <= ymax; y += fRaster_dy) {
         row++;
         int j = program_map->GetYaxis()->FindBin(y);
         int parray[width32];
         for (int n = 0; n < width32; ++n) {
            parray[n] = 0;
         }
         int pcount = 0;
         int col = 0;
         for (double x = xmin; x <= xmax; x += fRaster_dx) {
            col++;
            int i = program_map->GetXaxis()->FindBin(x);
            if (program_map->GetBinContent(i, j) >= pass) {
               int n32 = (col - 1) / 32;
               int b32 = (col - 1) % 32;
               parray[n32] |= (1 << b32);
               ++pcount;
            }
         }
         if (pcount > 0) {
            std::stringstream seqstr("");
            char seqline[100];
            sprintf(seqline, "%4d\t%6.3lf\t%6.3lf\t%6.3lf\t%6.3lf\t%6d",
                    row, y, y + fRaster_dy, xmin, xmax, width);
            seqstr << seqline;
            for (int n = 0; n < width32; ++n) {
               seqstr << " " << std::hex << std::setw(8) << parray[n];
            }
            ofs << seqstr.str() << std::endl;
            pulses_this_pass += pcount;
         }
      }
      std::cout << "writing " << pulses_this_pass
                << " pulses to file " << fname << std::endl;
   }
   return psum;
}

Map2D *ablator::detrench(const Map2D *program_map)
{
   // Takes a generated program map and modifies it to account for the
   // initial variation in laser pulse energy in a pulse sequence. The
   // input program_map is not modified, and the resulting program is
   // returned.
   //
   // Programmer's notes:
   // The algorithms used by the ablator class to generate a map of laser
   // pulses that produces a final surface with a desired shape profile
   // assume that the laser pulses are all identical. In actual fact, this
   // assumption is violated in two ways.
   //   1) The initial string of pulses generated after the laser has been
   //      off for a few-second period have significantly higher energy
   //      than their steady-state average.
   //   2) There are random pulse-to-pulse variations in the laser energy
   //      even after the output reaches its steady-state average.
   // The effects of (2) tend to cancel out on average and result only in
   // the surface being somewhat rougher than predicted by the model. The
   // effects of (1) are more serious because they all tend to pile up at
   // edges of the cut region, at the ends of the raster pattern where the
   // laser pulse sequence starts and ends. The effect is seen at both ends
   // because of the serpentine motion pattern, and appears as trenches up
   // to 20 microns deep that appear along the left and right boundaries 
   // of the cut region. The pulse energy is modeled as an exponential 
   // decay from an initial level to the steady-state level as
   //      E(t) = E0 (1 + fTaperFactor exp(-t / fTaperLength))
   // where E0 is the steady-state pulse energy and t is time in units
   // of pulse count (20ms at 50Hz), where the first pulse appears at
   // time t=0.
   //
   // To correct ("detrench") the program, some assumptions need to be made
   // regarding the sequence in which the program pulses are laid down. They
   // are stored in the program map as a total pulse count per pixel, but in
   // actual fact they are not laid down by sitting at a single pixel and
   // dropping N pulses, then moving to the next pixel. Instead, they are
   // laid down in a series of passes over the entire surface, where each
   // pixel receives one pulse per pass. This allows the control program to
   // move the sample at a fixed motor speed, while the laser is pulsing at
   // a fixed frequency, allowing several hundred pulses to be generated at
   // a time without pausing. Correcting the program to account for the power
   // taper effect requires taking into account that the t=0 point appears at
   // a different location on each pass. Let r(t) be the pulse repeat count
   // pattern for a given row of the program, and rcorr(t) be the effective
   // pulse repeat count that the program will produce due to the power taper
   // effect.
   //      rcorr(t) = (1 + fTaperFactor) r(t) - 
   //                 fTaperFactor exp(-t / fTaperLength) / fTaperLength
   //                 Integral[0,t] {r(t') exp(+t' / fTaperLength)} dt'
   //
   // To detrench the program, I take the input program to give rcorr(t) for
   // each pixel, and I invert the above integral equation to solve for the
   // r(t) that is needed in order to achieve it. The result is
   //     r(t) = s(t) exp(t / tau) exp(-t / fTaperLength)
   // where
   //     tau = fTaperLength (1 + fTaperFactor) / fTaperFactor
   // and
   //     s(t) = Integral[0,t] {f'(t') exp(-t' / tau)} dt' / (1 + fTaperFactor)
   // where
   //     f'(t) = d/dt f(t)
   // and 
   //     f(t) = rcorr(t) exp(t / fTaperLength)

   double effTaperFactor = fTaperFactor;
   Map2D *result = (Map2D*)program_map->Clone("program_detrenched");
   int nxbins = result->GetNbinsX();
   int nybins = result->GetNbinsY();
   double tauinv = effTaperFactor / (fTaperLength * (1 + effTaperFactor));
   for (int iy = 1; iy <= nybins; iy += 1) {
      // Take a serpentine sweep pattern across the sample, with the even
      // rows sweeping left-to-right and the odd rows right-to-left.
      int ix_start;
      int ix_stop;
      int ix_step;
      if (iy % 2 == 0) {
         ix_start = 1;
         ix_stop = nxbins + 1;
         ix_step = 1;
      }
      else {
         ix_start = nxbins;
         ix_stop = 0;
         ix_step = -1;
      }
      int t = 0;
      double lastf = 0;
      double sintegral = 0;
      for (int ix = ix_start; ix != ix_stop; ix += ix_step) {
         t += 1;
         double rcorr = program_map->GetBinContent(ix,iy);
         double f = rcorr * exp(t / fTaperLength);
         sintegral += (f - lastf) * (fTaperLength / effTaperFactor) *
                      (exp(-(t-1) * tauinv) - exp(-t * tauinv));
         double r = sintegral * exp(t * tauinv) * exp(-t / fTaperLength);
         r = (r < rcorr)? r : rcorr;
         result->SetBinContent(ix, iy, r);
         lastf = f;
      }
   }
   return result;
}

Map2D *ablator::simulate_mill(const Map2D *program_map, Map2D *part_map)
{
   // Simulate laser ablation with program_map by removing material
   // one pulse at a time from part_map. The input maps are not modified.
   // A clone of part_map containing the milled part surface profile is
   // returned.

   std::string partname(part_map->GetName());
   partname += "_silled";
   Map2D *result = (Map2D*)part_map->Clone(partname.c_str());

   // create a local image of the ablation spot map with the resolution of the
   // program_map, and store it in a local array for efficient access below.

   double dx = result->GetXaxis()->GetBinWidth(1);
   double dy = result->GetYaxis()->GetBinWidth(1);
   int ndx = (int)floor(fSpotHsize / dx);
   int ndy = (int)floor(fSpotVsize / dy);
   Map2D *spot = make_beam_spot(result);
   spot->SetName("sim_spot");
   spot->SetTitle("simulated milling tool spot");
   double pulse_depth[ndy + 1][ndx + 1];
   for (int jx = 0; jx <= ndx / 2; ++jx) {
      for (int jy = 0; jy <= ndy / 2; ++jy) {
         int ix1 = jx + 1;
         int ix2 = spot->GetNbinsX() - jx;
         int iy1 = jy + 1;
         int iy2 = spot->GetNbinsY() - jy;
         pulse_depth[jy][jx] = spot->GetBinContent(ix1,iy1);
         pulse_depth[jy][ndx-jx] = spot->GetBinContent(ix2,iy1);
         pulse_depth[ndy-jy][ndx-jx] = spot->GetBinContent(ix2,iy2);
         pulse_depth[ndy-jy][jx] = spot->GetBinContent(ix1,iy2);
      }
   }

   // now apply a series of pulses to the part_map and show the result

   int nxbins = result->GetNbinsX();
   int nybins = result->GetNbinsY();
   double count = 0;
   for (int py = 1; py <= program_map->GetNbinsY(); ++py) {
      double yrast = program_map->GetYaxis()->GetBinCenter(py);
      int iy0 = result->GetYaxis()->FindBin(yrast);
      double t = 0;
      double rintegral = 0;
      int px0 = (py % 2 == 1)? 1 : program_map->GetNbinsX();
      int px1 = (py % 2 == 1)? program_map->GetNbinsX() : 1;
      int dpx = (py % 2 == 1)? 1 : -1;
      for (int px = px0; px != px1 + dpx; px += dpx) {
         double xrast = program_map->GetXaxis()->GetBinCenter(px);
         int ix0 = result->GetXaxis()->FindBin(xrast);
         double reps = program_map->GetBinContent(px,py);

         // The laser pulse power is not constant across a single sweep 
         // of the diamond, but starts off high and tapers off exponentially
         // to a steady-state level after some number of pulses. The excess
         // power above the steady-state level is parameterized by the factor
         // fTaperFactor and the lifetime in pulses is fTaperLength. The pulse
         // energy model for a single sweep is given by the formula,
         //      E(t) = E0 (1 + fTaperFactor exp(-t / fTaperLength))
         // where E0 is the steady-state pulse energy and t is time in units
         // of pulse count (20ms at 50Hz).
         //
         // In the implementation of this model below, I take a serpentine
         // sweep pattern across the sample, with the even rows sweeping left-
         // to-right and the odd rows right-to-left. In one iteration over all
         // pixels in the sample, the algorithm removes material corresponding
         // to as many whole-integer passes as it can, which results in each
         // pixel seeing rep pulses. To implement the power-taper model, I need
         // to sum the above formula to account for the varying rep value over
         // the sample. The pulse power correction factor r(t) is
         //      rcorr(t) = (1 + fTaperFactor) r(t) - 
         //                 fTaperFactor exp(-t / fTaperLength) / fTaperLength
         //                 Integral[0,t] {r(t') exp(+t' / fTaperLength)} dt'

#if APPLY_DETRENCH_ALGORITHM
         t += 1;
         rintegral += reps * (exp(t / fTaperLength) -
                              exp((t-1) / fTaperLength));
         double power_factor =  1 + fTaperFactor -
                fTaperFactor * exp(-t / fTaperLength) * rintegral / reps;
#else
         double power_factor = 1;
#endif

         for (int jy = 0; jy <= ndy / 2; ++jy) {
            for (int jx = 0; jx <= ndx / 2; ++jx) {
               double depth1 = reps * pulse_depth[jy][jx] * power_factor;
               double depth2 = reps * pulse_depth[jy][ndx-jx] * power_factor;
               double depth3 = reps * pulse_depth[ndy-jy][ndx-jx] * power_factor;
               double depth4 = reps * pulse_depth[ndy-jy][jx] * power_factor;
               if (ix0 + jx <= nxbins && iy0 + jy <= nybins) {
                  result->SetBinContent(ix0 + jx, iy0 + jy, 
                  result->GetBinContent(ix0 + jx, iy0 + jy) - depth1);
               }
               if (ix0 - jx > 1 && jy + iy0 <= nybins) {
                  result->SetBinContent(ix0 - jx - 1, iy0 + jy,
                  result->GetBinContent(ix0 - jx - 1, iy0 + jy) - depth2);
               }
               if (ix0 - jx > 1 && iy0 - jy > 1) {
                  result->SetBinContent(ix0 - jx - 1, iy0 - jy - 1,
                  result->GetBinContent(ix0 - jx - 1, iy0 - jy - 1) - depth3);
               }
               if (ix0 + jx <= nxbins && iy0 - jy > 1) {
                  result->SetBinContent(ix0 + jx, iy0 - jy - 1,
                  result->GetBinContent(ix0 + jx, iy0 - jy - 1) - depth4);
               }
            }
         }
         count += reps;
      }
   }
   std::cout << "milling finished after " << count << " pulses"
             << std::endl;
   return result;
}