#ifndef __IMAGE_ALIGN__
#define __IMAGE_ALIGN__

#include <fftw3.h> 
#include "image2d.h"
#include "matrix.h"
#include <set>
#include <map>
//#define PMFSZ 64
#define PMFSZ 16

class imagealign
{
 public:

  struct searchspace
  {
    double xofflow, xoffhigh; ///< xoffset of imB normalized to [0,1]
    double yofflow, yoffhigh; ///< yoffset of imB normalized to [0,1]
  };

  struct transform
  {
    double xoff, yoff;
  };                              

  struct points
  {
    double x, y;
  };

  imagealign();
  ~imagealign();
  /// initialize
  /**
   * @param A image A
   * @param B image B
   * @param xoff xoffset to apply to image B
   */
  int init( image2d &A, image2d &B, int xoff );
  int setid( char Astack, int Aslice, char Bstack, int Bslice );
  int setoffset( int xoff );
  /// calculate mutual information score
  double miscore();
#ifndef _NOMAGICK_
  int write();
#endif
  int search( int min, int max, double *s );

  /// multi resultion search
  int searchmultires( transform &bt );

#ifndef _NOMAGICK_
  int write( image2d &A, image2d &B,
             transform ttt,
             Magick::Image &rim);
#endif

  double miscore( image2d &A, image2d &B,
                  image2d &mapa, image2d &mapb, 
                  int xoff, int yoff );

  double miscore( image2d &A, image2d &B,
                  int xoff, int yoff );

 private:
  /// calculate the pmf of an image
  /**
   * calculates pmf of box 
   * [0, swidth], [im.height, send]
   * @param hist histrogram to store
   * @param im image to calculate
   */
  int calcpmf( double *hist, image2d &im );
  int calcpmf( double *hist, image2d &im, image2d &map );

  /// calculates joint pmf with respect to offset
  int calcjpmf(double *hist, image2d &A,
               image2d &B);
  int calcjpmf(double *hist, image2d &A,
               image2d &B, image2d &map);

  double linearfit( points *pts, int npts, double *s, int degree );
  double jacobi( matrix &a, matrix &X, matrix &B );

  int search( image2d &A, image2d &B, 
              std::vector< transform > &mysp,
              std::vector< transform > &newsp,
              transform &btrans );
    
  image2d imA, imB;

  image2d imAres[4];
  image2d imBres[4];

  char Ast, Bst;
  int Asl, Bsl;
  double mi;
  int offset;
  std::string alignname;
};

#endif
