#include "tnt.h" /* This bit is to use the TNT package for arrays and matrices */
#define LOG_0 = -100;
#define LOGP(x,y) (((x)>(y))?(x)+log1p(exp((y)-(x))):(y)+log1p(exp((x)-(y))))
using namespace std;
/**********************************/
/* A useful elementary function ***/
/**********************************/
int fathernode(int n,int division,int grid_side) {
/* Computes the index of the father of a node in a series of nested histograms */
return( ((n/grid_side)/division)*(grid_side/division)+((n % grid_side)/division));
}
/***********************************
the following inputs are common to all functions, and are parameters:
double multires_g: the branching process prior, which controls the probability at each node of the tree to give birth to an offspring. We assume, as is the case in the original paper, that this prior is uniform over the tree.
int multires_d: the number of offsprings at each node. As is the case with the original paper, we assume that this number is constant (typically 4 if images are successively divided into 4 patches)
int multires_D: maximum depth of considered trees.
double temperature_t: width parameter for all kernels
double factor: factor is a not-so-useful parameter which can normalize counts to a certain value before using the Hilbertian metrics JD, H1, etc.. it overlaps a bit with temperature_t.
the two main inputs are :
TNT::Array2D< int > & data1,
TNT::Array2D< int > & data2,
are both arrays of integers which have d^D lines and a fixed number of columns which is equal to the total number of considered components ** +1 **. In other words, each line of data1 and data2 stand for the integral histogram of counts of each considered patch, plus the last column which counts their total contribution. As an example, if we divide images in 4 patchs, over a depth of 3, and only consider 512 colors, data1 and data2 will have
- 4^3 = 64 lines
- 512 + 1 = 513 columns, where the first 512 columns stand for each color's individual count in this patch, plus the last column which the is the sum of the first 512 coefficients of that line.
we have used the TNT::Array2D implementation for matrices but any other implementation of integral matrices should do.
**********************************/
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// JENSEN DIVERGENCE
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double mJD(TNT::Array2D< int > & data1, TNT::Array2D< int > & data2, double multires_g, int multires_d, int multires_D, double temperature_t, double factor){
const int p=data1.dim1(),N=data1.dim2()-1,D=(int) pow((double)multires_d,(double)multires_D);
int i,j,k,m;
double fact1=0,fact2=0;
for (i=0;i0) {fact1=1/factor; fact2=1/factor;}
TNT::Array2D< double > data1d(p,N);
TNT::Array2D< double > data2d(p,N);
TNT::Array2D< double > data1d2(p,N);
TNT::Array2D< double > data2d2(p,N);
for (i=0;i

lnP(D*D,0.0);
vector lnP2(D*D,0.0);
vector lnP2p(D*D,0.0);
double sum,sum1,sum2;
for (j=0;j<=multires_D;j++) {
const int iterations=(int) pow((double)multires_d,(multires_D-j));
const int iterations2=iterations*iterations;
for (i=0;i0 ? .5*(data1d[i][k]+data2d[i][k])*log(.5*(data1d[i][k]+data2d[i][k])) : 0);
sum1+= ( data1d[i][k]>0 ? data1d[i][k]*log(data1d[i][k]) : 0);
sum2+= ( data2d[i][k]>0 ? data2d[i][k]*log(data2d[i][k]) : 0);
}
lnP.at(i)=(sum-.5*(sum1+sum2))/temperature_t;
}
if (j>0) {
for (i=0;i & data1, TNT::Array2D< int > & data2, double multires_g, int multires_d, int multires_D, double temperature_t, double factor){
const int p=data1.dim1(),N=data1.dim2()-1,D=(int) pow((double)multires_d,(double)multires_D);
int i,j,k,m;
// The last column is the counter.
double fact1=0,fact2=0;
for (i=0;i0) {fact1=1/factor; fact2=1/factor;}
TNT::Array2D< double > data1d(p,N);
TNT::Array2D< double > data2d(p,N);
TNT::Array2D< double > data1d2(p,N);
TNT::Array2D< double > data2d2(p,N);
for (i=0;i

lnP(D*D,0.0);
vector lnP2(D*D,0.0);
vector lnP2p(D*D,0.0);
double sum;
for (j=0;j<=multires_D;j++) {
const int iterations=(int) pow((double)multires_d,(multires_D-j));
const int iterations2=iterations*iterations;
for (i=0;i0 ? pow(data1d[i][k]-data2d[i][k],2)/(data1d[i][k]+data2d[i][k]) : 0);
}
lnP.at(i)=-sum/temperature_t;
}
if (multires_g<0.00001) {for (i=0;i0) {
for (i=0;i & data1, TNT::Array2D< int > & data2, double multires_g, int multires_d, int multires_D, double a, double b, double rho, double c, double factor){
const int p=data1.dim1(),N=data1.dim2()-1,D=(int) pow((double)multires_d,(double)multires_D);
int i,j,k,m;
double fact1=0,fact2=0;
for (i=0;i0) {fact1=1/factor; fact2=1/factor;}
TNT::Array2D< double > data1d(p,N);
TNT::Array2D< double > data2d(p,N);
TNT::Array2D< double > data1d2(p,N);
TNT::Array2D< double > data2d2(p,N);
for (i=0;i

lnP(D*D,0.0);
vector lnP2(D*D,0.0);
vector lnP2p(D*D,0.0);
double sum;
for (j=0;j<=multires_D;j++) {
const int iterations=(int) pow((double)multires_d,(multires_D-j));
const int iterations2=iterations*iterations;
for (i=0;i0) {
for (i=0;i