#ifndef _NOMAGICK_
#include <Magick++.h>
#endif

#define GAUSSSIGMA 0.002

#include <iostream>
#include <string>
#include <fstream>
#include <math.h>
#include "imagealign.h"
#include "image2d.h"
#include <fftw3.h>
#include <tclap/CmdLine.h>
#include "bdbenv.h"
#include "dbimage.h"
#include "pso.h"
#include "histequalize.h"
#include "tile.h"
#include "tileimage.h"
#include "alignimage.h"
#include "dbtmpfeature.h"
#include "som.h"
#include "texturefeature.h"
#include <algorithm>

int calcfeatures( double low, double high,
                  image2d *mask, int aid,
                  std::set<int> *sl, int mych,
                  dbalnimage *adb, som &mysom);


int calchiststretch( double &avg, double &var, image2d *mymask, int aid,
                     std::set<int> *sl, int mych, dbalnimage *adb, bdbenv *bdb);

int calchiststretch( double &avg, double &var, image2d *mymask, int aid,
                     std::set<int> *sl, int mych, dbalnimage *adb );

int calchiststretch( double &avg, double &var, image2d *mymask, int aid,
                     std::set<int> *sl, int mych, dbalnimage *adb )
{
  image2d myim;
  avg = var = 0;
  int height, width;
  long int count;
  count = 0;
  std::vector< int > todo;
  height = mymask->getheight();
  width = mymask->getwidth();
  for(int i = 0; i < height*width; i++ )
    if( mymask->getimageptr()[i] > 0.5 )
      todo.push_back( i );
  for( std::set<int>::iterator it = sl->begin(); it != sl->end(); it++ )
    {
      count++;
      adb->get_aligned_image( aid, *it, mych, &myim );
      myim.gaussian( GAUSSSIGMA );
      for(unsigned int i = 0; i < todo.size(); i++ )
        avg = avg + myim.getimageptr()[i];
      std::cout << count << ": " << avg / ( (double)(todo.size()*count) ) << "\n";
    }
  avg = avg / ( (double)(todo.size()*sl->size()) );
  for( std::set<int>::iterator it = sl->begin(); it != sl->end(); it++ )
    {
      adb->get_aligned_image( aid, *it, mych, &myim );
      myim.gaussian( GAUSSSIGMA );
      for(unsigned int i = 0; i < todo.size(); i++ )
        var = var + pow( avg - myim.getimageptr()[i], 2 );
    }
  var = var / ( (double)(todo.size()*sl->size()) );
  return 0;
}

int calcfeatures( double avg, double var,
                  image2d *mask, int aid,
                  std::set<int> *sl, int mych,
                  dbalnimage *adb, som &mysom)
{
  std::vector< int > s;
  for( std::set<int>::iterator it = sl->begin(); it != sl->end(); it++ )
    {
      if( *it < 16 )
        s.push_back( *it );
    }
  std::sort( s.begin(), s.end() );
  image2d *ims[3];
  double feature[32+6];
  double scale;
  scale = 1.0/( var * 0.15 );
  std::cout << "starting calc features\n";
  for(int i = 0; i < 3; i++ )
    {
      ims[i] = new image2d;
      adb->get_aligned_image( aid, s[i], mych, ims[i] );
      ims[i]->gaussian( GAUSSSIGMA );
      ims[i]->add( -avg );
      ims[i]->multiply( scale );
    }
  std::cout << "first load complete\n";
  image2d ir, ig, ib;
  int width, height;
  texturefeature tf;
  width = ims[0]->getwidth();
  height =  ims[0]->getheight();
  ir.create( ims[0]->getheight(), ims[0]->getwidth() );
  ig.create( ims[0]->getheight(), ims[0]->getwidth() );
  ib.create( ims[0]->getheight(), ims[0]->getwidth() );
  Magick::Geometry geo;
  geo.width( ims[0]->getwidth() );
  geo.height( ims[0]->getheight() );
  Magick::Image mim( geo, Magick::Color("black") );
  mim.type( Magick::TrueColorType );
  
  for( unsigned int i = 1; i < s.size()-1; i++ )
    {
      std::cout << "initalizing features\n" << std::flush;
      tf.init( ims, mask );
      Magick::PixelPacket *pixel_cache 
        = mim.getPixels(0,0,width,height);
      Magick::ColorRGB rgb(0,0,0);
      int count, mcount;
      count = mcount = 0;
      std::cout << "calculating fatures\n" << std::flush;
      for( int x = 0; x < ims[0]->getheight(); x = x + 1 )
        {
          for( int y = 0; y < ims[0]->getwidth(); y = y + 1 )
            {
              int ret;
              ret = tf.calc( x, y, feature);
              if( ret == 0 )
                {
                  count++;
                  std::vector< int > myrgb;
                  int cl;
                  cl = mysom.classify( feature, myrgb );
                  if( cl < 0 )
                    {
                      rgb.red(0);
                      rgb.blue(0);
                      rgb.green(0);
                    }
                  else
                    {
                      rgb.red( ((double)myrgb[0]+1.0)/6.0 );
                      rgb.green( ((double)myrgb[1]+1.0)/6.0 );
                      rgb.blue( ((double)myrgb[2]+1.0)/6.0 );
                    }
                  pixel_cache[x*width+y] = rgb;
                }
              else if( ret == -2 )
                mcount++;
            }
        }
      std::cout << "contain " << count << "/" << mcount << " good pixels\n";
      mim.syncPixels();
      char buffer[512];
      sprintf( buffer, "class%i.tif", s[i] );
      mim.write( buffer );
      image2d *zz;
      zz = ims[0];
      ims[0] = ims[1];
      ims[1] = ims[2];
      ims[2] = zz;
      int next;
      next = i + 2;
      if( next < (int)s.size() )
        {
          adb->get_aligned_image( aid, s[next], mych, ims[2] );
          ims[2]->gaussian( GAUSSSIGMA );
          ims[2]->add( -avg );
          ims[2]->multiply( scale );
        }
    }
  delete ims[0];
  delete ims[1];
  delete ims[2];

  return 0;
}

int main( int argc, char *argv[] )
{
  try
    {
      ///set up the command line parsing
      TCLAP::CmdLine cmd("Image alignment program", ' ', "0.001");

      TCLAP::ValueArg<std::string> dbdir("d", "db", "database dirctory", true, "", "databaase");

      TCLAP::ValueArg<int> aid("a", "aid", "aligned image id", true, 0, "aligned image id");
      TCLAP::ValueArg<int> channel("c", "channel", "channel to get", true, 0, "channel" );
      TCLAP::ValueArg<std::string> somxml("s", "som", "som xml file", true, "", "som");
      TCLAP::ValueArg<std::string> maskim("m", "mask", "image mask", true, "", "imagemask");

#ifndef _NOMAGICK_
      Magick::Image mim;
#endif

      cmd.add(dbdir);
      cmd.add(aid);
      cmd.add(channel);
      cmd.add(somxml);
      cmd.add(maskim);

      cmd.parse(argc, argv);

      bdbenv bdb;
      dbimage imstore;
      dbalnimage dba;
      som mysom;

      mysom.init( NULL, somxml.getValue() );

      bdb.setopenflags( DB_CREATE | DB_THREAD |
                        DB_INIT_MPOOL | DB_PRIVATE );
      bdb.open( dbdir.getValue() );
      imstore.init( bdb.envptr() );
      dba.init( bdb.envptr() );

      bdb.adddb( &imstore );
      bdb.adddb( &dba );
      bdb.cleanlog();
      
      int myaid;
      std::cout << "checking aid\n";
      myaid = aid.getValue();
      //check for valid id
      if( dba.checkaid( myaid ) == 1 )
        {
          std::set<int> sl, ch;
          std::cout << "getting slices\n";
          dba.get_slices(myaid, sl );
          std::cout << "getting channels\n";
          dba.get_channels(myaid, ch );
          std::cout << "checking for valid channel\n";
          int mych;
          mych = channel.getValue();
          image2d mymask;
          if( ch.count( mych ) > 0 )
            {
              Magick::Image mim;
              mim.read( maskim.getValue().c_str() );
              mymask.load( mim, CH_GRAY );
              /*
              double low, high;
              //std::cout << "getting mask\n";
              //dba.get_aligned_mask( myaid, &mymask );
              std::cout << "calculating histogram\n";
              calchiststretch( low, high, &mymask, aid.getValue(),
                               &sl, mych, &dba);
              */
              calcfeatures( 0.00300228, 0.00384478, &mymask, aid.getValue(),
                            &sl, mych, &dba, mysom );
            }
          else
            {
              std::cout << "invalid channel\n";
            }
        }
      else
        {
          std::cout << "invalid aid\n";
        }


      dba.close();
      imstore.close();
      bdb.close();
    }
  catch (TCLAP::ArgException &e)
    {
      std::cout << "ERROR: " << e.error() 
                << " for arg " << e.argId() << "\n";
      return 0;
    }
  return 0;
}
