#include "mutualinformation.h"

mutualinformation::mutualinformation()
{
  ahist = new double[MIPMFSZ];
  bhist = new double[MIPMFSZ];
  jhist = new double[MIPMFSZ*MIPMFSZ];
}

mutualinformation::~mutualinformation()
{
  delete [] ahist;
  delete [] bhist;
  delete [] jhist;
}

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

  score = 0;
  double Ha, Hb;
  Ha = Hb = 0;

  subimage imAcut, imBcut;

  // 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;

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


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

  calcpmf( ahist, imAcut );
  calcpmf( bhist, imBcut );
  calcjpmf( jhist, imAcut, imBcut );
  
  double ha, hb, hab;
  ha = hb = hab = 0;
  for( int i = 0; i < MIPMFSZ; 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 < MIPMFSZ*MIPMFSZ; i++ )
    if( jhist[i] > 0 )
      hab = hab + jhist[i]*log(jhist[i]);

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

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

  score = 0;
  double Ha, Hb;
  Ha = Hb = 0;
  subimage imAcut, imBcut;
  subimage mapAcut, mapBcut;
  image2d map;

  // 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;

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

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


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

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

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

  calcpmf( ahist, imAcut, map );
  calcpmf( bhist, imBcut, map );
  calcjpmf( jhist, imAcut, imBcut, map );
  
  double ha, hb, hab;
  ha = hb = hab = 0;
  for( int i = 0; i < MIPMFSZ; 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 < MIPMFSZ*MIPMFSZ; i++ )
    if( jhist[i] > 0 )
      hab = hab + jhist[i]*log(jhist[i]);

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

int mutualinformation::calcpmf( double *hist, subimage &im )
{
  for( int i = 0; i < MIPMFSZ; i++ )
    hist[i] = 0;
  double ipart, fpart, sum;
  sum = 0;
  int w, h;
  h = im.getsizex();
  w = im.getsizey();
  for( int i = 0; i < h; i = i + MISTEP )
    {
      for( int j = 0; j < w; j = j + MISTEP )
        {
          double v;
          v = im.getpixel(i,j);
          fpart = modf(v, &ipart );
          hist[(int)ipart] += (1-fpart);
          if( ipart + 1 < MIPMFSZ )
            hist[(int)ipart+1] += fpart;            
        }
    }
  for( int i = 0; i < MIPMFSZ; i++ )
    sum = sum + hist[i]; 
  for( int i = 0; i < MIPMFSZ; i++ )
    hist[i] = hist[i] / sum;
  return 0;
}

int mutualinformation::calcpmf( double *hist, subimage &im, image2d &map )
{
  for( int i = 0; i < MIPMFSZ; i++ )
    hist[i] = 0;
  double ipart, fpart, sum;
  sum = 0;
  int w, h;
  h = im.getsizex();
  w = im.getsizey();
  for( int i = 0; i < h; i = i + MISTEP )
    {
      for( int j = 0; j < w; j = j + MISTEP )
        {
          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 < MIPMFSZ )
                hist[(int)ipart+1] += fpart;
            }
        }
    }
  for( int i = 0; i < MIPMFSZ; i++ )
    sum = sum + hist[i]; 
  for( int i = 0; i < MIPMFSZ; i++ )
    hist[i] = hist[i] / sum;
  return 0;
}

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

int mutualinformation::calcjpmf(double *hist, subimage &A, subimage &B, image2d &map)
{
  double sum, ipart1, ipart2, fpart1, fpart2;
  sum = 0;
  int h, w;
  h = A.getsizex();
  w = A.getsizey();
  if( h != B.getsizex() || w != B.getsizey() )
    {
      std::cout << "ERROR size in joint pmf calculation failed\n";
      exit(0);
    }
  for( int i = 0; i < MIPMFSZ*MIPMFSZ; i++ )
    hist[i] = 0;
  for( int i = 0; i < h; i = i + MISTEP )
    {
      for( int j = 0; j < w; j = j + MISTEP )
        {
          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*MIPMFSZ+(int)ipart2] += (1.0-fpart1)*(1.0-fpart2);
              if( ipart2 + 1 < MIPMFSZ )
                hist[(int)ipart1*MIPMFSZ+(int)ipart2+1] += (1.0-fpart1)*(fpart2);
              if( ipart1 + 1 < MIPMFSZ )
                hist[(int)(ipart1+1)*MIPMFSZ+(int)ipart2] += (fpart1)*(1.0-fpart2);
              if( ipart2 + 1 < MIPMFSZ && ipart1 + 1 < MIPMFSZ )
                hist[(int)(ipart1+1)*MIPMFSZ+(int)ipart2+1] += (fpart1)*(fpart2);
            }
        }
    }
  for( int i = 0; i < MIPMFSZ*MIPMFSZ; i++ )
    sum = sum + hist[i];
  for( int i = 0; i < MIPMFSZ*MIPMFSZ; i++ )
    hist[i] = hist[i] / sum;
  return 0;
}
