//
// Couple.C - implementation for Couple class.
//
// author: richard.t.jones at uconn.edu
// version: june 28, 2019
//
//////////////////////////////////////////////////////////////////////////
//                                                                      //
// Couple class                                                         //
//                                                                      //
// The Couple class provides a robust toolkit for relating pairs of     //
// rocking curve topographical maps within the root framework. It is    //
// based on the underlying Map2D class that extends the basic ROOT      //
// TH2D class with functionality related to treating 2D histograms as   //
// representations of 2D surfaces. A Couple consists of two rocking     //
// curve topographs of the same Bragg reflection from the same sample   //
// taken along orthogonal directions, within which a projection of the  //
// crystal strain pattern across the entire crystal onto a 2D surface   //
// is visualized.                                                       //
//                                                                      //
// The class currently supports the following on pairs of maps:         //
//        * association of complementary pairs of scans                 //
//        * arbitrary transformations (shifts, rotations etc.)          //
//        * masking of blank or random portions of the maps             //
//        * taking the curl, divergence in 2D                           //
//        * solving Poisson's equation in 2D                            //
//        * plotting pretty pictures                                    //
//                                                                      //
// This package was developed at the University of Connecticut by       //
// Richard T. Jones, as part of the analysis suite for the analysis     //
// of diamond crystal rocking curve scans at taken at CHESS and CLS,    //
// for the purpose of assessing coherent bremsstrahlung radiators for   //
// the GlueX experiment at Jefferson Lab.                               //
//                                                                      //
// This work is supported by the U.S. National Science Foundation       //
// under grant 1812415.                                                 //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "Couple.h"

#include <unistd.h>
#include <iostream>
#include <fstream>
#include <string>

#include <TROOT.h>
#include <TCanvas.h>
#include <TExec.h>
#include <TStyle.h>

int nsearchpaths = 2;
std::string searchpath[] = 
{
   "",
   "root://nod29.phys.uconn.edu//Gluex/beamline/diamonds/cls-6-2019/results/"
};

Couple *Couple::fPicking = 0;

ClassImp(Couple);

Couple::Couple()
 : Couple("", "")
{
   // default constructor
}

Couple::Couple(const char *name, const char *title)
 : TNamed(name, title)
{
   // barebones constructor
 
   for (int i=0; i<2; i++) {
      fXshift[i] = 0;
      fYshift[i] = 0;
      fPsirot[i] = 0;
      fXscale[i] = 1;
      fYscale[i] = 1;
      fXYskew[i] = 0;
      fChirot[i] = 0;
      fPhirot[i] = 0;
      fCleanMask[i] = 0;
      fMap[i] = 0;
   }
   fZeroMask = 0;
   fXrange[0] = 0;
   fXrange[1] = 8.5;
   fPixels[0] = 850;
   fYrange[0] = 0;
   fYrange[1] = 8.5;
   fPixels[1] = 850;
}

Couple::Couple(const Couple &src)
 : TNamed(src.GetName(), src.GetTitle())
{
   // copy constructor
 
   copy(src);
}

Couple::~Couple()
{
   // destructor
}

Couple &Couple::operator=(const Couple &src)
{
   // assignment operator

   copy(src);
   return *this;
}

void Couple::copy(const Couple &src)
{
   // safe copy of state from src to *this

   if (src.fZeroMask) {
      fZeroMask = new Map2D(*src.fZeroMask);
      fZeroMask->SetDirectory(0);
   }
   for (int p=0; p<2; p++) {
      fResultsPath[p] = src.fResultsPath[p];
      fResultsName[p] = src.fResultsName[p];
      fXshift[p] = src.fXshift[p];
      fYshift[p] = src.fYshift[p];
      fPsirot[p] = src.fPsirot[p];
      fXscale[p] = src.fXscale[p];
      fYscale[p] = src.fYscale[p];
      fXYskew[p] = src.fXYskew[p];
      fChirot[p] = src.fChirot[p];
      fPhirot[p] = src.fPhirot[p];
      fXrange[p] = src.fXrange[p];
      fYrange[p] = src.fYrange[p];
      fPixels[p] = src.fPixels[p];
      if (src.fCleanMask[p]) {
         fCleanMask[p] = new Map2D(*src.fCleanMask[p]);
         fCleanMask[p]->SetDirectory(0);
      }
      if (src.fMap[p]) {
         fMap[p] = new Map2D(*src.fMap[p]);
         fMap[p]->SetDirectory(0);
      }
      else {
         fMap[p] = 0;
      }
   }
}

void Couple::select(int p, const char *resultsfile)
{
   // Select a scan results file (XXX_results.root) for either
   // the first (p=0) or second (p=1) member of the pair.

   TFile *fresults = open(resultsfile);
   if (fresults) {
      fResultsPath[p] = TString(fresults->GetName());
      deletemap(p);
      if (fCleanMask[p]) {
         delete fCleanMask[p];
         fCleanMask[p] = 0;
      }
   }
}

void Couple::setresol(int xpixels, int ypixels)
{
   // Set the resolution of the images to be generated,
   // in terms of the width (xpixels) and height (ypixels) in pixels.

   fPixels[0] = xpixels;
   fPixels[1] = ypixels;
   deletemap(0);
   deletemap(1);
}

void Couple::setxrange(double xmax, double xmin)
{
   // Set the x coordinate range of the images to be generated,
   // in mm.
 
   fXrange[0] = xmin;
   fXrange[1] = xmax;
   deletemap(0);
   deletemap(1);
}

void Couple::setyrange(double ymax, double ymin)
{
   // Set the y coordinate range of the images to be generated,
   // in mm.

   fYrange[0] = ymin;
   fYrange[1] = ymax;
   deletemap(0);
   deletemap(1);
}

void Couple::shift(int p, double xshift, double yshift)
{
   // Set the offset in x (xshift) and y (yshift) of the images to be generated
   // for member p=0 or p=1, relative to the original topographs, in units of mm.

   fXshift[p] = xshift;
   fYshift[p] = yshift;
   deletemap(p);
}

void Couple::invert(int p, int ix, int iy)
{
   // Request that the generated images for member p=0 or p=1 be inverted
   // in x (ix != 0) and/or y (iy != 0) relative to how they appear in the
   // original topographs.

   if (ix) {
      fXscale[p] = -fabs(fXscale[p]);
      fXYskew[p] = -fabs(fXYskew[p]);
   }
   else {
      fXscale[p] = fabs(fXscale[p]);
      fXYskew[p] = fabs(fXYskew[p]);
   }
   if (iy) {
      fYscale[p] = -fabs(fYscale[p]);
      fXYskew[p] *= -1;
   }
   else
      fYscale[p] = fabs(fYscale[p]);
   deletemap(p);
}

void Couple::rotate(int p, double phi, int deg)
{
   // Request that the generated images for member p=0 or p=1 be rotated in
   // the xy plane by angle phi (radians) relative to how they appear in the
   // original topographs.
 
   if (deg)
      phi *= M_PI/180;
   fPsirot[p] = fmod(phi, 2*M_PI);
   deletemap(p);
}

void Couple::zoom(int p, double fzoom)
{
   // Request that the generated images for member p=0 or p=1 be zoomed by
   // factor fzoom the xy plane relative to how they appear in the
   // original topographs.
 
   fXscale[p] = fzoom;
   fYscale[p] = fzoom;
   deletemap(p);
}

void Couple::stretch(int p, double aspectratio)
{
   // Request that the generated images for member p=0 or p=1 be stretched
   // by aspect ratio aspectratio (width / height) relative to how they 
   // appear in the original topographs.
 
   double scalefactor = sqrt(fabs(aspectratio));
   fXscale[p] = (fXscale[p] > 0)? scalefactor : -scalefactor;
   fYscale[p] = (fYscale[p] > 0)? 1/scalefactor : -1/scalefactor;
   fXYskew[p] = 0;
   deletemap(p);
}

void Couple::setmask(const Map2D *mask)
{
   // Request that all zero pixels in the Map2D object mask be zeroed in
   // the generated images for either of the two members of this set.
   // The caller retains ownership of the passed object.
 
   if (fZeroMask) {
      delete fZeroMask;
      fZeroMask = 0;
   }
   if (mask == 0)
      return;
   TH2D hmap("hmask", "mask image", fPixels[0], fXrange[0], fXrange[1],
                                    fPixels[1], fXrange[0], fXrange[1]);
   fZeroMask = new Map2D(hmap);
   fZeroMask->SetDirectory(0);
   fZeroMask->Fill(mask);
   for (int p=0; p<4; ++p)
      deletemap(p);
}

void Couple::setclean(int p, Map2D *mask)
{
   // Specify a clean-up masking operation to be applied to the raw 
   // images before any transformations are applied. Separate masks
   // are kept for the two members of the Couple. The caller retains
   // ownership of the passed mask object.
   //
   // Note that the CleanMask is different from ZeroMask. ZeroMask
   // gets applied after the transformations take place, and matches
   // the final maps in pixel dimensions and axis extents. Just one
   // instance of ZeroMask is the kept for both members.
 
   if (fCleanMask[p]) {
      delete fCleanMask[p];
      fCleanMask[p] = 0;
      deletemap(p);
   }
   if (mask == 0)
      return;
   TFile *fresults = open(fResultsPath[p]);
   if (fresults == 0) {
      printf("Error in Couple::setclean - no raw images found!\n");
      return;
   }
   TH2D *h2d = (TH2D*)fresults->Get("hmu");
   if (h2d == 0) {
      printf("Error in Couple::setclean - no hmu image found!\n");
      return;
   }
   fCleanMask[p] = new Map2D(*h2d);
   fCleanMask[p]->SetDirectory(0);
   fCleanMask[p]->Reset();
   fCleanMask[p]->Fill(mask);
   deletemap(p);
}

Map2D *Couple::getmask()
{
   // Return a copy of the zero mask, if any.
   // The caller owns the returned object.

   if (fZeroMask)
      return new Map2D(*fZeroMask);
   else
      return 0; 
}

Map2D *Couple::getclean(int p)
{
   // Return a copy of the clean mask, if any, for member p.
   // The caller owns the returned object.

   if (fCleanMask[p])
      return new Map2D(*fCleanMask[p]);
   else
      return 0; 
}

Map2D *Couple::getmap(int p, const char *name)
{
   // Return a copy of the named image from the first (p=0) or second (p=1)
   // member of this Couple, if any. If argument name is defaulted, whatever
   // name was last requested is returned. The caller owns the returned object.

   if (name == 0) {
      if (fMap[p] == 0 && fResultsName[p].Length() == 0) {
         printf("Error in Couple::getmap - map name is missing!\n");
         return 0;
      }
      name = fResultsName[p].Data();
   }
   else if (fResultsName[p] != name) {
      deletemap(p);
   }
   if (fMap[p] == 0) {
      TFile *fresults = open(fResultsPath[p]);
      if (fresults == 0)
         return 0;
      TH2D *h2d = (TH2D*)fresults->Get(name);
      if (h2d == 0)
         return 0;
      Map2D m2d(*h2d);
      if (fCleanMask[p])
         m2d.Mask(fCleanMask[p]);
      fResultsName[p] = name;
      TH2D hmap("hmap", "", fPixels[0], fXrange[0], fXrange[1],
                            fPixels[1], fYrange[0], fYrange[1]);
      fMap[p] = new Map2D(hmap);
      fMap[p]->SetDirectory(0);
      fMap[p]->SetName(h2d->GetName());
      fMap[p]->SetTitle(h2d->GetTitle());
      fMap[p]->Fill(&m2d);
      if (fabs(fXYskew[p]) > 0.01) {
         printf("Warning in Couple::getmap - plotting with skew "
                "is not currently supported, ignoring x/y skew=%lf\n",
                fXYskew[p]);
      }
      fMap[p]->Rescale(fXscale[p], fYscale[p], 1);
      fMap[p]->Rotate(0, 0, fPsirot[p]);
      fMap[p]->Shift(fXshift[p], fYshift[p], 0);
      if (fZeroMask)
         fMap[p]->Mask(fZeroMask);
      fMap[p]->GetXaxis()->SetTitle(h2d->GetTitle());
      fMap[p]->GetXaxis()->SetTitle(h2d->GetXaxis()->GetTitle());
      fMap[p]->GetYaxis()->SetTitle(h2d->GetYaxis()->GetTitle());
      fMap[p]->GetYaxis()->SetTitleOffset(h2d->GetYaxis()->GetTitleOffset());
      fMap[p]->SetContour(100);
      fMap[p]->Shift(0,0,-fMap[p]->GetMinimum());
   }
   return new Map2D(*fMap[p]);
}

void Couple::lineup()
{
   // Performs a dynamic visual alignment image 1 relative to image 0
   // in the couple, adjusting the internal parameters to keep track of
   // the net transformation leading to the final alignment. It works
   // using the latest generated images from each of the two members.
   // Typically one performs a draw(), followed by a rough set of shifts
   // + rotations + inversions to get the two images close to alignment,
   // then finishes with a call to lineup().
 
   if (fMap[0] == 0 || fMap[1] == 0) {
      printf("Error in Couple::lineup - please update the two graphics "
             "displays with images you would like to align, "
             "then try again.\n");
      return;
   }
   printf(
   "Performing a dynamic visual alignment of two images. At the prompt,\n"
   "you can enter the following commands to tweak the alignment, where\n"
   "the factor M is an optional multiplier of the minimum change step:\n"
   " translation or rotation:\n"
   " [M]r : shift image 2 to the right\n"
   " [M]l : shift image 2 to the left\n"
   " [M]u : shift image 2 up\n"
   " [M]d : shift image 2 down\n"
   " [M]p : rotate image 2 counter-clockwise\n"
   " [M]n : rotate image 2 clockwise\n"
   " [M]i : zoom in (expand) image 2\n"
   " [M]o : zoom out (contract) image 2\n"
   " [M]s : expand x, contract y on image 2\n"
   " [M]t : contract x, expand y on image 2\n"
   " [M]b : make image 2 more bright\n"
   " [M]f : make image 2 more faint\n"
   " [M]B : make image 1 more bright\n"
   " [M]F : make image 1 more faint\n"
   " [M]q : quit\n"
   "    * : (or anything else) flip between the two images again\n"
   );
   TExec *pal1 = new TExec("pal1", "Couple::select_palette(0);");
   TExec *pal2 = new TExec("pal2", "Couple::select_palette(1);");
   while (true) {
      TCanvas *c1 = select_canvas(0);
      pal1->Draw();
      fMap[0]->Draw("colz");
      c1->Update();
      usleep(5e5);
      pal2->Draw();
      fMap[1]->Draw("col same");
      c1->Update();
      pal1->Draw();
      printf("r/l/u/d/p/n/i/o/s/t/b/f/B/F/q/*? ");
      std::cout << std::flush;
      std::string res;
      std::getline(std::cin, res);
      if (res.size() == 0)
         continue;
      char cmd = res.back();
      res.pop_back();
      int mult = 1;
      if (res.size() > 0)
         mult = std::atoi(res.c_str());
      if (cmd == 'r') {
         double dx = mult * fMap[1]->GetXaxis()->GetBinWidth(1);
         fMap[1]->Shift(dx, 0, 0);
         fXshift[1] += dx;
      }
      else if (cmd == 'l') {
         double dx = mult * fMap[1]->GetXaxis()->GetBinWidth(1);
         fMap[1]->Shift(-dx, 0, 0);
         fXshift[1] -= dx;
      }
      else if (cmd == 'u') {
         double dy = mult * fMap[1]->GetYaxis()->GetBinWidth(1);
         fMap[1]->Shift(0, dy, 0);
         fYshift[1] += dy;
      }
      else if (cmd == 'd') {
         double dy = mult * fMap[1]->GetYaxis()->GetBinWidth(1);
         fMap[1]->Shift(0, -dy, 0);
         fYshift[1] -= dy;
      }
      else if (cmd == 'p') {
         double dphi = mult * 2.0 / fMap[1]->GetNbinsX();
         fMap[1]->Rotate(0, 0, dphi);
         fPsirot[1] += dphi;
         double sx = fXshift[1];
         double sy = fYshift[1];
         fXshift[1] = sx * cos(dphi) - sy * sin(dphi);
         fYshift[1] = sy * cos(dphi) + sx * sin(dphi);
      }
      else if (cmd == 'm') {
         double dphi = mult * 2.0 / fMap[1]->GetNbinsX();
         fMap[1]->Rotate(0, 0, -dphi);
         fPsirot[1] -= dphi;
         double sx = fXshift[1];
         double sy = fYshift[1];
         fXshift[1] = sx * cos(dphi) + sy * sin(dphi);
         fYshift[1] = sy * cos(dphi) - sx * sin(dphi);
      }
      else if (cmd == 'i') {
         double dsize = 1 + mult * 2.0 / fMap[1]->GetNbinsX();
         fMap[1]->Rescale(dsize, dsize, 1);
         fXscale[1] *= dsize;
         fYscale[1] *= dsize;
      }
      else if (cmd == 'o') {
         double dsize = 1 - mult * 2.0 / fMap[1]->GetNbinsX();
         fMap[1]->Rescale(dsize, dsize, 1);
         fXscale[1] *= dsize;
         fYscale[1] *= dsize;
      }

// A linear coordinate transformation is used to map from from raw images
// from the camera to the maps that are returned by the getmap method.
//
//    r' = S + R I r
//
// where R and I are 2x2 matrices and S is a 2-vector.
//
//      / cos_phi  -sin_phi \
//  R = |                     |
//      \ sin_phi   cos_phi /
//
//       / Ix     Ixy \
//  I = |              |
//       \ 0      Iy  /
//
//  S = [ Sx  Sy ]
//
// The task in this code section is to update the current R, I, and S
// in such a way as to accumulate many small shifts, rotations, etc.
// to describe their net effect in a single final set R, I, and S.
// For displacement and rotation tweaks this is straightforward, as
// coded above.
//
// shifts dS: Sx -> Sx + dSx, Sy -> Sy + dSy
//
// tilts d_phi: phi -> phi + d_phi, 
//              Sx -> Sx cos_dphi - Sy sin_dphi
//              Sy -> Sy cos_dphi + Sx sin_dphi
//
// But for tweaks that affect the strain matrix I (stretch/squeeze), some work 
// is needed to see how these aggregate together into accumulated transform
// matrices. Consider a stretch by factor s along x, 1/s along y described by
//
//      / s   0  \
// D = |          |
//      \ 0  1/s /
//
// The task is to see how (S + D R I) gets rewritten as (S + R' I') such that
// the effect of D is incorporated into updates R -> R' and I -> I', further
// stipulating that I' must remain upper-diagonal as required in a minimal
// description of a general linear coordinate transform. The solution to this
// linear algebra problem is:
//
// R' = R At, where At is the transpose of an orthogonal matrix A, and
// I' = A Rt D R I
//
// where the new matrix A is introduced for the sole purpose of returning I'
// to upper-diagonal form.
//
//      / 1 -a \
// A = |        |
//      \ a  1 /
// 
// The following lines work out the solution for I' and R' to leading order
// in the small tweak parameter s.
//
//           / s cos_phi^2 + 1/s sin_phi^2      cos_phi sin_phi (1/s - s) \
// Rt D R = |                                                              |
//           \ cos_phi sin_phi (1/s - s)   s sin_phi^2 + (1/s) cos_phi^2  /
//
//                                       / cos_2phi  -sin_2phi \
//        = 1/2 [(s + 1/s) + (s - 1/s)  |                       | ]
//                                       \-sin_2phi  -cos_2phi /
//
//                        / cos_2phi  -sin_2phi \
//        ~= 1 + (s - 1) |                       | 
//                        \-sin_2phi  -cos_2phi /
//
//                          / Ix cos_2phi      Ixy cos_2phi - Iy sin_2phi \
// Rt D R I ~= I + (s - 1) |                                               |
//                          \ -Ix sin_2phi    -Ixy sin_2phi - Iy cos_2phi /
//
// This result is not ready to be adopted as I' because it is not in upper-
// diagonal form yet. For That, we need to find the right matrix A that can
// make this happen. Consider the general form,
//
//      / 1  -a \   / u   v \      / u - a w   v - a x \
//     |         | |         |  = |                     |
//      \ a   1 /   \ w   x /      \ w + a u   x + a v /
//
// To make the rhs upper-diagonal, we need to let
//
//  a = -w/u = (s - 1) sin_2phi / [1 + (s - 1) cos_2phi]
//   ~= (s - 1) sin_2phi 
//
// to leading order in the small parameter s-1. Putting all of this together,
//
//                    / Ix cos_2phi   Ixy cos_2phi - 2 Iy sin_2phi \
// I' ~= I + (s - 1) |                                              |
//                    \    0                          -Iy cos_2phi /
//
//       / cos_phi  -sin_phi \   / 1 -a \      / cos_phi'   -sin_phi' \
// R' = |                     | |        |  = |                        |
//       \ sin_phi   cos_phi /   \ a  1 /      \ sin_phi'    cos_phi' /
//
// where phi' = phi + (s - 1) sin_2phi, and
//
// S' = [s Sx, (1/s) Sy ]

      else if (cmd == 's') {
         double dsize = 1 + mult * 2.0 / fMap[1]->GetNbinsX();
         fMap[1]->Rescale(dsize, 1/dsize, 1);
         fXshift[1] *= dsize;
         fYshift[1] /= dsize;
         double cos2phi = cos(2 * fPsirot[1]);
         double sin2phi = sin(2 * fPsirot[1]);
         fXscale[1] *= 1 + (dsize - 1) * cos2phi;
         fYscale[1] *= 1 - (dsize - 1) * cos2phi;
         fXYskew[1] *= 1 + (dsize - 1) * cos2phi;
         fXYskew[1] -= 2 * (dsize - 1) * sin2phi * fYscale[1];
         fPsirot[1] += (dsize - 1) * sin2phi;
      }
      else if (cmd == 't') {
         double dsize = 1 - mult * 2.0 / fMap[1]->GetNbinsX();
         fMap[1]->Rescale(dsize, 1/dsize, 1);
         fXshift[1] *= dsize;
         fYshift[1] /= dsize;
         double cos2phi = cos(2 * fPsirot[1]);
         double sin2phi = sin(2 * fPsirot[1]);
         fXscale[1] *= 1 + (dsize - 1) * cos2phi;
         fYscale[1] *= 1 - (dsize - 1) * cos2phi;
         fXYskew[1] *= 1 + (dsize - 1) * cos2phi;
         fXYskew[1] -= 2 * (dsize - 1) * sin2phi * fYscale[1];
         fPsirot[1] += (dsize - 1) * sin2phi;
      }
      else if (cmd == 'b') {
         double zmax = fMap[1]->GetMaximum();
         fMap[1]->SetMaximum(zmax * (1 - 0.01 * mult));
      }
      else if (cmd == 'f') {
         double zmax = fMap[1]->GetMaximum();
         fMap[1]->SetMaximum(zmax * (1 + 0.01 * mult));
      }
      else if (cmd == 'B') {
         double zmax = fMap[0]->GetMaximum();
         fMap[0]->SetMaximum(zmax * (1 - 0.01 * mult));
      }
      else if (cmd == 'F') {
         double zmax = fMap[0]->GetMaximum();
         fMap[0]->SetMaximum(zmax * (1 + 0.01 * mult));
      }
      else if (cmd == 'q') {
         break;
      }
   }
}

void Couple::select_palette(int p)
{
   if (p == 0)
      gStyle->SetPalette(kBird, 0, 1.0);
   else
      gStyle->SetPalette(kBird, 0, 0.5);
}

Map2D *Couple::cleanup(int p)
{
   // Examine the pixels of the rocking curve topographs from the first
   // member (p=0) or second member (p=1) of this couple, and distinguish
   // the ones that show a legitimate diffraction peak from those that do
   // not. Zero any pixels that do not have a legitimate diffraction peak
   // and return a pointer to the image as a Map2D. The caller becomes
   // the owner of the returned object.
   //
   // This algorithm has been tested on the images taken during run
   // cls-6-2019. It may need to be tweaked to work properly on images
   // taken from other runs.

   TFile *fresults = open(fResultsPath[p]);
   if (fresults == 0)
      return 0;
   TH2D *hbase = (TH2D*)fresults->Get("hbase");
   if (hbase == 0)
      return 0;
   TH2D *hamp = (TH2D*)fresults->Get("hamp");
   if (hamp == 0)
      return 0;
   hamp->Add(hbase);
   int col1 = 1;
   int colN = hamp->GetNbinsX();
   int row1 = 1;
   int rowN = hamp->GetNbinsY();
   TH1D *hedge[4];
   hedge[0] = hamp->ProjectionX("col1", 1, rowN);
   if (hedge[0]->Integral() == 0) {
      printf("Error in Couple::cleanup - map is empty!\n");
      return 0;
   }
   while (hedge[0]->GetBinContent(col1) == 0)
      ++col1;
   while (hedge[0]->GetBinContent(colN) == 0)
      --colN;
   hedge[1] = hamp->ProjectionY("row1", 1, colN);
   while (hedge[1]->GetBinContent(row1) == 0)
      ++row1;
   while (hedge[1]->GetBinContent(rowN) == 0)
      --rowN;
   hedge[0] = hamp->ProjectionX("row1", row1+1, row1+1);
   hedge[1] = hamp->ProjectionY("col1", col1+1, col1+1);
   hedge[2] = hamp->ProjectionX("rowN", rowN-1, rowN-1);
   hedge[3] = hamp->ProjectionY("colN", colN-1, colN-1);
   double bgheight = 9999;
   for (int side=0; side<4; ++side) {
      double height = hedge[side]->Integral() / hedge[side]->GetNbinsX();
      if (height < bgheight)
         bgheight = height;
      delete hedge[side];
   }
   Map2D *silh = new Map2D(*hamp);
   Map2D silh0(*silh);
   int speck0 = silh0.Despeckle();
   silh->ZeroSuppress(bgheight);
   Map2D silh1(*silh);
   int speck1 = silh1.Despeckle();
   double bgheight_step[2];
   if (bgheight > 100) {
      bgheight_step[0] = 1;
      bgheight_step[1] = 0.1;
   }
   else {
      bgheight_step[0] = 0.1;
      bgheight_step[1] = 0.05;
   }
   while (speck1 < speck0 * 1.2) {
      bgheight += bgheight_step[0];
      silh->ZeroSuppress(bgheight);
      Map2D silh2(*silh);
      speck1 = silh2.Despeckle();
   }
   int speck2 = speck1;
   while (speck2 >= speck1) {
      bgheight += bgheight_step[1];
      silh->ZeroSuppress(bgheight);
      speck1 = speck2;
      Map2D silh3(*silh);
      speck2 = silh3.Despeckle();
   }
   silh->Despeckle(10000);
   silh->SetDirectory(0);
   select_canvas(p);
   silh->Draw("colz");
   setclean(p, silh);
   deletemap(p);
   return silh;
}

void Couple::polycrop(int p, const char *name)
{
   // Display the topograph image with the given name belonging
   // to the first (p=0) or second (p=1) member of this couple
   // as a topological color map and ask the user to select the
   // vertices of a polygon to use for cropping the extents of
   // the image. All pixels outside the selected polygon are set
   // to matte (zero).
   //
   // WARNING: Only the copy of the named map in memory is modified
   // by this method. If you want to use the adopt the result to
   // update the mask for this couple, take the object returned
   // by getmap(p,name) after this call, mask it with *getmask()
   // if any, and register the result using setmask().

   if (name == 0 && fMap[p] == 0) {
      printf("Error in Couple::polycrop - no image found!\n");
      return;
   }
   else if (getmap(p, name) == 0) {
      printf("Error in Couple::polycrop - image not available!\n");
      return;
   }
   TCanvas *can = select_canvas(p);
   fPolyVertex.clear();
   can->DeleteExec("polypicker");
   if (p == 0)
      can->AddExec("polypicker", "Couple::polypicker(0);");
   else 
      can->AddExec("polypicker", "Couple::polypicker(1);");
   fMap[p]->SetStats(0);
   fMap[p]->Draw("colz");
   can->Update();
   printf("Use the mouse to pick the vertices of the cropping polygon\n"
          "in window c%d, then click the middle mouse button to terminate.\n"
          "The polygon will close itself by connecting the last vertex\n"
          "to the first.\n", p+1);
   std::cout << "waiting..." << std::flush;
   fPicking = this;
}

void Couple::polypicker(int p)
{
   // This callback is called whenever a mouse event (move, click, etc.)
   // happens inside one of the open ROOT graphics windows. This one
   // activates when the user left-clicks inside the active region of a
   // displayed TH2D histogram any number of times to select a polygonal
   // region of interest, then right-clicks in the image to execute the
   // cropping action.

   TCanvas *can;
   if (p == 0)
      can = (TCanvas*)gROOT->FindObject("c1");
   else
      can = (TCanvas*)gROOT->FindObject("c2");
   TObject *select = can->GetSelected();
   if (!select) {
      return;
   }
   if (!select->InheritsFrom("Map2D")) {
      return;
   }
   Map2D *refmap = (Map2D*)select;
   int event = can->GetEvent();
   if (event == 1) {
      double x = can->AbsPixeltoX(can->GetEventX());
      double y = can->AbsPixeltoY(can->GetEventY());
      int px = refmap->GetXaxis()->FindBin(x);
      int py = refmap->GetYaxis()->FindBin(y);
      fPicking->fPolyVertex.push_back(std::pair<int,int>(px,py));
      printf("cropping polygon vertex %ld (%d,%d)\n",
             fPicking->fPolyVertex.size(), px, py);
   }
   else if (event == 12) {
      std::cout << "polygon closed, cropping..." << std::flush;
      can->DeleteExec("polypicker");
      fPicking->polycropper(p);
      std::cout << "done." << std::endl;
      fPicking = 0;
   }
}

void Couple::polycropper(int p)
{
   // Clear all pixels in the topo map that lie outside the cropping region
   // defined as the interior of a polygon with vertices stored in vector
   // fPolyVertex, that have just been selected by the user with the mouse.
   // This algorithm assumes that the sides of the polygon are touching only
   // at the vertices.

   TCanvas *can;
   if (p == 0)
      can = (TCanvas*)gROOT->FindObject("c1");
   else
      can = (TCanvas*)gROOT->FindObject("c2");
   TObject *select = can->GetSelected();
   TH2D *refmap = (TH2D*)select;
   Map2D mask(*refmap);
   mask.Flood(1);
   int nv = fPolyVertex.size();
   int nborder = 0;
   for (int iv=0; iv<nv; ++iv) {
      int i = fPolyVertex[iv].first;
      int j = fPolyVertex[iv].second;
      mask.SetBinContent(i,j,1);
      int ie = fPolyVertex[(iv+1)%nv].first;
      int je = fPolyVertex[(iv+1)%nv].second;
      mask.SetBinContent(ie,je,0);
      nborder += 1;
      int dx = (ie > i)? +1 : -1;
      int dy = (je > j)? +1 : -1;
      double x = (i + 0.5) * dx;
      double y = (j + 0.5) * dy;
      double slope = abs((je - j)/(ie - i + 1e-9));
      while (mask.GetBinContent(i,j) > 0) {
         mask.SetBinContent(i,j,0);
         nborder += 1;
         double xb = floor(x + 1) - x;
         double yb = floor(y + 1) - y;
         if (yb > xb * slope) {
            x += xb;
            y += xb * slope;
            i += dx;
         }
         else if (yb < xb * slope) {
            x += yb / slope;
            y += yb;
            j += dy;
         }
         else {
            x += xb;
            y += yb;
            i += dx;
            j += dy;
         }
      }
   }
   int nx = mask.GetNbinsX();
   int ny = mask.GetNbinsY();
   mask.Despeckle((nx*ny - nborder - 1) / 2);
   if (mask.GetBinContent(1,1) > 0 || mask.GetBinContent(nx,1) > 0 ||
                                      mask.GetBinContent(1,ny) > 0 ||
                                      mask.GetBinContent(nx,ny) > 0)
   {
      Map2D mask2(mask);
      mask2.Flood(2).Mask(&mask,1);
      mask.Add(&mask2);
      mask.Shift(0,0,-1);
   }
   for (int p=0; p<2; ++p) {
      if (refmap == fMap[p]) {
         fMap[p]->Mask(&mask);
         fMap[p]->Draw("colz");
         can->Update();
         if (fZeroMask)
            mask.Mask(fZeroMask);
         setmask(&mask);
      }
   }
}

Map2D *Couple::zerocurl(const Map2D *gx, const Map2D *gy, Map2D **silhouette)
{
   // Run Map2D::TotalCurl(gx,gy) but then take the result and
   // apply a smooth correction to the resulting contour path
   // integral that makes it sum to zero around the closed loop.
   // The caller becomes owner of the returned object.
 
   Map2D *curlmap = new Map2D(*gx);
   curlmap->TotalCurl(gx, gy, silhouette);

   // TotalCurl generates some 1D histograms as by-products,
   // use them to find a correction that zeros the total curl.

   TH1D *loopcurl = (TH1D*)gROOT->FindObject("loopcurl");
   TH1I *looppixi = (TH1I*)gROOT->FindObject("looppixi");
   TH1I *looppixj = (TH1I*)gROOT->FindObject("looppixj");
   int nsteps = loopcurl->GetNbinsX();
   int ctotal = loopcurl->GetBinContent(nsteps);
   for (int istep=1; istep <= nsteps; ++istep) {
      int i = looppixi->GetBinContent(istep);
      int j = looppixj->GetBinContent(istep);
      double ci = loopcurl->GetBinContent(istep);
      curlmap->SetBinContent(i+1, j+1, ci - ctotal * i/nsteps);
   }
   return curlmap;
}

Double_t Couple::uncurl(const Map2D *gx, const Map2D *gy, Map2D **silhouette)
{
   // Rotate the gradient vector (gx,gy) by angle chi, then
   // evaluate Map2D::TotalCurl(gx,gy) and allow the user to
   // adjust chi to zero out the curl.

   double curl=0;
   double dchi=0;
   Map2D *mcurl;
   Map2D *grad[2];
   mcurl = new Map2D(*gx);
   grad[0] = new Map2D(*gx);
   grad[1] = new Map2D(*gx);
   Map2D gxm(*gx), gym(*gy);
   gxm.Mask(&gym);
   gym.Mask(&gxm);
   while (true) {
      double chi[2];
      for (int p=0; p<2; ++p) {
         chi[p] = dchi * M_PI/180 + fChirot[p];
         grad[p]->Reset();
      }
      printf("chi[0], chi[1] = %lf, %lf\n", chi[0], chi[1]);
      grad[0]->Add(&gxm,-cos(chi[0]));
      grad[0]->Add(&gym,-sin(chi[0]));
      grad[1]->Add(&gym,-cos(chi[1]));
      grad[1]->Add(&gxm,+sin(chi[1]));
      curl = mcurl->TotalCurl(grad[0],grad[1],silhouette);
      TH1D *loopcurl = (TH1D*)gROOT->FindObject("loopcurl");
      TCanvas *c1 = (TCanvas*)gROOT->FindObject("c1");
      if (c1) {
         c1->cd();
         loopcurl->Draw();
         c1->Update();
      }
      printf("chi=%lf, curl=%lf, new chi (degrees, or enter to accept)? ",
             dchi, curl);
      std::cout << std::flush;
      std::string res;
      std::getline(std::cin, res);
      if (res.size() == 0) {
         fChirot[0] = chi[0];
         fChirot[1] = chi[1];
         break;
      }
      dchi = std::atof(res.c_str());
   }
   mcurl->SetName("mcurl");
   grad[0]->SetName("gradx");
   grad[1]->SetName("grady");
   mcurl->SetStats(0);
   grad[0]->SetStats(0);
   grad[1]->SetStats(0);
   return curl;
}

Map2D *Couple::divergence(const Map2D *gx, const Map2D *gy)
{
   // Computes the 2d divergence of the vector field (gx,gy)
   // where the gradient is evaluated in a basis that is
   // rotated by chi with respect to the camera.
 
   Map2D *grad[2];
   grad[0] = new Map2D(*gx);
   grad[1] = new Map2D(*gx);
   Map2D gxm(*gx), gym(*gy);
   gxm.Mask(&gym);
   gym.Mask(&gxm);
   grad[0]->Reset();
   grad[0]->Add(&gxm,-cos(fChirot[0]));
   grad[0]->Add(&gym,-sin(fChirot[0]));
   grad[1]->Reset();
   grad[1]->Add(&gym,-cos(fChirot[1]));
   grad[1]->Add(&gxm,+sin(fChirot[1]));
   grad[0]->SetName("gradx");
   grad[1]->SetName("grady");
   grad[0]->SetStats(0);
   grad[1]->SetStats(0);
   Map2D *mdiv = new Map2D(*gx);
   mdiv->Divergence(&gxm, &gym);
   mdiv->SetName("mdiv");
   mdiv->SetDirectory(0);
   mdiv->SetStats(0);
   return mdiv;
}
 
TFile *Couple::open(const char *filename)
{
   if (strlen(filename) == 0)
      return 0;
   int n;
   TString path;
   for (n=0; n<nsearchpaths; ++n) {
      path = TString(searchpath[n]) + TString(filename);
      ifstream fin(path.Data());
      if (fin.good())
         break;
   }
   if (n < nsearchpaths) {
      TDirectory *cwd(gDirectory);
      TFile *file = TFile::Open(path.Data());
      cwd->cd();
      return file;
   }
   printf("Error in Couple::open - unable to open input file %s\n", filename);
   return 0;
}

void Couple::dump()
{
   printf("Contents of Couple %s: %s\n", GetName(), GetTitle());
   printf("   images are %d x %d pixels\n", fPixels[1], fPixels[0]);
   printf("              [%lf, %lf] mm horizontal\n", fXrange[0], fXrange[1]);
   printf("              [%lf, %lf] mm vertical\n", fYrange[0], fYrange[1]);
   if (fZeroMask)
      printf("   mask is enabled (%s)\n", fZeroMask->GetName());
   else
      printf("   mask is disabled\n");
   for (int p=0; p<2; ++p) {
      if (fResultsPath[p].Length() == 0)
         printf("   couple member %d: (none)\n", p+1);
      else {
         printf("   couple member %d: %s ", p+1, fResultsPath[p].Data());
         if (fResultsName[p].Length() == 0)
            printf("\n");
         else
            printf("(%s)\n", fResultsName[p].Data());
         printf("      cleaning mask ");
         if (fCleanMask[p])
            printf("%s (%s)\n", fCleanMask[p]->GetName(),
                                fCleanMask[p]->GetTitle());
         else
            printf("(none)\n");
         printf("      shifted (%lf, %lf) mm\n", fXshift[p], fYshift[p]);
         printf("      rotated %lf degrees\n", fPsirot[p] * 180/M_PI);
         printf("      rescaled by (%lf, %lf)\n", fXscale[p], fYscale[p]);
         printf("      x/y skew factor %lf\n", fXYskew[p]);
         printf("      chi angle %lf degrees\n", fChirot[p] * 180/M_PI);
         printf("      phi angle %lf degrees\n", fPhirot[p] * 180/M_PI);
         printf("      current image: ");
         if (fMap[p])
            printf("%s (%s)\n", fMap[p]->GetName(), fMap[p]->GetTitle());
         else
            printf("(none)\n");
      }
   }
}

void Couple::draw(const char *name)
{
   // Draw the two images of the Couple with the given name, or the
   // last requested images if name is not given.

   TCanvas *c1 = select_canvas(0);
   if (getmap(0, name))
      fMap[0]->Draw("colz");
   TCanvas *c2 = select_canvas(1);
   if (getmap(1, name))
      fMap[1]->Draw("colz");
}

TCanvas *Couple::select_canvas(int p)
{
   // Select one of two square drawing areas for subsequent drawing.

   gStyle->SetCanvasPreferGL(true);
   TCanvas *can;
   if (p == 0) {
      can = (TCanvas*)gROOT->FindObject("c1");
      if (can == 0) {
         can = new TCanvas("c1", "c1", 50, 50, 560, 500);
         can->SetRightMargin(0.15);
      }
   }
   else {
      can = (TCanvas*)gROOT->FindObject("c2");
      if (can == 0) {
         can = new TCanvas("c2", "c2", 1300, 50, 560, 500);
         can->SetRightMargin(0.15);
      }
   }
   can->cd();
   return can;
}