// walsh-nd.c           Section: 8.5
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
// Image Algebra C++ implementation of the N-D fast
//      Walsh transform using spatial transforms
//

#include "walsh.h"
#include "Closure.h"
#include "ClosureI.h"


//==============================================================
// Utility routines here
//==============================================================


////////////////////////////////////////////////////////////////////
// idiv2 - integer divide x by 2^n

inline int  idiv2( int x, int n ) { return x >> n; }


////////////////////////////////////////////////////////////////////
// ipow2 - integer calculate  2^n

inline int  ipow2( int n ) { return 1 << n; }


////////////////////////////////////////////////////////////////////
// ilog2 - calculate the log base 2 of an integer power of 2

int  ilog2(int i)
{
    int j=-1;
    while( i ) {
	j++;
	i = i >> 1;
    }
    return j;
}


////////////////////////////////////////////////////////////////////
// pointsetDimenWidth -
//    Calculate the  WIDTH  for a specific dimension, of the
//    bounding box for a given pointset.
//    NOTE that there is no guarantee that the  inf()  or  sup()
//    of a pointset are actually members of the pointset.

inline int
pointsetDimenWidth(const IA_Set<IA_Point<int> > &ips, int dimen) {
    return (ips.sup()-ips.inf()+1)[dimen];
}


//==============================================================
//
// walsh specific utility routines follow
//
//==============================================================


////////////////////////////////////////////////////////////////////
// rho  - permutation function used in  reordering of input that the
//      Fast Fourier Transform requires  

int
rho(int i, int n)
{
    // Reverse the bit representation of i in field log2(n) bits
    int j=0, m = i;
    for( int k=ilog2(n); k!=0; k-- ) {
	j = (j<<1) + (m&1);
	m = m >> 1;
    }
    return j;
}


////////////////////////////////////////////////////////////////////
// class R_m  - row permutation function, implemented as an
//     iac++ closure.  Instances provide a point to point mapping
//     object that looks like a function due to the definition of
//     operator().
//     m  should be a power of 2.
//     d  is the dimension utilized for mapping purposes.

class R_m: public IA_Closure<IA_Point<int>,IA_Point<int> > {
private:
    int  m;   // dimension width
    int  d;   // Which dimension are we dealing with ?
public:
    R_m( int dimenWidth, int whichDimen )
    : m(dimenWidth), d(whichDimen) {}

    IA_Point<int>  operator() (const IA_Point<int> &ip) const {
	IA_Point<int> rmip = ip;
	rmip[d] = rho(rmip[d], m);
	return  rmip;
    }

    IA_Closure<IA_Point<int>,IA_Point<int> > *clone_self() const {
	return new R_m(*this);
    }
};


////////////////////////////////////////////////////////////////////
// class F_i -
//     An iac++ closure used to implement the f() function as
//     defined for the walsh.
//     Instances of this class can be  COMPOSED  with an image,
//     in a similar fashion as a function.
//     Note that the quantity 2^(i-1) is computed once, when
//     an instance of this class is created.

class F_i: public IA_Closure<IA_Point<int>,IA_Point<int> > {
private:
    int  i_1;               // i-1
    int  ipow2_i_1;         // pre-compute 2^(i-1)
    int  d;                 // What dimension are we dealing with ?
public:
    F_i( int i, int whichDimen )
    : i_1(i-1), ipow2_i_1(ipow2(i-1)), d(whichDimen) {}

    IA_Point<int>  operator() (const IA_Point<int> &ip) const {
	IA_Point<int> fiip = ip;
	int val = idiv2(fiip[d], i_1) %2 ? 0 : ipow2_i_1;
	fiip[d] += val;
	return fiip;
    }

    IA_Closure<IA_Point<int>,IA_Point<int> > *clone_self() const {
        return new F_i(*this);
    }
};


////////////////////////////////////////////////////////////////////
// class G_i -
//     Instances represent the g() function as defined for the walsh.
//     See f() above for other notes.

class G_i: public IA_Closure<IA_Point<int>,IA_Point<int> > {
private:
    int  i_1;               // i-1
    int  ipow2_i_1;         // pre-compute 2^(i-1)
    int  d;                 // What dimension are we dealing with ?
public:
    G_i( int i, int whichDimen )
    : i_1(i-1), ipow2_i_1(ipow2(i-1)), d(whichDimen) {}

    IA_Point<int>  operator() (const IA_Point<int> &ip) const {
	IA_Point<int> giip = ip;
	int val = idiv2(giip[d], i_1) %2 ? ipow2_i_1 : 0;
	giip[d] -= val;
	return giip;
    }

    IA_Closure<IA_Point<int>,IA_Point<int> > *clone_self() const {
        return new G_i(*this);
    }
};


////////////////////////////////////////////////////////////////////
// class w_image -
//     An image closure implementation of the omega function
//     associated with the walsh.

class  w_image: public IA_ClosureI<IA_Point<int>,float> {
private:
    int     w_i;    // The index this instance is concerned with.
    int     d;      // Which dimension are we working on ?
public:
    w_image( IA_Set<IA_Point<int> > ps, int i, int whichDimen )
    :IA_ClosureI<IA_Point<int>,float>(ps), w_i(i-1), d(whichDimen) {};

    float operator() (const IA_Point<int>& xy) const {
	int evenodd = (idiv2(xy[d],w_i) % 2); // 0 if even

	return (evenodd ? -1 : 1);
    }
    
    IA_ClosureI<IA_Point<int>,float> *clone_self() const {
        return new w_image(*this);
    }
};


////////////////////////////////////////////////////////////////////
// walsh - Calcuate the Fast Fourier Transform.
//   This version is very similar to the one presented in the
//   IA Handbook, except that has been genealized to N dimensions.
//   Because of this the transpose() is no longer needed.

IA_Image<IA_Point<int>,float>
walsh(const IA_Image<IA_Point<int>,float> & image) {

    IA_Image<IA_Point<int>,float>  a(image);     // Make a local copy

    // Calculate the width of each dimension.
    // It should be a power of 2.
    IA_Point<int> a_dim_width( a.domain().sup() - a.domain().inf() + 1 );


    // We will do the forward transform for each dimension of the image.
    // iac++ library first dimension index is zero (0), not one (1).

    IA_Set<IA_Point<int> >  aps = a.domain();

    for ( int d=0; d< a_dim_width.dim(); d++) {
	R_m  rm(a_dim_width[d], d);
	a = compose(a,rm,aps);  // permutate rows of the image

	for ( int i=1;
	      i<=ilog2( a_dim_width[d] );
	      i++ ) {

	    w_image w(aps,i,d);  // Create an omega image
	    F_i fi(i,d);         // Create a  f()  `function`
	    G_i gi(i,d);         // Create a  g()  `function`

	    a = compose(a,gi,aps) + w*compose(a,fi,aps);
	}
    }
    return a;
}


////////////////////////////////////////////////////////////////////
// walsh_inv - Calculate the inverse walsh image using the forward walsh

IA_Image<IA_Point<int>,float>
walsh_inv(const IA_Image<IA_Point<int>,float> & a) {
    IA_Point<int>
	a_dim_width( a.domain().sup() - a.domain().inf() + 1 );

    double nm = product( a_dim_width );

    return walsh( a ) / nm;
}
