head	1.3;
access;
symbols;
locks; strict;
comment	@// @;


1.3
date	95.01.06.06.18.28;	author rjj;	state Exp;
branches;
next	1.2;

1.2
date	94.06.05.02.40.29;	author rjj;	state Exp;
branches;
next	1.1;

1.1
date	94.02.20.04.22.53;	author rjj;	state Exp;
branches;
next	;


desc
@Implemention of the walsh transform, based on the fft-nd.c code.
This version utilizes composition.  A faster version could
be produced based on the fft-nd-iter.c code.
@


1.3
log
@Added standard header.
Changed names of functions to mirror fft-nd.c
Inlined usage of F(), G() and W() objects.
@
text
@// walsh-nd.c           Section: 8.5, Walsh Transform
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
// Image Algebra C++ implementation of the N-dimensional
//      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
//      Walsh 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  - 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: public IA_Closure<IA_Point<int>,IA_Point<int> > {
private:
    int  m;   // dimension width
    int  d;   // Which dimension are we dealing with ?
public:
    R( 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(*this);
    }
};


////////////////////////////////////////////////////////////////////
// class F -
//     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: 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( 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(*this);
    }
};


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

class G: 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( 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(*this);
    }
};


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

class  W: 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( int i, int whichDimen, IA_Set<IA_Point<int> > ps )
    :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(*this);
    }
};


////////////////////////////////////////////////////////////////////
// walsh - Calcuate the Walsh 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++) {
	// permutate rows of the image
	a = compose( a, R(a_dim_width[d], d), aps );

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

	    a = compose( a,  G(i,d), aps ) +
		W(i,d,aps) * compose( a, F(i,d), 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;
}
@


1.2
log
@Type names harmonized to v1.0 standards.
Standard header adopted.
@
text
@d1 1
a1 1
// walsh-nd.c           Section: 8.5
d4 2
a5 2
// Image Algebra C++ implementation of the N-D fast
//      Walsh transform using spatial transforms
a15 2


d64 1
a64 1
//      Fast Fourier Transform requires  
d66 1
a66 2
int
rho(int i, int n)
d79 1
a79 1
// class R_m  - row permutation function, implemented as an
d86 1
a86 1
class R_m: public IA_Closure<IA_Point<int>,IA_Point<int> > {
d91 1
a91 1
    R_m( int dimenWidth, int whichDimen )
d101 1
a101 1
	return new R_m(*this);
d107 1
a107 1
// class F_i -
d115 1
a115 1
class F_i: public IA_Closure<IA_Point<int>,IA_Point<int> > {
d121 1
a121 1
    F_i( int i, int whichDimen )
d132 1
a132 1
        return new F_i(*this);
d138 1
a138 1
// class G_i -
d142 1
a142 1
class G_i: public IA_Closure<IA_Point<int>,IA_Point<int> > {
d148 1
a148 1
    G_i( int i, int whichDimen )
d159 1
a159 1
        return new G_i(*this);
d165 1
a165 1
// class w_image -
d169 1
a169 1
class  w_image: public IA_ClosureI<IA_Point<int>,float> {
d174 1
a174 1
    w_image( IA_Set<IA_Point<int> > ps, int i, int whichDimen )
d184 1
a184 1
        return new w_image(*this);
d190 1
a190 1
// walsh - Calcuate the Fast Fourier Transform.
d211 2
a212 2
	R_m  rm(a_dim_width[d], d);
	a = compose(a,rm,aps);  // permutate rows of the image
d218 2
a219 5
	    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);
@


1.1
log
@Initial revision
@
text
@d1 5
a5 3
// Emacs: -*- C++ -*-
// walsh-nd.c   Image Algebra C++ implementation of the N-D fast
//            Walsh transform using spatial transforms
a6 1
// Cross-references: IA with Apps. pp. x-X; IA Handbook pp. y-Y
d8 3
a10 1
static char rcsid[] = "$Id$";
a11 1
#include "imagealgebra.h"
a12 1

d52 1
a52 1
pointsetDimenWidth(const IA_IntPointSet &ips, int dimen) {
a56 12
////////////////////////////////////////////////////////////////////
// coordProduct - Calculate the product of a point\'s coordinates.
double
coordProduct( const IA_IntPoint &ip ) {
    double ipprod = 1;
    for ( int d=0; d < ip.dim(); d++ ) {
	ipprod *= ip[d];
    }
    return ipprod;
}


a77 11

//==============================================================
// Sample code as suggested in the IA Handbook
//    int h, j=0, m=i;
//    for( int l=0; l<=ilog2(n)-1 ; l++ ) {
//	h = int( floor((double(m))/2.0) );
//	j = 2*j + m - 2*h;
//	m = h;
//    }
//    return j;
//==============================================================
d89 1
a89 1
class R_m: public IA_Closure<IA_IntPoint,IA_IntPoint> {
d97 2
a98 2
    IA_IntPoint  operator() (const IA_IntPoint &ip) const {
	IA_IntPoint rmip = ip;
d103 1
a103 1
    IA_Closure<IA_IntPoint,IA_IntPoint> *clone_self() const {
d118 1
a118 1
class F_i: public IA_Closure<IA_IntPoint,IA_IntPoint> {
d127 2
a128 2
    IA_IntPoint  operator() (const IA_IntPoint &ip) const {
	IA_IntPoint fiip = ip;
d134 1
a134 1
    IA_Closure<IA_IntPoint,IA_IntPoint> *clone_self() const {
d145 1
a145 1
class G_i: public IA_Closure<IA_IntPoint,IA_IntPoint> {
d154 2
a155 2
    IA_IntPoint  operator() (const IA_IntPoint &ip) const {
	IA_IntPoint giip = ip;
d161 1
a161 1
    IA_Closure<IA_IntPoint,IA_IntPoint> *clone_self() const {
d172 1
a172 1
class  w_image: public IA_ClosureI<IA_IntPoint,float> {
d177 2
a178 2
    w_image( IA_IntPointSet ps, int i, int whichDimen )
    :IA_ClosureI<IA_IntPoint,float>(ps), w_i(i), d(whichDimen) {};
d180 2
a181 2
    float operator() (const IA_IntPoint& xy) const {
	int evenodd = (idiv2(xy[d],w_i-1) % 2); // 0 if even
d186 1
a186 1
    IA_ClosureI<IA_IntPoint,float> *clone_self() const {
d198 2
a199 2
IA_FloatDiscreteImage
walsh(const IA_FloatDiscreteImage & image) {
d201 1
a201 1
    IA_FloatDiscreteImage  a(image);     // Make a local copy
d205 1
a205 1
    IA_IntPoint a_dim_width( a.domain().sup() - a.domain().inf() + 1 );
d211 1
a211 1
    IA_IntPointSet  aps = a.domain();
d215 1
a215 1
	a = compose(a,rm).restrict(aps);  // permutate rows of the image
d225 1
a225 7
	    // We want to say: a = compose(a,gi) + w*compose(a,fi);

	    // However, compose() creates a comprehensive image,
	    // it must be restricted to make it an extensive image.

	    a = compose(a,gi).restrict(aps)
		+ w*compose(a,fi).restrict(aps);
d235 4
a238 3
IA_FloatDiscreteImage
walsh_inv(const IA_FloatDiscreteImage & a) {
    IA_IntPoint a_dim_width( a.domain().sup() - a.domain().inf() + 1 );
d240 1
a240 1
    double nm = coordProduct( a_dim_width );
@
