#include "imagealign.h"

struct tcmp
{
  bool operator()( const imagealign::transform &a, 
                   const imagealign::transform &b ) const
  {
    if( a.xoff < b.xoff )
      return 1;
    else if( a.xoff == b.xoff )
      {
        if( a.yoff < b.yoff )
          return 1;
      }
    return 0;
  }
};

imagealign::imagealign()
{
}

imagealign::~imagealign()
{
}

int imagealign::init( image2d &A, image2d &B,
                      int xoff )
{
  offset = xoff;
  imA.load( &A );
  imB.load( &B );
  //  imA.normalizeminmax();
  //  imB.normalizeminmax();
  return 0;
}

int imagealign::setid( char Astack, int Aslice, char Bstack, int Bslice )
{
  Ast = Astack;
  Asl = Aslice;
  Bst = Bstack;
  Bsl = Bslice;
  char buf[512];
  sprintf( buf, "%c%i%c%i_%i", Ast, Asl, Bst, Bsl, offset );
  alignname = buf;
  return 0;
}

int imagealign::setoffset( int xoff )
{
  offset = xoff;
  char buf[512];
  sprintf( buf, "%c%i%c%i_%i", Ast, Asl, Bst, Bsl, offset );
  alignname = buf;
  return 0;
}

int imagealign::searchmultires( transform &bt )
{
  std::cout << "loading low res...\n";
  imAres[0].loadlowres( &imA, 2 );
  imBres[0].loadlowres( &imB, 2 );
  std::cout << "0 ";
  for( int i = 1; i < 4; i++ )
    {
      imAres[i].loadlowres( &(imAres[i-1]), 2 );
      imBres[i].loadlowres( &(imBres[i-1]), 2 );
      std::cout << i << " ";
    }
  std::cout << "\n";
  std::vector< transform > mysp;
  std::vector< transform > mynewsp;

  //int lh = imA.getheight();
  //int lw = imA.getwidth();

  int lh = imAres[2].getheight();
  int lw = imAres[2].getwidth();

  transform ttt;
  for( double i = 0; i < lh; i = i + 1 )
    {
      for( double j = 0; j < lw; j = j + 1 )
        {
          //ttt.xoff = i/( (double)lh );
          ttt.xoff = 0;
          ttt.yoff = j/( (double)lw );
          mysp.push_back( ttt );
          ttt.xoff = -ttt.xoff;
          mysp.push_back( ttt );
          ttt.yoff = -ttt.yoff;
          mysp.push_back( ttt );
          ttt.xoff = -ttt.xoff;
          mysp.push_back( ttt );
        }
    }
  
  for( int i = 2; i >= 0; i-- )
    {
      std::cout << "seearching level: " << i
                << " : " << mysp.size() << " "
                << imAres[i].getheight() 
                << ":" 
                << imBres[i].getwidth() << "\n";
      search( imAres[i], imBres[i], mysp, mynewsp, ttt );
      mysp = mynewsp;
    }
  
  search( imA, imB, mysp, mynewsp, ttt );
  std::cout << "best transform: " << ttt.xoff*imA.getheight() << " " 
            << ttt.yoff*imA.getwidth() << "\n";
  bt = ttt;
  bt.xoff = bt.xoff * imA.getheight();
  bt.yoff = bt.yoff * imA.getwidth();
  return 0;
}

int imagealign::search( image2d &A, image2d &B, 
                        std::vector< transform > &mysp,
                        std::vector< transform > &newsp,
                        transform &btrans )
{
  std::map< transform, double, tcmp >::iterator st;
  std::map< transform, double, tcmp > ssp, snsp;
  for( unsigned i = 0; i < mysp.size(); i++ )
    ssp[ mysp[i] ] = 0;
  int h, w;
  h = A.getheight();
  w = A.getwidth();
  st = ssp.begin();
  double avgsc, stdvar, max, count;
  count = 0;
  avgsc = stdvar = max = 0;
  for( ; st != ssp.end(); st++ )
    {
      double s;
      if( st->first.yoff < 0)
        s = 0;
      else
        s = miscore( A, B, (int)(st->first.xoff*h), (int)(st->first.yoff*w) );
      if( s > 0 )
        count = count + 1;
      std::cout << "calculating... "
                << h << " " << w << " "
                << st->first.xoff << ":" 
                << st->first.yoff << " "
                << st->first.xoff*h << "::"
                << st->first.yoff*w << " "
                << s << "\n";;
      st->second = s;
      avgsc += st->second;
    }
  std::cout << "number of alignments: " << count << " avg: ";
  avgsc = avgsc/ ( count );
  std::cout << avgsc << "\n";
  st = ssp.begin();
  for( ; st != ssp.end(); st++ )
    {
      if( st->second > 0 )
        
        stdvar = stdvar + pow( avgsc - st->second, 2 );
      if( st->second > max )
        {
          max = st->second;
          btrans = st->first;
        }
    }
  stdvar = sqrt( stdvar/count );
  std::cout << "max: " << max << " ";
  max = max - stdvar;
  std::cout << "cutoff: " << max << "\n";
  double dtx, dty;
  dtx = 1.0/( (double)h*2.0 );
  dty = 1.0/( (double)w*2.0 );
  st = ssp.begin();
  for( ; st != ssp.end(); st++ )
    {
      if( st->second > max )
        {
          snsp[ st->first ] = 0;
          transform ttt;

          ttt.yoff = st->first.yoff;
          ttt.xoff = st->first.xoff+dtx;
          snsp[ttt] = 0;

          ttt.yoff = st->first.yoff+dty;
          ttt.xoff = st->first.xoff+dtx;
          snsp[ttt] = 0;

          ttt.yoff = st->first.yoff-dty;
          ttt.xoff = st->first.xoff+dtx;
          snsp[ttt] = 0;

          ttt.yoff = st->first.yoff-dty;
          ttt.xoff = st->first.xoff;
          snsp[ttt] = 0;

          ttt.yoff = st->first.yoff+dty;
          ttt.xoff = st->first.xoff;
          snsp[ttt] = 0;

          ttt.yoff = st->first.yoff-dty;
          ttt.xoff = st->first.xoff-dtx;
          snsp[ttt] = 0;

          ttt.yoff = st->first.yoff;
          ttt.xoff = st->first.xoff-dtx;
          snsp[ttt] = 0;

          ttt.yoff = st->first.yoff+dty;
          ttt.xoff = st->first.xoff-dtx;
          snsp[ttt] = 0;
        }
    }
  st = snsp.begin();
  for( ; st != snsp.end(); st++ )
    newsp.push_back( st->first );
  return 0;
}

double imagealign::miscore( image2d &A, image2d &B, 
                            int xoff, int yoff )
{
  double score = 0;

  score = 0;
  double Ha, Hb;
  Ha = Hb = 0;
  double ahist[PMFSZ], bhist[PMFSZ];
  double jhist[PMFSZ*PMFSZ];
  image2d imAcut, imBcut;
  subimage simg;

  // top left and bottom right corners for image A
  int Atlx, Atly;
  int Btlx, Btly;

  if( xoff >= 0 )
    {
      Atlx = xoff;
      Btlx = 0;
    }
  else
    {
      Btlx = -xoff;
      Atlx = 0;
    }
  if( yoff > 0 )
    {
      Atly = yoff;
      Btly = 0;
    }
  else
    {
      Atly = 0;
      Btly = -yoff;
    }

  // max sizes with the top corner defined
  int Amaxw, Amaxh, Bmaxh, Bmaxw;
  Amaxw = A.getwidth() - Atly;
  Amaxh = A.getheight() - Atlx;

  Bmaxw = B.getwidth() - Btly;
  Bmaxh = B.getheight() - Btlx;
  int wwidth, wheight;
  if( Amaxw < Bmaxw )
    wwidth = Amaxw;
  else
    wwidth = Bmaxw;
  if( Amaxh < Bmaxh )
    wheight = Amaxh;
  else
    wheight = Bmaxh;

  // force 5% overlap min
  if( wheight * wwidth < .05 *A.getwidth()*A.getheight() )
    return 0;

  simg.init( A.getimageptr(), Atlx, Atlx+wheight, Atly, Atly+wwidth,
             A.getheight(), A.getwidth() );
  imAcut.load( &simg );

  simg.init( B.getimageptr(), Btlx, Btlx+wheight, Btly, Btly+wwidth,
             B.getheight(), B.getwidth() );
  imBcut.load( &simg );;

  calcpmf( ahist, imAcut );
  calcpmf( bhist, imBcut );
  calcjpmf( jhist, imAcut, imBcut );
  
  double ha, hb, hab;
  ha = hb = hab = 0;
  for( int i = 0; i < PMFSZ; i++ )
    {
      if( ahist[i] > 0 )
        ha = ha + ahist[i]*log(ahist[i]);
      if( bhist[i] > 0 )
        hb = hb + bhist[i]*log(bhist[i]);
    }
  ha = -ha;
  hb = -hb;  
  
  for( int i = 0; i < PMFSZ; i++ )
    for( int j = 0; j < PMFSZ; j++ )
      if( ahist[i] > 0 && bhist[j] > 0 )
        {
          score = score + jhist[i*PMFSZ+j] *
            log(jhist[i*PMFSZ+j]/(ahist[i]*bhist[j]));
          hab = hab + jhist[i*PMFSZ+j]*log(jhist[i*PMFSZ+j]);
        }

  score = (ha + hb ) / ( -hab );
  
  return score;
}

double imagealign::miscore( image2d &A, image2d &B, 
                            image2d &mapa, image2d &mapb, 
                            int xoff, int yoff )
{
  double score = 0;

  score = 0;
  double Ha, Hb;
  Ha = Hb = 0;
  double ahist[PMFSZ], bhist[PMFSZ];
  double jhist[PMFSZ*PMFSZ];
  image2d imAcut, imBcut;
  image2d mapAcut, mapBcut;
  subimage simg;

  // top left and bottom right corners for image A
  int Atlx, Atly;
  int Btlx, Btly;

  if( xoff >= 0 )
    {
      Atlx = xoff;
      Btlx = 0;
    }
  else
    {
      Btlx = -xoff;
      Atlx = 0;
    }
  if( yoff > 0 )
    {
      Atly = yoff;
      Btly = 0;
    }
  else
    {
      Atly = 0;
      Btly = -yoff;
    }

  // max sizes with the top corner defined
  int Amaxw, Amaxh, Bmaxh, Bmaxw;
  Amaxw = A.getwidth() - Atly;
  Amaxh = A.getheight() - Atlx;

  Bmaxw = B.getwidth() - Btly;
  Bmaxh = B.getheight() - Btlx;
  int wwidth, wheight;
  if( Amaxw < Bmaxw )
    wwidth = Amaxw;
  else
    wwidth = Bmaxw;
  if( Amaxh < Bmaxh )
    wheight = Amaxh;
  else
    wheight = Bmaxh;

  // force 5% overlap min
  if( wheight * wwidth < .05 *A.getwidth()*A.getheight() )
    return 0;

  simg.init( A.getimageptr(), Atlx, Atlx+wheight, Atly, Atly+wwidth,
             A.getheight(), A.getwidth() );
  imAcut.load( &simg );
  simg.init( mapa.getimageptr(), Atlx, Atlx+wheight, Atly, Atly+wwidth,
             A.getheight(), A.getwidth() );
  mapAcut.load( &simg );

  simg.init( B.getimageptr(), Btlx, Btlx+wheight, Btly, Btly+wwidth,
             B.getheight(), B.getwidth() );
  imBcut.load( &simg );
  simg.init( mapb.getimageptr(), Btlx, Btlx+wheight, Btly, Btly+wwidth,
             B.getheight(), B.getwidth() );
  mapBcut.load( &simg );

  int numpixels;
  numpixels = mapAcut.logicaland( mapBcut );
  if( numpixels < 26214 )
    return 0;

  calcpmf( ahist, imAcut, mapAcut );
  calcpmf( bhist, imBcut, mapAcut );
  calcjpmf( jhist, imAcut, imBcut, mapAcut );
  
  double ha, hb, hab;
  ha = hb = hab = 0;
  for( int i = 0; i < PMFSZ; i++ )
    {
      if( ahist[i] > 0 )
        ha = ha + ahist[i]*log(ahist[i]);
      if( bhist[i] > 0 )
        hb = hb + bhist[i]*log(bhist[i]);
    }
  ha = -ha;
  hb = -hb;  
  
  for( int i = 0; i < PMFSZ; i++ )
    for( int j = 0; j < PMFSZ; j++ )
      if( ahist[i] > 0 && bhist[j] > 0 )
        {
          score = score + jhist[i*PMFSZ+j] *
            log(jhist[i*PMFSZ+j]/(ahist[i]*bhist[j]));
          hab = hab + jhist[i*PMFSZ+j]*log(jhist[i*PMFSZ+j]);
        }

  score = (ha + hb ) / ( -hab );
  
  return score;
}

int imagealign::search( int min, int max, double *s )
{
  double mscore = 0;
  int goodoff = 0;
  //  std::cout << "[";
  for( int i = min; i < max; i = i + 1 )
    {
      setoffset( i );
      double m;
      m = miscore();
      s[ i-min] =  m;
      //      std::cout << "score: " << m << "\n";
    }
  double yint;
  yint = s[0];
  double slope;
  slope = ( s[max-min-1] - s[0] )/ ((double)(max-min));
  //  std::cout << "y=" << slope << "x + " << yint << "\n";
  for( int i = 0; i < max-min; i++ )
    {
      //      s[i] = s[i]- ( ((double)i) * slope - yint );
      if( s[i] < 0 )
        s[i] = 0;
      double m;
      m = s[i];
      if( m > mscore )
        {
          mscore = m;
          goodoff = i+min;
        }
    }
  /*
  std::cout << "]\n";
  std::cout << "best offset: " << goodoff 
            << " " << mscore << "\n";
  */
  if( goodoff > min + 10 && goodoff < max-10 )
    {
      double sum;
      sum = 0;
      for( int i = 0; i < max-min; i++ )
        sum += s[i];
      for( int i = 0; i < max-min; i++ )
        s[i] = s[i] / sum;
    }
  else
    for( int i = 0; i < max-min; i ++ )
      s[i] = 0;
  return goodoff;
}

double imagealign::linearfit( points *pts, int npts, double *a, int degree )
{
  double res = 0;
  matrix X, sx, tx, Y, A;
  sx.init( npts, degree+1 );
  Y.init( npts, 1 );
  for( int i = 0; i < npts; i++ )
    {
      Y.set( pts[i].y, i, 0 );
      sx.set( 1, i, 0 );
      for( int j = 1; j <= degree; j++ )
        sx.set( pow(pts[i].x, j), i, j );
    }
  tx = sx.transpose();
  X = tx * sx;
  Y = tx * Y;
  if( X.rows() != degree+1 || X.cols() != degree+1 ||
      Y.rows() != degree+1 || Y.cols() != 1 )
    {
      std::cout << "ERROR matrix size\n";
      std::cout << X.rows() << "x" << X.cols() << "\n";
      std::cout << Y.rows() << "x" << Y.cols() << "\n";
      exit(0);
    }
  A.init( degree+1, 1 );
  res = jacobi( X, A, Y );
  for( int i = 0; i <= degree; i++ )
    a[i] = A.get(i, 0);
  return res;
}

double imagealign::jacobi( matrix &A, matrix &X, matrix &B )
{
  double error;
  int flag = 1;
  double *x1, *x2, *xt;
  x1 = new double[X.rows()];
  x2 = new double[X.rows()];
  for( int i = 0; i < X.rows(); i++ )
    x1[i] = 0;
  while( flag < 20 )
    {
      error = 0;
      double t;
      for( int i = 0; i < A.rows(); i++ )
        {
          t = 0;
          for( int j = 0; j < A.cols(); j++ )
            t = t + A.get(i,j)*x1[j];
          t = (B.get(i,0) - t )/(A.get(i,i));
          error = error + fabs(t);
          x2[i] = x1[i] + t;
        }
      if( error < 0.0001 )
        flag++;
      xt = x1;
      x1 = x2;
      x2 = xt;
    }
  for( int i = 0; i < X.rows(); i++ )
    X.set( x1[i], i, 0 );
  delete [] x1;
  delete [] x2;
  return error;
}

double imagealign::miscore()
{
  mi = 0;
  double Ha, Hb;
  Ha = Hb = 0;
  double ahist[PMFSZ], bhist[PMFSZ];
  double jhist[PMFSZ*PMFSZ];
  image2d imAcut, imBcut;
  subimage simg;

  simg.init( imA.getimageptr(), 0, 512, offset, 512, 512, 512 );
  imAcut.load( &simg );
  //  imAcut.normalizeminmax();
  calcpmf( ahist, imAcut );

  simg.init( imB.getimageptr(), 0, 512, 0, 512-offset, 512, 512 );
  imBcut.load( &simg );
  //  imBcut.normalizeminmax();
  calcpmf( bhist, imBcut );

  calcjpmf( jhist, imAcut, imBcut );

  for( int i = 0; i < PMFSZ; i++ )
    for( int j = 0; j < PMFSZ; j++ )
      if( ahist[i] > 0 && bhist[j] > 0 )
        mi = mi + jhist[i*PMFSZ+j] *
          log(jhist[i*PMFSZ+j]/(ahist[i]*bhist[j]));
  for( int i = 0; i < PMFSZ; i++ )
    {
      if( ahist[i] > 0 )
        Ha = Ha - ahist[i]*log( ahist[i] );
      if( bhist[i] > 0 )
        Hb = Hb - bhist[i]*log( bhist[i] );
    }
  /*
  if( Ha < Hb )
    mi = mi + Ha;
  else
    mi = mi + Hb;
  */
  return mi;
}

int imagealign::calcpmf( double *hist, image2d &im )
{
  for( int i = 0; i < PMFSZ; i++ )
    hist[i] = 0;
  double ipart, fpart, sum;
  sum = 0;
  int w, h;
  h = im.getheight();
  w = im.getwidth();
  for( int i = 0; i < h; i++ )
    {
      for( int j = 0; j < w; j++ )
        {
          double v;
          v = im.getpixel(i,j);
          fpart = modf(v, &ipart );
          hist[(int)ipart] += (1-fpart);
          if( ipart + 1 < PMFSZ )
            hist[(int)ipart+1] += fpart;            
        }
    }
  for( int i = 0; i < PMFSZ; i++ )
    sum = sum + hist[i]; 
  for( int i = 0; i < PMFSZ; i++ )
    hist[i] = hist[i] / sum;
  return 0;
}

int imagealign::calcpmf( double *hist, image2d &im, image2d &map )
{
  for( int i = 0; i < PMFSZ; i++ )
    hist[i] = 0;
  double ipart, fpart, sum;
  sum = 0;
  int w, h;
  h = im.getheight();
  w = im.getwidth();
  for( int i = 0; i < h; i++ )
    {
      for( int j = 0; j < w; j++ )
        {
          double v;
          if( map.getpixel(i,j) > 0.5 )
            {
              v = im.getpixel(i,j);
              fpart = modf(v, &ipart );
              hist[(int)ipart] += (1-fpart);
              if( ipart + 1 < PMFSZ )
                hist[(int)ipart+1] += fpart;
            }
        }
    }
  for( int i = 0; i < PMFSZ; i++ )
    sum = sum + hist[i]; 
  for( int i = 0; i < PMFSZ; i++ )
    hist[i] = hist[i] / sum;
  return 0;
}

int imagealign::calcjpmf(double *hist, image2d &A, image2d &B)
{
  double sum, ipart1, ipart2, fpart1, fpart2;
  sum = 0;
  int h, w;
  h = A.getheight();
  w = A.getwidth();
  if( h != B.getheight() || w != B.getwidth() )
    {
      std::cout << "ERROR size in joint pmf calculation failed\n";
      exit(0);
    }
  for( int i = 0; i < PMFSZ*PMFSZ; i++ )
    hist[i] = 0;
  for( int i = 0; i < h; i++ )
    {
      for( int j = 0; j < w; j++ )
        {
          double v1, v2;
          v1 = A.getpixel(i, j);
          v2 = B.getpixel(i, j);
          fpart1 = modf(v1, &ipart1 );
          fpart2 = modf(v2, &ipart2 );
          hist[(int)ipart1*PMFSZ+(int)ipart2] += (1.0-fpart1)*(1.0-fpart2);
          if( ipart2 + 1 < PMFSZ )
            hist[(int)ipart1*PMFSZ+(int)ipart2+1] += (1.0-fpart1)*(fpart2);
          if( ipart1 + 1 < PMFSZ )
            hist[(int)(ipart1+1)*PMFSZ+(int)ipart2] += (fpart1)*(1.0-fpart2);
          if( ipart2 + 1 < PMFSZ && ipart1 + 1 < PMFSZ )
            hist[(int)(ipart1+1)*PMFSZ+(int)ipart2+1] += (fpart1)*(fpart2);
        }
    }
  for( int i = 0; i < PMFSZ*PMFSZ; i++ )
    sum = sum + hist[i];
  for( int i = 0; i < PMFSZ*PMFSZ; i++ )
    hist[i] = hist[i] / sum;
  return 0;
}

int imagealign::calcjpmf(double *hist, image2d &A, image2d &B, image2d &map)
{
  double sum, ipart1, ipart2, fpart1, fpart2;
  sum = 0;
  int h, w;
  h = A.getheight();
  w = A.getwidth();
  if( h != B.getheight() || w != B.getwidth() )
    {
      std::cout << "ERROR size in joint pmf calculation failed\n";
      exit(0);
    }
  for( int i = 0; i < PMFSZ*PMFSZ; i++ )
    hist[i] = 0;
  for( int i = 0; i < h; i++ )
    {
      for( int j = 0; j < w; j++ )
        {
          if( map.getpixel(i,j) > 0.5 )
            {
              double v1, v2;
              v1 = A.getpixel(i, j);
              v2 = B.getpixel(i, j);
              fpart1 = modf(v1, &ipart1 );
              fpart2 = modf(v2, &ipart2 );
              hist[(int)ipart1*PMFSZ+(int)ipart2] += (1.0-fpart1)*(1.0-fpart2);
              if( ipart2 + 1 < PMFSZ )
                hist[(int)ipart1*PMFSZ+(int)ipart2+1] += (1.0-fpart1)*(fpart2);
              if( ipart1 + 1 < PMFSZ )
                hist[(int)(ipart1+1)*PMFSZ+(int)ipart2] += (fpart1)*(1.0-fpart2);
              if( ipart2 + 1 < PMFSZ && ipart1 + 1 < PMFSZ )
                hist[(int)(ipart1+1)*PMFSZ+(int)ipart2+1] += (fpart1)*(fpart2);
            }
        }
    }
  for( int i = 0; i < PMFSZ*PMFSZ; i++ )
    sum = sum + hist[i];
  for( int i = 0; i < PMFSZ*PMFSZ; i++ )
    hist[i] = hist[i] / sum;
  return 0;
}

#ifndef _NOMAGICK_
int imagealign::write( image2d &A, image2d &B,
                       transform ttt,
                       Magick::Image &rim)
{
  Magick::Geometry geo;
  int oimw, oimh;
  oimh = A.getheight()+B.getheight();
  oimw = A.getwidth()+B.getwidth();
  geo.height( A.getheight()+B.getheight() );
  geo.width( A.getwidth()+B.getwidth() );
  Magick::Image tim(geo, Magick::Color("black") );
  rim = tim;
  rim.type( Magick::TrueColorType );
  transform atrans, btrans;  // transforms for image a and b
  if( ttt.xoff >= 0 )
    {
      atrans.xoff = 0;
      btrans.xoff = ttt.xoff;
    }
  else
    {
      atrans.xoff = -ttt.xoff;
      btrans.xoff = 0;
    }
  if( ttt.yoff >= 0 )
    {
      atrans.yoff = 0;
      btrans.yoff = ttt.yoff;
    }
  else
    {
      atrans.yoff = -ttt.yoff;
      btrans.yoff = 0;
    }
  /*
  std::cout << A.getheight() << " " << B.getheight() << ":"
            << A.getwidth() << " " << B.getwidth() << "::"
            << oimh << " " << oimw << "\n";
  std::cout << atrans.xoff << ":" << atrans.yoff << " "
            << btrans.xoff << ":" << btrans.yoff << "\n";
  */
  Magick::PixelPacket *pc;
  pc = rim.getPixels(0, 0, oimw, oimh );
  Magick::ColorRGB rgb(0,0,0);
  for( int i = 0; i < oimh; i++ )
    {
      for( int j = 0; j < oimw; j++ )
        {
          int ti, tj;
          ti = (int)(i - atrans.xoff);
          tj = (int)(j - atrans.yoff);
          rgb.blue(0);
          rgb.green(0);
          rgb.red(0);
          if( ti >= 0 && ti < A.getheight() &&
              tj >= 0 && tj < A.getwidth() )
            {
              rgb.red( A.getpixel( ti, tj ) );
              //              std::cout << "A";
            }

          ti = (int)(i - btrans.xoff);
          tj = (int)(j - btrans.yoff);
          if( ti >= 0 && ti < B.getheight() &&
              tj >= 0 && tj < B.getwidth() )
            {
              rgb.blue( B.getpixel( ti, tj ) );
              //              std::cout << "B";
            }
          pc[i*oimw+j] = rgb;
        }
    }
  //std::cout << "\n";
  rim.syncPixels();
  return 0;
}


int imagealign::write()
{
  Magick::Geometry geo;
  geo.height( imA.getheight() );
  geo.width( imA.getwidth()+offset );
  int oh, ow;
  oh = geo.height();
  ow = geo.width();
  Magick::Image im(geo, Magick::Color("black") );
  im.type(Magick::TrueColorType);
  Magick::ColorRGB rgb(0,0,0);
  Magick::PixelPacket *pc;
  pc = im.getPixels(0,0,ow,oh);
  for( int i = 0; i < imA.getheight(); i++ )
    {
      for( int j = 0; j < imA.getwidth(); j++)
        {
          rgb = pc[i*ow+j];
          rgb.red( imA.getpixel(i, j) );
          pc[i*ow+j] = rgb;
          rgb = pc[i*ow+j+offset];
          rgb.blue( imB.getpixel(i, j) );
          pc[i*ow+j+offset] = rgb;
        }
    }
  im.write( alignname+".jpg" );
  return 0;
}
#endif
