#include "alignimage.h"

void *blendthread_start( void *d )
{
  alignimage_blend_td *atd;
  atd = (alignimage_blend_td *) d;
  blend_calctodo( atd );
  pthread_barrier_wait( atd->sync );
  blend_mainloop( atd );
  pthread_barrier_wait( atd->sync );
  blend_cleanup( atd );
  return NULL;
}

int blend_calctodo( alignimage_blend_td *atd )
{
  int sz;
  sz =  (int)atd->cmaps->size();
  atd->todo = new std::vector< std::pair<int, int> > [ sz ];
  atd->imh = atd->overlap->getheight();
  atd->imw = atd->overlap->getwidth();
  for( int i = 0; i < sz; i++ )
    {
      if( atd->imh != (*(atd->cmaps))[i]->getheight() ||
          atd->imw != (*(atd->cmaps))[i]->getwidth() )
        {
          std::cout << "ERROR im size not match\n";
          exit(0);
        }
    }
  atd->sz = sz;
  double *map;
  map = atd->overlap->getimageptr();
  atd->dptrs.clear();
  for( int i = 0; i < sz; i++ )
    {
      atd->dptrs.push_back( (*(atd->cmaps))[i]->getimageptr() );
      atd->lowlockidx.push_back(0);
      atd->highlockidx.push_back(0);
    }

  for( int i = atd->lowrow; i < atd->highrow; i++ )
    {
      for( int j = 0; j < atd->imw; j++ )
        {
          if( map[ i*atd->imw + j ] > 0.5 )
            {
              std::pair< int, int > pr;
              pr.first = i;
              pr.second = j;
              for( int k = 0; k < sz; k++ )
                if( atd->dptrs[k][i*atd->imw + j] > 0.5 )
                  atd->todo[k].push_back( pr );
            }
        }
      if( i == atd->lowrow+1 )
        for( int k = 0; k < sz; k++ )
          atd->lowlockidx[k] = atd->todo[k].size();
      if( i == atd->highrow-2 )
        for( int k = 0; k < sz; k++ )
          atd->highlockidx[k] = atd->todo[k].size();
    }

  pthread_mutex_lock( atd->clock );
  std::cout << "thread: " << (unsigned int)atd->tid << "\n";
  std::cout << "\t" << atd->lowrow << " " << atd->highrow << "\n";
  for( int k = 0; k < sz; k++ )
    {
      if( atd->lowlockidx[k] == 0 &&  atd->highlockidx[k] == 0 )
        atd->lowlock[k] = atd->highlock[k] = NULL;
      std::cout << "\t" << atd->todo[k].size() << " "
                << atd->lowlockidx[k] << " " << atd->highlockidx[k] 
                << " " << (unsigned long int)atd->lowlock[k]
                << " " << (unsigned long int)atd->highlock[k] << "\n";
    }
  pthread_mutex_unlock( atd->clock );

  return 0;
}

int blend_cleanup( alignimage_blend_td *atd )
{
  delete [] atd->todo;
  return 0;
}

int blend_mainloop( alignimage_blend_td *atd )
{
  std::vector< std::pair<int, int> > *todo;
  todo = atd->todo;
  for( int itr = 0; itr < 1000; itr++ )
    {
      double err;
      err = 0;
      for( int k = 0; k < atd->sz; k++ )
        {
          if( atd->lowlock[k] != NULL )
            {
              pthread_mutex_lock( atd->lowlock[k] );
            }
          for( unsigned i = 0; i < todo[k].size(); i++ )
            {
              
              if( i == (unsigned int)atd->lowlockidx[k] )
                {
                  if( atd->lowlock[k] != NULL )
                    {
                      pthread_mutex_unlock( atd->lowlock[k] );
                    }
                }
              if( i == (unsigned int)atd->highlockidx[k] )
                {
                  if( atd->highlock[k] != NULL )
                    {
                      pthread_mutex_lock( atd->highlock[k] );
                    }
                }
              
              double value;
              int x, y;
              x = todo[k][i].first;
              y = todo[k][i].second;
              value = blend_calc( x, y, atd->dptrs[k], atd->imh, atd->imw );
              if( atd->dptrs[k][x*atd->imw+y] > 0.000001 )
                {
                  if( value > 0.00001 )
                    {
                      err = err + fabs(atd->dptrs[k][x*atd->imw+y] - value );
                      atd->dptrs[k][x*atd->imw+y] = value;
                    }
                }
            }
          if( todo[k].size() == (unsigned int)atd->lowlockidx[k] )
            {
              if( atd->lowlock[k] != NULL )
                {
                  pthread_mutex_unlock( atd->lowlock[k] );
                }
            }
          if( atd->highlock[k] != NULL )
            {
              pthread_mutex_unlock( atd->highlock[k] );
            }
        }
      if( itr % 50 == 0 )
        {
          std::cout << "thread: " << (unsigned int)atd->tid 
                    << " " << itr << " " << err << "\n";
          pthread_barrier_wait( atd->sync );
        }
    }
  return 0;
}

double blend_calc( int x, int y, double *d, int imh, int imw )
{
  int i = x-3;
  int j = y-3;
  if( i < 0 )
    i = 0;
  if( j < 0 )
    j = 0;
  int k = x+4;
  int l = y+4;
  if( k > imh )
    k = imh;
  if( l > imw )
    l = imw;
  double z;
  z = 0;
  double v;
  v = 0;
  for( ; i < k; i++ )
    {
      for(int jj = j ; jj < l; jj++ )
        {
          z += 1;
          v = v + d[i*imw+jj];
        }
    }
  if( z == 0 )
    return 0;
  return v/z;
}

alignimage::alignimage()
{
}

alignimage::~alignimage()
{
  del();
}

int alignimage::del()
{
  for( unsigned i = 0; i < imsets.size(); i++ )
    {
      if( imsets[i] != NULL )
        {
          for( std::map<int, alignedset* >::iterator it = imsets[i]->sets.begin();
               it != imsets[i]->sets.end(); it++)
            {
              delete it->second;
            }
          delete imsets[i];
        }
    }
  imsets.clear();
  return 0;
}

int alignimage::init( dbalnimage *d )
{
  db = d;
  return 0;
}

int alignimage::deliteration( int iter )
{
  std::vector< int > myaids;
  db->get_align_sets( iter, myaids );
  for( unsigned int i = 0; i < myaids.size(); i++ )
    {
      std::cout << "deleting stack: " << myaids[i] << "\n";
      delaid( myaids[i] );
    }
  db->deleteconnection( iter );
  return 0;
}

int alignimage::delaid( int daid )
{
  return db->deleteaid( daid );
}

int alignimage::load()
{
  del();
  int flag = 1;
  int iter = 0;
  std::vector< int > myids;
  while( flag == 1 )
    {
      myids.clear();
      db->get_align_sets( iter, myids );
      if( myids.size() == 0 )
        flag = 0;
      else
        {
          imsets.push_back( new imageset );
          imsets[iter]->level = iter;
          for( unsigned i = 0; i < myids.size(); i++ )
            {
              if( imsets[iter]->sets.count( myids[i] ) > 0 )
                {
                  std::cout << "duplicated aid: " << myids[i] << "\n";
                  for( unsigned j = 0; j < myids.size(); j++ )
                    {
                      std::cout << myids[j] << " ";
                    }
                  std::cout << "\n";
                  exit(0);
                }
              tile mytile;
              mytile.init( myids[i], db );
              imsets[iter]->sets[ myids[i] ] = new alignedset;
              imsets[iter]->sets[ myids[i] ]->aid = myids[i];
              mytile.getsize( imsets[iter]->sets[ myids[i] ]->sizex,
                              imsets[iter]->sets[ myids[i] ]->sizey );
              // load data describing which tiles were used
              std::set< tileinfo, tileinfo > subtiles;
              int ttiter;
              std::vector< int > ttids, ttxtrans, ttytrans;
              ttids.clear();
              ttxtrans.clear();
              ttytrans.clear();
              db->get_alignedset_data( myids[i], ttiter, ttids,
                                       ttxtrans, ttytrans );
              std::cout << "aid: " << myids[i] << " " << ttiter;
              for( unsigned int z = 0; z < ttids.size(); z++ )
                {
                  if( ttiter > 0 )
                    {
                      imageset *tas = imsets[ttiter-1];
                      std::set< tileinfo, tileinfo >::iterator iter;
                      for( iter = tas->sets[ ttids[z] ]->tiles.begin();
                           iter != tas->sets[ ttids[z] ]->tiles.end(); iter++ )
                        {
                          tileinfo ttf = *iter;
                          ttf.transx = ttf.transx + ttxtrans[z];
                          ttf.transy = ttf.transy + ttytrans[z];
                          subtiles.insert( ttf );
                        }
                    }
                  else
                    {
                      tileinfo ttf;
                      db->decodeto2d( ttids[z], ttf.x, ttf.y );
                      ttf.transx = ttxtrans[z];
                      ttf.transy = ttytrans[z];
                      subtiles.insert( ttf );
                    }
                }
              std::set< tileinfo, tileinfo >::iterator myiter;
              for( myiter = subtiles.begin(); myiter != subtiles.end(); myiter++ )
                {
                  std::cout << " (" << myiter->x << ", " << myiter->y << " : "
                            << myiter->transx << " " << myiter->transy << ")";
                }
              std::cout << "\n";
              imsets[iter]->sets[ myids[i] ]->tiles = subtiles;
            }
          // get image links
          std::vector< dbalnimage::image_conn > connections;
          db->get_imconnection( iter, connections );
          setconn sc;
          for( unsigned i = 0; i < connections.size(); i++ )
            {
              sc.aid = connections[i].b;
              sc.con.lowtx = connections[i].lowx;
              sc.con.lowty = connections[i].lowy;
              sc.con.hightx = connections[i].highx;
              sc.con.highty = connections[i].highy;
              imsets[iter]->sets[ connections[i].a ]->connections.push_back( sc );
            }
          iter++;
        }
    }
  dotout();
  return 0;
}

int alignimage::showbox( int aid, Magick::Image &im, dbimage *imdb, Magick::Image &noblend )
{
  alignedset *myas;
  myas = NULL;
  for( unsigned int i = 0; i < imsets.size() && myas == NULL; i++ )
    {
      if(imsets[i]->sets.count( aid ) > 0 )
        {
          myas = imsets[i]->sets[aid];
        }
    }
  if( myas == NULL )
    {
      std::cout << "did not fin anyting\n";
      std::cout << "aid = " << aid << "\n";
      return -1;
    }
  std::cout << "aid: " << myas->aid << " (" << myas->sizex << " "
            << myas->sizey << ")\n";
  Magick::Geometry geo( myas->sizey, myas->sizex );
  im.type( Magick::TrueColorType );
  im.size( geo );
  noblend.type( Magick::TrueColorType );
  noblend.size( geo );
  Magick::PixelPacket *pixel_cache = noblend.getPixels(0,0,myas->sizey, myas->sizex);
  Magick::ColorRGB rgb(0, 0, 0);
  std::list<Magick::Drawable> drawlist;
  std::set<tileinfo, tileinfo>::iterator titer;
  int minx, miny;
  minx = 10000;
  miny = 10000;
  for( titer = myas->tiles.begin(); titer != myas->tiles.end(); titer++ )
    {
      if( minx > titer->transx )
        minx = titer->transx;
      if( miny > titer->transy )
        miny = titer->transy;
    }
  for( titer = myas->tiles.begin(); titer != myas->tiles.end(); titer++ )
    {
      int tx, ty;
      tx = titer->transx - minx;
      ty = titer->transy - miny;
      std::cout << "drawing box at " << ty << " " << tx << "\n";
      drawlist.push_back( Magick::DrawableRectangle( ty, tx, ty+512, tx+512 ) );
      image2d chr, chb;
      imdb->getimage( titer->x, titer->y, 10, 3, chr );
      imdb->getimage( titer->x, titer->y, 10, 2, chb );
      chr.normalize( 0, 256.0/4096.0 );
      chb.normalize( 0, 256.0/4096.0 );
      for( int i = 0; i < 512; i++ )
        {
          for( int j = 0; j < 512; j++) 
            {
              double val;
              val = chr.getpixel( i, j );
              rgb.red( val );
              val = chb.getpixel( i, j );
              rgb.blue( val );
              int ii, jj;
              ii = i + tx;
              jj = j + ty;
              pixel_cache[ii*myas->sizey+jj] = rgb;
            }
        }
    }

  
  noblend.syncPixels();
  noblend.write( "noblend.png" );

  im.strokeColor("red"); // Outline color
  Magick::ColorRGB alpha;
  alpha.alpha(0.75);
  alpha.green( 1 );
  alpha.blue( 0 );
  alpha.red( 0 );
  im.fillColor(alpha); 
  im.strokeWidth(5);
  im.draw(drawlist);
  return 0;
}

int alignimage::mergesets( std::vector< tile* > &tiles,
                           std::vector< setconn2 > &conn,
                           int iteration )
{
  int ntiles = (int) tiles.size();
  std::vector<int> nodeflag;
  std::vector<int> nodedegree;
  nodeflag.resize( ntiles, 0 );
  nodedegree.resize( ntiles, 0 );
  std::map< std::pair<int, int>, transrange > edges;
  std::cout << "copying edges\n";
  for( unsigned i = 0; i < conn.size(); i++ )
    {
      std::pair<int, int> t;
      t.first = conn[i].a;
      t.second = conn[i].b;
      edges[t] = conn[i].con;
      nodedegree[ t.second ] += 1;
    }
  int ncolors;
  ncolors = colorsegments( nodeflag, nodedegree, edges );
  std::cout << "found " << ncolors << " segments\n";
  // for each connected section of the graph create a new image
  for( int c = 1; c < ncolors; c++ )
    {
      std::cout << "processing color " << c << "\n";
      std::vector< std::pair<tile*, transrange> > ts;
      ts.clear();
      int start;
      start = -1;
      for( unsigned i = 0; i < nodeflag.size(); i++ )
        {
          if( nodeflag[i] == c )
            {
              if( nodedegree[i] == 0 )
                {
                  if( start == -1 )
                    start = i;
                  else
                    {
                      std::cout << "error multiple starts found\n";
                      exit(0);
                    }
                }
            }
        }
      if( start == -1 ) 
        {
          std::cout << "error no start found\n";
          exit(0);
        }
      std::cout << "calculating transformations\n";
      caltransforms( start, tiles, edges, ts );
      int aid;     // database align id

      // create data structure and initial structures in database
      std::cout << "creating db entry\n";
      aid = createalignset( iteration, ts );
      std::cout << "aid: " << aid << "\n";
      std::cout << "adding images to set\n";
      addtoimageset( aid, ts, imsets[iteration]->sets[aid]->sizex,
                     imsets[iteration]->sets[aid]->sizey );
    }
  imsets[iteration]->level = iteration;
  // create links between the images
  std::cout << "creating links " << iteration << "\n";
  createsetlinks( iteration );
  savelinks( iteration );
  dotout();
  return 0;
}

int alignimage::createsets( std::vector< tile* > &tiles,
                            std::vector< setconn2 > &conn )
{
  int ntiles = (int) tiles.size();
  std::vector<int> nodeflag;
  std::vector<int> nodedegree;
  nodeflag.resize( ntiles, 0 );
  nodedegree.resize( ntiles, 0 );
  std::map< std::pair<int, int>, transrange > edges;
  std::cout << "copying edges\n";
  for( unsigned i = 0; i < conn.size(); i++ )
    {
      std::pair<int, int> t;
      t.first = conn[i].a;
      t.second = conn[i].b;
      edges[t] = conn[i].con;
      nodedegree[ t.second ] += 1;
    }
  int ncolors;
  ncolors = colorsegments( nodeflag, nodedegree, edges );
  std::cout << "found " << ncolors << " segments\n";
  // for each connected section of the graph create a new image
  for( int c = 1; c < ncolors; c++ )
    {
      std::cout << "processing color " << c << "\n";
      std::vector< std::pair<tile*, transrange> > ts;
      ts.clear();
      int start;
      start = -1;
      for( unsigned i = 0; i < nodeflag.size(); i++ )
        {
          if( nodeflag[i] == c )
            {
              if( nodedegree[i] == 0 )
                {
                  if( start == -1 )
                    start = i;
                  else
                    {
                      std::cout << "error multiple starts found\n";
                      exit(0);
                    }
                }
            }
        }
      if( start == -1 ) 
        {
          std::cout << "error no start found\n";
          exit(0);
        }
      std::cout << "calculating transformations\n";
      caltransforms( start, tiles, edges, ts );
      int aid;     // database align id

      // create data structure and initial structures in database
      std::cout << "creating db entry\n";
      aid = createalignset( 0, ts );
      std::cout << "aid: " << aid << "\n";
      std::cout << "adding images to set\n";
      addtoimageset( aid, ts, imsets[0]->sets[aid]->sizex,
                     imsets[0]->sets[aid]->sizey );
    }
  imsets[0]->level = 0;
  // create links between the images
  std::cout << "creating links\n";
  createsetlinks( 0 );
  savelinks( 0 );
  dotout();
  return 0;
}

int alignimage::getrowalnjobs( std::vector< alignjob > &joblist,
                               int iteration)
{
  joblist.clear();
  int mt;
  if( iteration == -1 )
    mt = imsets.size()-1;
  else
    mt = iteration;
  alignjob aj;
  if( mt > -1 )
    {
      std::map<int, alignedset* >::iterator it;
      for( it = imsets[mt]->sets.begin();
           it != imsets[mt]->sets.end(); it++ )
        {
          aj.aid1 = it->second->aid;
          for( unsigned i = 0; i < it->second->connections.size(); i++ )
            {
              aj.aid2 = it->second->connections[i].aid;
              aj.lowx = (int) it->second->connections[i].con.lowtx;
              aj.lowy = (int) it->second->connections[i].con.lowty;
              aj.highx = (int) it->second->connections[i].con.hightx;
              aj.highy = (int) it->second->connections[i].con.highty;
              joblist.push_back( aj );
            }
        }
    }
  return 0;
}

int alignimage::savelinks( int iteration )
{
  std::vector< dbalnimage::image_conn > data;
  dbalnimage::image_conn td;
  std::map<int, alignedset* >::iterator it;
  for(it = imsets[iteration]->sets.begin(); it != imsets[iteration]->sets.end(); it++)
    {
      td.a = it->first;
      for( unsigned i = 0; i < it->second->connections.size(); i++ )
        {
          td.b = it->second->connections[i].aid;
          td.lowx = (int) it->second->connections[i].con.lowtx;
          td.lowy = (int) it->second->connections[i].con.lowty;
          td.highx = (int) it->second->connections[i].con.hightx;
          td.highy = (int) it->second->connections[i].con.highty;
          data.push_back( td );
        }
    }
  db->set_imconnection( iteration, data );
  return 0;
}

int alignimage::createsetlinks( int iteration )
{
  std::vector< int > aids;
  std::map<int, alignedset* >::iterator it;
  for( it = imsets[iteration]->sets.begin(); it != imsets[iteration]->sets.end(); it++ )
    aids.push_back( it->first );
  for( unsigned int i = 0; i < aids.size(); i++ )
    {
      for( unsigned int j = 0; j < aids.size(); j++ )
        {
          if( i != j )
            {
              //check if two aligned sets of images are next to each other and create the propoer connections
              createlinks(iteration, aids[i], aids[j] );
            }
        }
    }
  return 0;
}

int alignimage::createlinks( int iteration, int aid1, int aid2 )
{
  int flag;
  flag = 0;
  std::set< tileinfo >::iterator it1, it2;
  for( it1 = imsets[iteration]->sets[aid1]->tiles.begin(); 
       it1 != imsets[iteration]->sets[aid1]->tiles.end(); it1++ )
    {
      for( it2 = imsets[iteration]->sets[aid2]->tiles.begin(); 
           it2 != imsets[iteration]->sets[aid2]->tiles.end(); it2++ )
        {
          int rowdif;
          rowdif = it2->x - it1->y;
          if( abs( rowdif )  == 1 ) //on adjcent row
            {
              std::cout << "image " << aid1 << " and " << aid2 << " row adjcent\n";
              transrange tr;
              // aid2 above aid1
              if( rowdif < 0 )
                {
                  tr.lowtx = -1*round( imsets[iteration]->sets[aid2]->sizex * 0.9 );
                  tr.hightx = -1*round( imsets[iteration]->sets[aid2]->sizex * 0.1 );
                }
              // aid1 above aid2
              else
                {
                  tr.lowtx = round( imsets[iteration]->sets[aid1]->sizex * 0.1 );
                  tr.hightx = round( imsets[iteration]->sets[aid1]->sizex * 0.9 );
                }
              tr.lowty = -1*round( imsets[iteration]->sets[aid2]->sizey * 0.6 );
              tr.highty = round( imsets[iteration]->sets[aid1]->sizey * 0.6 );
              addconnection( imsets[iteration]->sets[aid1], aid2, tr );
            }
          else if( rowdif == 0 )     //on same row
            {
              int coldif;
              coldif = it2->y - it1->y;
              if( abs(coldif) == 1 )
                {
                  std::cout << "image " << aid1 << " and " << aid2 << " col adjcent\n";
                  transrange tr;
                  if( coldif < 0 )
                    {
                      tr.lowty = -1*round(imsets[iteration]->sets[aid2]->sizey * 0.9 );
                      tr.highty = -1*round(imsets[iteration]->sets[aid2]->sizey*0.1);
                    }
                  else
                    {
                      tr.lowty = round(imsets[iteration]->sets[aid1]->sizey*0.1 );
                      tr.highty = round(imsets[iteration]->sets[aid1]->sizey*0.9);
                    }
                  tr.lowtx = -32;
                  tr.hightx = 32;
                  addconnection( imsets[iteration]->sets[aid1], aid2, tr );
                }
            }
        }
    }
  return 0;
}

int alignimage::addconnection( alignedset *aset, int aid2, transrange tr )
{
  int flag = -1;
  for( unsigned i = 0; i < aset->connections.size(); i++ )
    {
      if( aset->connections[i].aid == aid2 )
        flag = i;
    }
  if( flag == -1 )
    {
      setconn sc;
      sc.aid = aid2;
      sc.con = tr;
      aset->connections.push_back( sc );
      return 0;
    }

  if( aset->connections[flag].con.lowtx > tr.lowtx )
    aset->connections[flag].con.lowtx = tr.lowtx;
  if( aset->connections[flag].con.lowty > tr.lowty )
    aset->connections[flag].con.lowty = tr.lowty;

  if( aset->connections[flag].con.hightx < tr.hightx )
    aset->connections[flag].con.hightx = tr.hightx;
  if( aset->connections[flag].con.highty < tr.highty )
    aset->connections[flag].con.highty = tr.highty;

  return 0;
}

int alignimage::addtoimageset( int aid, std::vector< std::pair<tile*, transrange> > &ts,
                               int &imh, int &imw)
{
  std::vector< int > channels, slices;
  ts[0].first->getids( slices, channels );
  std::cout << "number of channels: " << channels.size()
            << " number of slices: " << slices.size() << "\n";
  image2d immap, imoverlap;
  std::vector< image2d* > imcontribute;
  std::pair<int, double> maxtx;
  std::pair<int, double> maxty;
  std::cout << "fixing transformations\n";
  fixtransforms( ts, maxtx, maxty );
  std::cout << "preparing images\n";
  preparemerge( ts, immap, imoverlap, maxtx, maxty, imcontribute );
  std::cout << "setting image mask\n";
  db->set_aligned_mask( aid, &immap );
#ifndef _NOMAGICK_
  for( unsigned i = 0; i < imcontribute.size(); i++ )
    {
      Magick::Image mim;
      imcontribute[i]->getimage( mim );
      char buffer[128];
      sprintf( buffer, "contribute%i_%i.png", aid, i );
      mim.write( buffer );
    }
#endif
  for( unsigned int i = 0; i < channels.size(); i++ )
    for( unsigned int j = 0; j < slices.size(); j++ )
      createimage( aid, slices[j], channels[i], ts,
                   imcontribute );
  db->cleanlog();
  imh = imcontribute[0]->getheight();
  imw = imcontribute[0]->getwidth();
  for( unsigned i = 0; i < imcontribute.size(); i++ )
    delete imcontribute[i];
  imcontribute.clear();
  return 0;
}

int alignimage::preparemerge( std::vector< std::pair<tile*, transrange> > &ts,
                              image2d &immap, image2d &imoverlap,
                              std::pair<int, double> &maxtx,
                              std::pair<int, double> &maxty,
                              std::vector< image2d* > &imcontribute )
{
  immap.create( (int)maxtx.second, (int)maxty.second );
  imoverlap.create( (int)maxtx.second, (int)maxty.second );
  immap.clear();
  imoverlap.clear();

  for( unsigned i = 0; i < ts.size(); i++ )
    {
      image2d *nim;
      nim = new image2d;
      nim->load( &immap );
      imcontribute.push_back( nim );
    }
  // mark where each subimage goes and the overall reagion
  for( unsigned int i = 0; i < ts.size(); i++ )
    {
      image2d origmap;
      std::cout << "marking... " << i << "\n";
      ts[i].first->getmap( &origmap );
      addimage_mark( origmap, *imcontribute[i], immap, imoverlap, ts[i].second,
                     ts[i].first->getheight(), ts[i].first->getwidth() );
    }

#ifndef _NOMAGICK_
  Magick::Image mim;
  immap.getimage( mim );
  char buff[256];
  sprintf(buff, "immap.png" );
  mim.write( buff );
  imoverlap.getimage( mim );
  
  sprintf(buff, "imoverlap.png" );
  mim.write( buff );
#endif
  if( imoverlap.sum() > 0.5 )
    addimag_calcblend( imcontribute, immap, imoverlap );
  return 0;
}

int alignimage::fixtransforms( std::vector< std::pair<tile*, transrange> > &ts,
                               std::pair<int, double> &maxtx,
                               std::pair<int, double> &maxty )
{
  image2d im;
  // calculate the size of image
  std::pair< int, double > mintx, minty;
  mintx.first = minty.first = maxtx.first = maxty.first = 0;
  for( unsigned int i = 0; i < ts.size(); i++ )
    {
      if( ts[i].second.lowtx < mintx.second )
        {
          mintx.second = ts[i].second.lowtx;
          mintx.first = i;
        }
      if( ts[i].second.lowty < minty.second )
        {
          minty.second = ts[i].second.lowty;
          minty.first = i;
        }
    }
  // rewrite all transformations such that the min x, y is 0,0
  for( unsigned int i = 0; i < ts.size(); i++ ) 
    {
      ts[i].second.lowtx = ts[i].second.hightx = ts[i].second.lowtx - mintx.second;
      ts[i].second.lowty = ts[i].second.highty = ts[i].second.lowty - minty.second;
    }
  for( unsigned int i = 0; i < ts.size(); i++ )
    {
      if( ts[i].second.lowtx + ts[i].first->getheight() > maxtx.second )
        {
          maxtx.second = ts[i].second.lowtx + ts[i].first->getheight();
          maxtx.first = i;
        }
      if( ts[i].second.lowty + ts[i].first->getwidth() > maxty.second )
        {
          maxty.second = ts[i].second.lowty + ts[i].first->getwidth();
          maxty.first = i;
        }
    }
  return 0;
}

int alignimage::createimage( int aid, int slice, int channel,
                             std::vector< std::pair<tile*, transrange> > &ts,
                             std::vector< image2d* > &imcontribute )
{
  // create the image
  image2d im;
  im.create( imcontribute[0]->getheight(),
             imcontribute[0]->getwidth() );
  im.clear();
  std::vector< image2d* > imvec;
  imvec.clear();

  for( unsigned i = 0; i < ts.size(); i++ )
    {
      image2d *im;
      im = new image2d;
      ts[i].first->getimage( slice, channel, im );
      imvec.push_back( im );
    }

  for( unsigned int i = 0; i < ts.size(); i++ )
    addimage( imvec[i], imcontribute[i], ts[i].second, im );
  
  for( unsigned i = 0; i < imvec.size(); i++ )
    delete imvec[i];
    
  imvec.clear();
  // add to database
  db->add_aligned_image( aid, slice, channel, &im );
  im.del();
  db->get_aligned_image( aid, slice, channel, &im );
#ifndef _NOMAGICK_
  Magick::Image mim;
  im.getimage_normalized( mim );
  mim.equalize();
  char buff[256];
  sprintf(buff, "%i_%i_%i.png", aid, channel, slice );
  mim.write( buff );
#endif
  return 0;
}

int alignimage::calcblending( int x, int y, std::vector< image2d* > &cmaps,
                              image2d &immap, int h, int w,
                              std::vector< double > &blending )
{
  for( unsigned i = 0; i < blending.size(); i++ )
    blending[i] = 0;
  double counts;
  counts = 0;
  for( int i = x-1; i < x+2; i++ )
    {
      if( i > 0 && h > i )
        {
          for(int j = y-1; j < y+2; j++ )
            {
              if( j > 0 && w > j )
                {
                  double value;
                  value = immap.getpixel(i, j);
                  if( value > 0.5 )
                    {
                      for( unsigned k = 0; k < cmaps.size(); k++ )
                        blending[k] = blending[k] + cmaps[k]->getpixel(i,j);
                      counts += 1;
                    }
                }
            }
        }
    }
  if( counts == 0 )
    {
      std::cout << "ERROR blending 0 pixels\n";
      exit(0);
    }
  for( unsigned k = 0; k < cmaps.size(); k++ )
    blending[k] = blending[k] / counts;
  for( unsigned k = 0; k < cmaps.size(); k++ )
    if( blending[k] < 0.001 )
      blending[k] = 0.001;
  return 0;
}

int alignimage::addimage_mark( image2d &origmap,
                               image2d &immark, image2d &immap,
                               image2d &imoverlap, transrange p,
                               int imh, int imw)
{
  immark.clear();
  for( int i = 0; i < imh; i++ )
    {
      for( int j = 0; imw > j; j++ )
        {
          int ti, tj;
          if( origmap.getpixel( i, j) > 0.5 )
            {
              ti = (int) p.lowtx+i;
              tj = (int) p.lowty+j;
              immark.setpixel( ti, tj, 1 );
              if( immap.getpixel( ti, tj ) > 0.5 )
                imoverlap.setpixel( ti, tj, 1 );
              else
                immap.setpixel( ti, tj, 1 );
            }
        }
    }
  return 0;
}

int alignimage::addimage( image2d* in, image2d *contri, transrange p, image2d &im )
{
  for( int i = 0; i < in->getheight(); i++ )
    {
      for( int j = 0; j < in->getwidth(); j++ )
        {
          int ti, tj;
          ti = (int)p.lowtx+i;
          tj = (int)p.lowty+j;
          double value;
          value = in->getpixel(i,j) * contri->getpixel( ti, tj );
          value = value + im.getpixel( ti, tj );
          im.setpixel( ti, tj, value );
        }
    }

  return 0;
}

int alignimage::createalignset( int iteration,
                                std::vector<std::pair<tile*, transrange> > &ts )
{
  std::vector< int > row, col, transx, transy;
  alignedset *as;
  as = new alignedset;
  std::cout << "creating set in iteration: " << iteration << "\n";
  for( unsigned i = 0; i < ts.size(); i++ )
    {
      if( iteration == 0 )
        {
          row.push_back( ts[i].first->xcoor() );
          col.push_back( ts[i].first->ycoor() );
        }
      else
        {
          row.push_back( ts[i].first->getaid() );
        }
      transx.push_back( (int) ts[i].second.lowtx );
      transy.push_back( (int) ts[i].second.lowty );
      tileinfo i2;
      i2.transx = ts[i].second.lowtx;
      i2.transy = ts[i].second.lowty;
      if( iteration == 0 )
        {
          i2.x = ts[i].first->xcoor();
          i2.y = ts[i].first->ycoor();
        }
      else
        {
          i2.x = ts[i].first->getaid();
          i2.y = -1;
        }
      as->tiles.insert( i2 );
    }
  int aid;
  // create new database stuff
  if( iteration == 0 )
    aid = db->create_firstaln_set( row, col, transx, transy );
  else
    aid = db->create_align( iteration, row, transx, transy );
  as->aid = aid;
  as->connections.clear();
  while( (int)imsets.size() <= iteration )
    imsets.push_back(NULL);
  if( imsets[iteration] == NULL )
    imsets[iteration] = new imageset;
  imsets[ iteration ]->sets[aid] = as;
  return aid;
}

int alignimage::caltransforms( int start,
                               std::vector< tile* > &tiles,
                               std::map< std::pair<int, int>, transrange > &edges,
                               std::vector<std::pair<tile*, transrange> > &ts )
{
  int flag;
  flag = 1;
  ts.clear();
  int curnode = start;
  transrange tr;
  tr.lowtx = tr.lowty = tr.hightx = tr.highty = 0;
  std::pair< tile*, transrange > tp;
  tp.second = tr;
  tp.first = tiles[curnode];
  // push back transform and tile pair for a image set
  ts.push_back( tp );
  addtile( curnode, tiles, edges, ts, tr );
  
  return 0;
}

int alignimage::addtile( int curnode,
                         std::vector< tile* > &tiles,
                         std::map< std::pair<int, int>, transrange > &edges,
                         std::vector< std::pair<tile*, transrange> > &ts,
                         transrange basetrans)
{
  transrange myt;
  std::pair<tile*, transrange> tp;
  for(unsigned i = 0; i < tiles.size(); i++ )
    {
      if( (int)i != curnode )
        {
          std::pair<int, int> t;
          t.first = curnode;
          t.second = i;
          // if there is a connection
          if( edges.count( t ) > 0 )
            {
              int flag;
              flag = 0;
              // check if inserted already
              for( unsigned j = 0; j < ts.size(); j++ )
                if( ts[j].first == tiles[i] )
                  flag = 1;
              if( flag == 0 )
                {
                  myt.lowtx = basetrans.lowtx + edges[t].lowtx;
                  myt.lowty = basetrans.lowty + edges[t].lowty;
                  myt.hightx = basetrans.hightx + edges[t].hightx;
                  myt.highty = basetrans.highty + edges[t].highty;
                  tp.first = tiles[i];
                  tp.second = myt;
                  ts.push_back( tp );
                  addtile( i, tiles, edges, ts, myt );
                }
            }
        }
    }
  return 0;
}

int alignimage::colorsegments( std::vector<int> &flag, std::vector<int> &degree,
                               std::map< std::pair<int, int>, transrange > &edge )
{
  int color = 1;
  for( unsigned i = 0; i < degree.size(); i++ )
    {
      if( degree[i] == 0 )
        {
          if( flag[i] != 0 )
            std::cout << "merge error\n";
          graphdfs( flag, edge, i, color );
          color++;
        }
    }
  return color;
}

int alignimage::graphdfs( std::vector<int> &flag,
                          std::map< std::pair<int, int>, transrange > &edge,
                          int i, int color )
{
  if( flag[i] != 0 )
    return 0;
  flag[i] = color;
  for( std::map< std::pair<int, int>, transrange >::iterator it = edge.begin();
       it != edge.end(); it++ )
    {
      if( it->first.first == i )
        graphdfs(flag, edge, it->first.second, color );
    }
  return 0;
}

int alignimage::dotout()
{
  std::cout << "digraph G{\n";
  for( unsigned l = 0; l < imsets.size(); l++ )
    {
      if( imsets[l] != NULL )
        {
          std::cout << "subgraph cluster_level" << l << " {\n";
          std::cout << "label=\"iteration " << imsets[l]->level << "\";\n";
          for( std::map<int, alignedset* >::iterator it = imsets[l]->sets.begin();
               it != imsets[l]->sets.end(); it++ )
            {
              std::cout << "S" << it->first << " [shape=\"record\" label=\"{ {"
                        << it->first << " | size: " << it->second->sizex 
                        << "x" << it->second->sizey << "} | { ";
              for( std::set< tileinfo >::iterator st = it->second->tiles.begin();
                   st != it->second->tiles.end(); st++ )
                {
                  std::cout << st->x << " " << st->y << " | ";
                }
              std::cout << "} }\"];\n";
            }
          for( std::map<int, alignedset* >::iterator it = imsets[l]->sets.begin();
               it != imsets[l]->sets.end(); it++ )
            {
              for( unsigned i = 0; i < it->second->connections.size(); i++ )
                {
                  if( it->first < it->second->connections[i].aid )
                    {
                      std::cout << "TS" << it->first << "_" << it->second->connections[i].aid 
                                << " [shape=\"circle\" label=\"(" 
                                << it->second->connections[i].con.lowtx
                                << " " << it->second->connections[i].con.hightx << ")\\n("
                                << it->second->connections[i].con.lowty << " "
                                << it->second->connections[i].con.highty << ")\"];\n";
                      std::cout << "S" << it->first << " -> TS" << it->first << "_"
                                << it->second->connections[i].aid << ";\n";
                      std::cout << "TS" << it->first << "_" 
                                << it->second->connections[i].aid << " -> S"
                                << it->second->connections[i].aid << ";\n";
                      /*std::cout << "S" << it->first << " -> S" << it->second->connections[i].aid
                            << " [label=\"[" << it->second->connections[i].con.lowtx
                            << " " << it->second->connections[i].con.hightx << "]\\n["
                            << it->second->connections[i].con.lowty << " "
                            << it->second->connections[i].con.highty << "]\"];\n";*/
                    }
                }
            }
          std::cout << "}\n";
        }
    }
  std::cout << "}\n";
  return 0;
}


int alignimage::addimag_calcblend( std::vector< image2d* > cmaps, 
                                   image2d &immap, image2d &imoverlap )
{
  alignimage_blend_td *abt;
  int numthreads;
  int h = immap.getheight();
  int w = immap.getwidth();
  numthreads = 2;
  std::cout << "allocating memory to start " << numthreads << " threads\n";
  abt = new alignimage_blend_td[ numthreads ];
  std::vector< pthread_mutex_t * > lcs;
  pthread_mutex_t coutlock;
  pthread_mutex_init( &coutlock, NULL );
  pthread_barrier_t pbt;
  pthread_barrier_init( &pbt, NULL, numthreads );
  for( unsigned i = 0; i < cmaps.size(); i++ )
    {
      pthread_mutex_t *tlcs;
      tlcs = new pthread_mutex_t [numthreads-1];
      for( int j = 0; j < numthreads-1; j++ )
        pthread_mutex_init( tlcs+j, NULL );
      lcs.push_back( tlcs );
    }
  int amount;
  amount = immap.getheight()/numthreads;
  for( int i = 0; i < numthreads; i++ )
    {
      abt[i].lowlock.clear();
      abt[i].highlock.clear();
      if( i > 0 )
        for( unsigned j = 0; j < cmaps.size(); j++ )
          abt[i].lowlock.push_back( &(lcs[j][i-1]) );
      else
        for( unsigned j = 0; j < cmaps.size(); j++ )
          abt[i].lowlock.push_back( NULL );
      if( i != (numthreads-1) )
        for( unsigned j = 0; j < cmaps.size(); j++ )
          abt[i].highlock.push_back( &(lcs[j][i]) );
      else
        for( unsigned j = 0; j < cmaps.size(); j++ )
          abt[i].highlock.push_back( NULL );
      abt[i].sync = &pbt;
      abt[i].overlap = &imoverlap;
      abt[i].cmaps = &cmaps;
      abt[i].lowrow = i*amount;
      abt[i].highrow = (i+1)*amount;
      abt[i].clock = &coutlock;
      if( i == numthreads-1 )
        abt[i].highrow = immap.getheight();
    }
  std::cout << "data set, starting threads\n";
  void *ret;
  for( int i = 0; i < numthreads; i++ )
    pthread_create( &(abt[i].tid), NULL, blendthread_start, &(abt[i]) );
  for( int i = 0; i < numthreads; i++ )
    pthread_join( abt[i].tid, &ret );
  pthread_mutex_destroy( &coutlock);
  pthread_barrier_destroy( &pbt );
  for( unsigned j = 0; j < lcs.size(); j++ )
    {
      for( int i = 0; i < numthreads-1; i++ )
        pthread_mutex_destroy( lcs[j]+i );
      delete [] lcs[j];
    }
  lcs.clear();
  delete [] abt;

  // normalize the pixels
  std::cout << "normalizing...\n";
  std::vector< double > blending;
  blending.resize( cmaps.size(), 0 );
  for( int i = 0; i < h; i++ )
    {
      for( int j = 0; j < w; j++ )
        {
          if( immap.getpixel(i,j) > 0.5 )
            {
              for( unsigned k = 0; k < cmaps.size(); k++ )
                blending[k] = cmaps[k]->getpixel( i, j );
              double sum;
              sum = 0;    
              for( unsigned k = 0; k < cmaps.size(); k++ )
                sum = sum + blending[k];
              if( sum == 0 )
                {
                  std::cout << "sum = 0\n";
                  exit(0);
                }
              for( unsigned k = 0; k < cmaps.size(); k++ )
                blending[k] = blending[k]/sum;
              for( unsigned k = 0; k < cmaps.size(); k++ )
                cmaps[k]->setpixel( i, j, blending[k] );
            }
        }
    }

  std::cout << "normalization done\n";
  // assert that every point either sums to zero or one
  for( int i = 0; i < h; i++ )
    {
      for( int j = 0; j < w; j++ )
        {
          double sum;
          sum = 0;
          for( unsigned k = 0; k < cmaps.size(); k++ )
            sum = sum + cmaps[k]->getpixel( i, j );
          if( immap.getpixel(i,j) > 0.5 )
            {
              if( fabs(sum - 1 ) > 0.000001 )
                {
                  std::cout << "ERROR2 contribution calc\n";
                  std::cout << "sum = " << sum << "\n";
#ifndef _NOMAGICK_
                  Magick::Image mim;
                  immap.getimage(mim);
                  mim.write("immap.png");
#endif
                  std::cout << "i=" << i << " j=" << j << "\n";
                  exit(0);
                }
            }
          else if( sum > 0.0001 )
            {
              std::cout << "ERROR1 contributeion calc\n";
              exit(0);
            }
        }
    }
  std::cout << "contribution check passed\n";
  return 0;
}
