/************************************************************************

   File:		DivergenceMap.cpp

   Author(s):		Pavel Dimitrov

   Created:		26 Jun 2002

   Last Revision:	$Date: 2002/07/23 21:02:53 $

   Description:	

   $Revision: 1.3 $

   $Log: DivergenceMap.cpp,v $
   Revision 1.3  2002/07/23 21:02:53  pdimit
   The branch segmentation of the thinned DivMap is improved --
   the two pixel branches are now included. Still, the degenerate case
   of the square of junction-points is not taken care of.

   Also, a DiscreteDivergenceSkeleton is created which still does not know
   of left and right for the branches, i.e. the contour has not been cut.

   Revision 1.2  2002/06/30 01:22:31  pdimit
   Added DivergenceSkeletonMaker and a test for it. Incomplete functionality
   and the testing is mostly in DivergenceSkeletonMaker.cpp. Only the
   thinning of the DivergenceMap is implemented.

   Revision 1.1  2002/06/26 11:51:37  pdimit
   Implemented the DivergenceMap class. It is supposed to be a base
   class for other implementations. It has very dumb algorithms, but
   which seem to work just fine. Look at testDivergenceMap to see how
   to use it. Also, testSimpleShapeMaker has a much nicer interface --
   exactly the same as for testDivergenceMap...


	
   Copyright (c) 2002 by Pavel Dimitrov, Centre for Intelligent Machines,
   McGill University, Montreal, QC.  Please see the copyright notice
   included in this distribution for full details.

 ***********************************************************************/

#define DIVERGENCEMAP_C

#include "DivergenceMap.h"
#include "Vector.h"
#include <math.h>

using namespace sg;

// eps is the half side of the square
double DivergenceMap::getValue(double x, double y, double eps) const{
  double ret=0;

//    double sqrt_2 = sqrt(2.0)/2.0;
  // The 8 normals of the square
  // 1 2 3 
  // 8   4
  // 7 6 5
//    Vector normals[] = {Vector(-sqrt_2, sqrt_2), Vector(0,1),  // 1 2
//  		      Vector(sqrt_2, sqrt_2), Vector(1,0),   // 3 4
//  		      Vector(sqrt_2, -sqrt_2), Vector(0,-1), // 5 6
//  		      Vector(-sqrt_2, -sqrt_2), Vector(-1,0) // 7 8
//    };


//    double deps = sqrt_2 * eps; // on the diagonal entries

//    double gradeps = 0.001*eps;

//    ret += normals[0].dot(dt->getGradValue(x-deps, y+deps, gradeps)); // 1
//    ret += normals[1].dot(dt->getGradValue(x, y+eps, gradeps));     // 2
//    ret += normals[2].dot(dt->getGradValue(x+deps, y+deps, gradeps)); // 3
//    ret += normals[3].dot(dt->getGradValue(x+eps, y, gradeps));     // 4
//    ret += normals[4].dot(dt->getGradValue(x+deps, y-deps, gradeps)); // 5
//    ret += normals[5].dot(dt->getGradValue(x, y-eps, gradeps));     // 6
//    ret += normals[6].dot(dt->getGradValue(x-deps, y-deps, gradeps)); // 7
//    ret += normals[7].dot(dt->getGradValue(x-eps, y, gradeps));     // 8

  const int n = 10;
  const double Pi_n = 2.0 * 3.141592654 / (double) n;
  for(int i = 0; i < n; i++){
    
    double xx,yy;
    
    xx = cos(((double)i)*Pi_n);
    yy = sin(((double)i)*Pi_n);

    Vector N;
    N.x = xx; N.y = yy; //(xx, yy);
    
    xx = x + eps*xx;
    yy = y + eps*yy;

    Vector v = dt->getGradValue(xx, yy, 0.001*eps);

    ret += N.dot(v);
  }


  return Pi_n * ret;
}
