//
// 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 <TRandom.h>
#include <TFFTComplexReal.h>
#include <TDecompSVD.h>

#include "ablator.h"

#define SIMULATE_MILLING 0
#define SMOOTH_REVERSIBILITY_CHECK 0

Map2D *ablator::mill(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.

   double cut_depth = (fSpotDepth*2*M_PI)*(fSpotHsize/fRaster_dx)
                                         *(fSpotVsize/fRaster_dy);
   std::cout << "ablator::mill - per-pass cut depth is " 
             << cut_depth
             << std::endl;

   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 nxbins = (int)floor((xmax - xmin)/fRaster_dx + 0.5);
   int nybins = (int)floor((ymax - ymin)/fRaster_dy + 0.5);
   xmax = xmin + fRaster_dx * nxbins;
   ymax = ymin + fRaster_dy * nybins;
   Map2D *shape = new Map2D(TH2D("shape","ablation map",
                                 nxbins, xmin, xmax, nybins, ymin, ymax));
 
   // form the shape of the beam spot and take its DFT

   Map2D *spot = spot_filter(shape);
   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

   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-6) {
            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)? par - tar : 0);
      }
   }
   Map2D *program = (Map2D*)shape->Clone("program");
   program->SetName("program");
   program->SetTitle("milling program");

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

   Map2D *programFFTr = (Map2D*)program->Clone("programFFTr");
   program->FFT(programFFTr,"RE R2C P");
   Map2D *programFFTi = (Map2D*)program->Clone("programFFTi");
   program->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*)program->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);

#if SIMULATE_MILLING
   result->Fill(part_map);
   Map2D *milled = simulate_mill(program, result);
   result->Fill(milled);
   delete milled;
#else
   result->Fill(part_sharp);
   result->Add(program, -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;
#endif

   // rescale program to units of laser pulses

   program->Rescale(1, 1, 1/cut_depth);
   delete spot;
   delete spotFFTr;
   delete spotFFTi;
   delete invFFT;

   return result;
}

Map2D *ablator::mill_sparse(Map2D *part_map, Map2D *target_map)
{
   // DEPRECATED - RTJ, December 29, 2015
   // This is a revision of ablator::mill for use when the part_map and/or
   // target_map is sparse, that is, contains many pixels with unknown
   // height values.  This routine is meant to handle maps generated from
   // surface profile images which generally have incomplete coverage of
   // the surface.  Similar to ablator::mill, the goal is to generate an
   // ablation pulse sequence that mills the part_map surface down to
   // match the target_map as closely as possible, considering only the
   // regions of the two maps that are filled.  Unknown (matte) pixels in
   // these maps are those with zero values.  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 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::string name(part_map->GetName());
   name += "_milled";
   Map2D *result = (Map2D*)part_map->Clone(name.c_str());

   // Repeat the following simple whack-a-mole algorithm until there
   // are no longer any high spots above the target_map level.  All
   // pulses are required to align with the nodes of the raster grid.
   // If any pulses would create a dip below the target_map level in
   // any pixel then they are forbidden, so use a double-buffer strategy
   // for material removal.

   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 nxbins = (int)floor((xmax-xmin)/fRaster_dx + 0.5);
   int nybins = (int)floor((ymax-ymin)/fRaster_dy + 0.5);
   xmax = xmin + fRaster_dx*nxbins;
   ymax = ymin + fRaster_dy*nybins;
   Map2D *shape = new Map2D(TH2D("shape","ablation map",
                                 nxbins,xmin,xmax,nybins,ymin,ymax));
   Map2D *trial = new Map2D(TH2D("trial_shape","trial ablation map",
                                 nxbins,xmin,xmax,nybins,ymin,ymax));

   // To sample part_map and target_map,
   // first rebin them onto the raster grid

   Map2D *part_remap = new Map2D(TH2D("part_remap","part map rebinned",
                                      nxbins,xmin,xmax,nybins,ymin,ymax));
   Map2D *target_remap = new Map2D(TH2D("target_remap","target map rebinned",
                                   nxbins,xmin,xmax,nybins,ymin,ymax));
   Map2D *hcount = new Map2D(TH2D("hcount","",
                                  nxbins,xmin,xmax,nybins,ymin,ymax));
   hcount->Flood(1e-30);
   for (int ix = 1; ix <= part_map->GetNbinsX(); ++ix) { 
      for (int iy = 1; iy <= part_map->GetNbinsY(); ++iy) { 
         double x = part_map->GetXaxis()->GetBinCenter(ix);
         double y = part_map->GetYaxis()->GetBinCenter(iy);
         double height = part_map->GetBinContent(ix,iy);
         if (height != 0) {
            part_remap->Fill(x,y,height);
            hcount->Fill(x,y);
         }
      }
   }
   part_remap->Divide(hcount);
   delete hcount;
   for (int ix = 1; ix <= target_map->GetNbinsX(); ++ix) { 
      for (int iy = 1; iy <= target_map->GetNbinsY(); ++iy) { 
         double height = target_map->GetBinContent(ix,iy);
         double x = target_map->GetXaxis()->GetBinCenter(ix);
         double y = target_map->GetYaxis()->GetBinCenter(iy);
         int bin = target_remap->FindBin(x,y);
         if (height != 0 && height > target_remap->GetBinContent(bin)) {
            target_remap->SetBinContent(bin,height);
         }
      }
   }

   for (int ix = 1; ix <= nxbins; ++ix) { 
      for (int iy = 1; iy <= nybins; ++iy) { 
         double par = part_remap->GetBinContent(ix,iy);
         double tar = target_remap->GetBinContent(ix,iy);
         if (par != 0 && tar != 0) {
            shape->SetBinContent(ix,iy,(tar < par-fSpotDepth)? par-tar : 0);
         }
      }
   }

   // save the pulse sequence in a list, for sorting in raster order later

   std::list<Pulse> pulses;
   int pulse_count = 0;

   while (1 > 0) {
      Pulse p;
      double top_peak = -1;
      for (int iy = 1; iy <= nybins; ++iy) {
         for (int ix = 1; ix <= nxbins; ++ix) {
            double height = shape->GetBinContent(ix,iy);
            trial->SetBinContent(ix,iy,height);
            if (height != 0 && height > top_peak) {
               top_peak = height;
               p.ix = ix;
               p.iy = iy;
            }
         }
      }
      if (top_peak < fSpotDepth) {
         break;
      }

      int n5dx = (int)floor(5./fRaster_dx+0.5);
      int n5dy = (int)floor(5./fRaster_dy+0.5);
      double dx = fRaster_dx/fSpotHsize;
      double dy = fRaster_dy/fSpotVsize;
      bool pulse_forbidden = false;
      int bogging_down = 0;
      for (int ix = -n5dx; ix <= +n5dx && !pulse_forbidden; ++ix) {
         for (int iy = -n5dy; iy <= +n5dy && !pulse_forbidden; ++iy) {
            double height = shape->GetBinContent(ix+p.ix,iy+p.iy);
            if (height != 0) {
               height -= fSpotDepth * exp(-0.5*(dx*dx*ix*ix + dy*dy*iy*iy));
               trial->SetBinContent(ix+p.ix,iy+p.iy,height);
               if (height < 0) {
                  ++bogging_down;
               }
            }
            if (bogging_down > 12) {
               pulse_forbidden = true;
            }
         }
      }
      if (pulse_forbidden) {
         double height = shape->GetBinContent(p.ix,p.iy);
         shape->SetBinContent(p.ix,p.iy,height-fSpotDepth);
      }
      else {
         pulses.push_back(p);
         ++pulse_count;
         Map2D *tmp = shape;
         shape = trial;
         trial = tmp;
      }
      if (pulse_count/10000 * 10000 == pulse_count) {
         std::cout << pulse_count << " pulses defined, latest at "
                   << p.ix << "," << p.iy << std::endl;
      }
   }

   std::cout << pulse_count << " pulses defined, latest at "
             << pulses.back().ix << "," << pulses.back().iy << std::endl;
   pulses.sort();

   // now apply the ordered series of pulses to the part_map and
   // return the result

   std::list<Pulse>::iterator iter;
   pulse_count = 0;
   int ix = -1;
   int iy = -1;
   for (iter = pulses.begin(); iter != pulses.end(); ++iter) {
      if (iter->ix == ix && iter->iy == iy) {
         ++ pulse_count;
         continue;
      }
      if (pulse_count > 0) {
         double x = shape->GetXaxis()->GetBinCenter(ix);
         double y = shape->GetYaxis()->GetBinCenter(iy);
         int jx0 = result->GetXaxis()->FindBin(x);
         int jy0 = result->GetYaxis()->FindBin(y);
         double dx = result->GetXaxis()->GetBinWidth(jx0)/fSpotHsize;
         double dy = result->GetYaxis()->GetBinWidth(jy0)/fSpotVsize;
         int j5dx = (int)floor(5./dx+0.5);
         int j5dy = (int)floor(5./dy+0.5);
         for (int jx = -j5dx; jx <= +j5dx; ++jx) {
            for (int jy = -j5dy; jy <= +j5dy; ++jy) {
               double height = result->GetBinContent(jx+jx0,jy+jy0);
               if (height != 0) {
                  height -= pulse_count * fSpotDepth *
                            exp(-0.5*(dx*dx*jx*jx + dy*dy*jy*jy));
                  result->SetBinContent(jx+jx0,jy+jy0,height);
               }
            }
         }
      }
      pulse_count = 1;
      ix = iter->ix;
      iy = iter->iy;
   }
   std::cout << "milling finished after " << pulse_count << " pulses "
             << std::endl;
   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 the
   // Gaussian beam focal spot profile.  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.

   // form the shape of the beam spot and take its DFT
   // Here I use a slightly larger smoothing function than the actual
   // beam spot dimensions. This stabilizes the deconvolution with the
   // beam spot, which otherwise becomes very noisy and subject to 
   // ringing effects (Gibbs phenomenon) close to sharp edges.

   int nxbins = target->GetNbinsX();
   int nybins = target->GetNbinsY();
   double sigma_x = getSpotHsize();
   double sigma_y = getSpotVsize();
   setSpotSize(sigma_x * smear_factor, sigma_y * smear_factor);
   Map2D *fatspot = spot_filter(target);
   setSpotSize(sigma_x, sigma_y);
   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 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 *spot = spot_filter(target);
   Map2D *spotFFTr = (Map2D*)spot->Clone("spotFFTr");
   spot->FFT(spotFFTr,"RE R2C P");
   Map2D *spotFFTi = (Map2D*)spot->Clone("spotFFTi");
   spot->FFT(spotFFTi,"IM R2C P");
   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 = spotFFTr->GetBinContent(ix,iy);
         double spot_i = spotFFTi->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-3) {
            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 spot;
   delete spotFFTr;
   delete spotFFTi;
   delete resultFFTr;
   delete resultFFTi;

#endif

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

   return result;
}

Map2D *ablator::unsmooth(Map2D *part, double min_eigenvalue)
{
   // 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 = spot_filter(part);
   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,1+colrel);
      }
   }
   TDecompSVD xsvd(xspot);
   if (! xsvd.Decompose()) {
      std::cerr << "ablator::unsmooth - 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();
   const TMatrixD &xV = xsvd.GetV();
   TMatrixD xfilter(nxbins,nxbins);
   xfilter *= 0;
   for (int eig = 0; eig < xsvd.GetNrows(); ++eig) {
      if (xsig[eig] < min_eigenvalue * xsig[0]) 
         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::unsmooth - 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();
   const TMatrixD &yV = ysvd.GetV();
   TMatrixD yfilter(nybins,nybins);
   yfilter *= 0;
   for (int eig = 0; eig < ysvd.GetNrows(); ++eig) {
      if (ysig[eig] < min_eigenvalue * ysig[0]) 
         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;
}

Map2D *ablator::fillvoids(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 / fSpotHsize;
            double zleft = exp(-0.5*sleft*sleft);
            double sright = (xmax - ixright) * fRaster_dx / fSpotHsize;
            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 / fSpotHsize;
                  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 / fSpotHsize;
                  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 / fSpotVsize;
            double zlow = exp(-0.5*slow*slow);
            double shigh = (ymax - iyhigh) * fRaster_dy / fSpotVsize;
            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 / fSpotVsize;
                  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 / fSpotVsize;
                  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::spot_filter(const Map2D *shape) const
{
   int nxbins = shape->GetNbinsX();
   int nybins = shape->GetNbinsY();
   double delx = shape->GetXaxis()->GetBinWidth(1);
   double dely = shape->GetYaxis()->GetBinWidth(1);
   Map2D *spot = (Map2D*)shape->Clone("spot");
   spot->SetTitle("beam spot");
   double norm = 0;
   for (int ix = 0; ix <= nxbins/2; ++ix) {
      for (int iy = 0; iy <= nybins/2; ++iy) {
         double r = pow(delx*ix/fSpotHsize,2) + pow(dely*iy/fSpotVsize,2);
         double s = exp(-r/2);
         spot->SetBinContent(ix+1,iy+1,s);
         norm += s;
         if (ix > 0) {
            spot->SetBinContent(nxbins-ix+1,iy+1,s);
            norm += s;
         }
         if (iy > 0) {
            spot->SetBinContent(ix+1,nybins-iy+1,s);
            norm += s;
         }
         if (ix > 0 && iy > 0) {
            spot->SetBinContent(nxbins-ix+1,nybins-iy+1,s);
            norm += s;
         }
      }
   }
   spot->Rescale(1, 1, 1/norm);
   return spot;
}

Map2D *ablator::discretize(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);
         if (z > fSpotDepth) {
            double zint = floor(z + column_res[ix] + row_res + 0.5);
            zint = (zint > 0)? zint : 0;
            result->SetBinContent(ix, iy, zint);
            column_res[ix] += (z - zint)/2;
            row_res += (z - zint)/2;
         }
         else {
            result->SetBinContent(ix, iy, 0);
         }
      }
   }
   return result;
}

int ablator::rasterize(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(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(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*)program_map->Clone(partname.c_str());
   int nxbins = program_map->GetNbinsX();
   int nybins = program_map->GetNbinsY();
   result->Fill(part_map);

   // build a local model of the pulse depth profile
   // only one quadrant is needed because of symmetry

   double dx = result->GetXaxis()->GetBinWidth(1) / fSpotHsize;
   double dy = result->GetYaxis()->GetBinWidth(1) / fSpotVsize;
   int ndx = (int)floor(8/dx + 0.5);
   int ndy = (int)floor(8/dy + 0.5);
   double pulse_depth[ndy + 1][ndx + 1];
   for (int jx = 0; jx <= ndx; ++jx) {
      for (int jy = 0; jy <= ndy; ++jy) {
         double dr2 = dx*dx * jx*jx + dy*dy * jy*jy;
         double depth = fSpotDepth * exp(-0.5 * dr2);
         pulse_depth[jy][jx] = depth;
      }
   }

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

   int count = 0;
   int pass = 0;
   int lastcount = -1;
   Map2D *residue = (Map2D*)program_map->Clone("residue");
   while (count > lastcount) {
      lastcount = count;
      for (int iy = 1; iy <= nybins; iy += 1) {

         // 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'

         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;
         }
         double t = 0;
         double rintegral = 0;
         for (int ix = ix_start; ix != ix_stop; ix += ix_step) {
            double pulses_to_go = residue->GetBinContent(ix,iy);
            if (pulses_to_go > 1.0e-10) {
               float reps = (int)floor(pulses_to_go + 0.5);
               reps = pulses_to_go; // taking out digitization here
               reps = (reps > fReps)? fReps : reps;
               double x = residue->GetXaxis()->GetBinCenter(ix);
               double y = residue->GetYaxis()->GetBinCenter(iy);
               residue->SetBinContent(ix, iy, pulses_to_go - reps);
               int jx0 = result->GetXaxis()->FindBin(x);
               int jy0 = result->GetYaxis()->FindBin(y);
               t += 1;
#if APPLY_DETRENCH_ALGORITHM
               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; ++jy) {
                  for (int jx = 0; jx <= ndx; ++jx) {
                     double depth = reps * pulse_depth[jy][jx] * power_factor;
                     if (jx0 + jx <= nxbins && jy0 + jy <= nybins) {
                        result->SetBinContent(jx0 + jx, jy0 + jy, 
                        result->GetBinContent(jx0 + jx, jy0 + jy) - depth);
                     }
                     if (jx > 0 && jx0 - jx > 0 && jy + jy0 <= nybins) {
                        result->SetBinContent(jx0 - jx, jy0 + jy,
                        result->GetBinContent(jx0 - jx, jy0 + jy) - depth);
                     }
                     if (jx > 0 && jx0 - jx > 0 && jy > 0 && jy0 - jy > 0) {
                        result->SetBinContent(jx0 - jx, jy0 - jy,
                        result->GetBinContent(jx0 - jx, jy0 - jy) - depth);
                     }
                     if (jx0 + jx <= nxbins && jy > 0 && jy0 - jy > 0) {
                        result->SetBinContent(jx0 + jx, jy0 - jy,
                        result->GetBinContent(jx0 + jx, jy0 - jy) - depth);
                     }
                  }
               }
               count += reps;
            }
         }
      }
      ++pass;
      std::cout << "completed pass " << pass
                << " after " <<  count << " pulses" << std::endl;
   }
   std::cout << "milling finished after " << count << " pulses, "
             << pass << " passes." << std::endl;

   return result;
}