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


1.4
date	95.01.06.05.28.07;	author rjj;	state Exp;
branches;
next	1.3;

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

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

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


desc
@Iterative version of N-dimension fft.
@


1.4
log
@Changed name of f(), g(), w_image() functions/class to F, G, W,
  and inlined their usage.
@
text
@// fft-nd-iter.c        Section: 8.3  Fast Fourier Transform
//      Copyright 1994, Center for Computer Vision and Visualization,
//      University of Florida.  All rights reserved.
// Image Algebra C++ implementation of the N-dimensional Fast
//       Fourier Transform using spatial transforms.
// This version is implemented with an iterative algorithm for
//       the F() and G() functions, rather than using compose().
//

#include "fft.h"
#include "PSIter.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];
}



//==============================================================
//
// FFT 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;
}


////////////////////////////////////////////////////////////////////
// permute_rows - uses p function above to permute the
//  dimension  whichDimen  rows.

IA_Image<IA_Point<int>,complex>
permute_rows(IA_Image<IA_Point<int>,complex> img, int whichDimen) {

    int dWidth = pointsetDimenWidth(img.domain(), whichDimen);
    IA_Image<IA_Point<int>,complex> tmp(img.domain(),0.0);
    IA_PSIter<IA_Point<int> > psi(img.domain());

    IA_Point<int> ip, pip;
    while(psi(ip)) {
	pip = ip;
	pip[whichDimen] = rho( pip[whichDimen], dWidth );
	tmp[ip] = img( pip );
    }
    return tmp;
}


////////////////////////////////////////////////////////////////////
// F -  An iterative version of the F() function.
//     It handles the even indices.

IA_Image<IA_Point<int>,complex>
F (int i, IA_Image<IA_Point<int>,complex> img, int d) {

    IA_Image<IA_Point<int>,complex> tmp( img );
    IA_PSIter<IA_Point<int> > psi(img.domain());

    IA_Point<int> ip, fip;
    int evenodd, ipow2i_1(ipow2(i-1));

    while(psi(ip)) {
	fip = ip;
	evenodd = idiv2( ip[d], i-1 ) % 2;
	if(evenodd == 0) {
	    fip[d] += ipow2i_1;
	    tmp[ip] = img( fip );
	}
    }
    return tmp;
}


////////////////////////////////////////////////////////////////////
// G -  An iterative version of the G() function.
//     It handles the odd indices cases.

IA_Image<IA_Point<int>,complex>
G (int i,IA_Image<IA_Point<int>,complex> img, int d) {

    IA_Image<IA_Point<int>,complex> tmp(img);
    IA_PSIter<IA_Point<int> > psi(img.domain());

    IA_Point<int> ip, gip;
    int evenodd, ipow2i_1(ipow2(i-1));

    while(psi(ip)) {
	gip = ip;
	evenodd = idiv2(ip[d], i-1) % 2;
	if(evenodd) {
	    gip[d] -= ipow2i_1;
	    tmp[ip] = img( gip );
	}
    }
    return tmp;
}


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

class W: public IA_ClosureI<IA_Point<int>,complex> {
private:
    int         w_i;
    int         u;
    double      piU;
    int         d;
public:
    W( int i, int whichDimen, IA_Set<IA_Point<int> > ps )
    :IA_ClosureI<IA_Point<int>,complex>(ps), 
     w_i(i-1), u(ipow2(i-1)), d(whichDimen) {
	 piU = PI/double(u);
     }

    complex operator() (const IA_Point<int>& xy) const {
        int evenodd( idiv2(xy[d],w_i) % 2 );
	double v = piU * double(xy[d] % u);

	if(evenodd) {
	    return  complex( -cos(v), -sin(v) );
	} else {
	    return  complex(  cos(v),  sin(v) );
	}
    }

    IA_ClosureI<IA_Point<int>,complex> *clone_self() const {
        return new W(*this);
    }
};


////////////////////////////////////////////////////////////////////
// fft -  Calcualte the Fast Fourier Transform.
// This version utilizes iterative versions of the f() and g()
//   functions, it does not utilize composition.

IA_Image<IA_Point<int>,complex>
fft(const IA_Image<IA_Point<int>,complex> & image) {

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

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

    // Note: By adding the dimension variable  d  we are able to
    //       avoid the text solution of applying a transpose to
    //       the image for the second dimension.  Also, it makes
    //       code more general, since it now is applicable to
    //       images over and dimension point set.
    for ( int d=0; d < a.domain().dim(); d++ ) {
	a = permute_rows(a, d);

	// The following is the acutal FFT code for single dimension.
	for( int i = 1;
	     i <= ilog2( pointsetDimenWidth(aps,d) );
	     i++ ) {

	    
	    a = G(i, a, d)  +  W(i, d, aps) * F(i, a, d);
	}
    }
    return a;
}


////////////////////////////////////////////////////////////////////
// fft_inv - Calculate the inverse FFT image using the forward FFT

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

    double nm = product( a_dim_width );  // area of point set

    return conj( fft( conj(a) )/nm );
}
@


1.3
log
@Added algorithm name to header.
Added return 0; to test files.
@
text
@d4 4
a7 3
// Image Algebra C++ implementation of the ND fast
//       Fourier transform using spatial transforms.
// This version is implemented with an iterative algorithm.
a17 2


d69 1
a69 2
int
rho(int i, int n) {
d103 1
a103 1
// f -  An iterative version of the f() function.
d107 1
a107 1
f (int i, IA_Image<IA_Point<int>,complex> img, int d) {
d128 1
a128 1
// g -  An iterative version of the g() function.
d132 1
a132 1
g(int i,IA_Image<IA_Point<int>,complex> img, int d) {
d153 1
a153 1
// class w_image -
d157 1
a157 1
class w_image: public IA_ClosureI<IA_Point<int>,complex> {
d164 1
a164 1
    w_image( IA_Set<IA_Point<int> > ps, int i, int whichDimen )
d182 1
a182 1
        return new w_image(*this);
d215 2
a216 3
	    w_image w(aps, i, d);

	    a = g(i, a, d) + w * f(i, a, d);
@


1.2
log
@Sanitized IA type names to be v1.0 compliant.
Added/modifed to have standard header.
@
text
@d1 1
a1 1
// fft-nd-iter.c    Section: 8.3
a7 2

static char rcsid[] = "$Id: fft-nd-iter.c,v 1.1 1994/02/20 04:05:09 rjj Exp rjj $";
@


1.1
log
@Initial revision
@
text
@d1 6
a6 4
// Emacs: -*- C++ -*-
// fft-nd-iter.c   Image Algebra C++ implementation of the ND fast
//         Fourier transform using spatial transforms.
//         This version is implemented with an iterative algorithm.
a7 1
// Cross-references: IA with Apps. pp. x-X; IA Handbook pp. y-Y
d9 1
a9 1
static char rcsid[] = "$Id$";
d11 3
a13 1
#include "imagealgebra.h"
a14 1
typedef IA_PSIter<IA_IntPoint> IA_IPSIter;
a15 1

d55 1
a55 1
pointsetDimenWidth(const IA_IntPointSet &ips, int dimen) {
a59 25
////////////////////////////////////////////////////////////////////
// coordSum -  Calculate the sum of a point\'s coordinates.

int
coordSum(const IA_IntPoint &ip) {
    int ipsum = 0;
    for ( int d=0; d < ip.dim(); d++ ) {
	ipsum += ip[d];
    }
    return ipsum;
}


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

a68 22
// conjugate - Calculate the conjuate of a complex-valued image.

IA_ComplexDiscreteImage
conjugate(const IA_ComplexDiscreteImage &img) {

#if IA_MAJOR_VERSION_NUMBER<1 && IA_MINOR_VERSION_NUMBER<10
    IA_ComplexDiscreteImage tmp(img.domain(),complex(0.0,0.0));
    IA_IPSIter psi(img.domain());

    IA_IntPoint ip;
    while(psi(ip)) {
	tmp[ip] = conj(img(ip));
    }
    return tmp;
#else
    // iac++ library now supports this directly.
    return conj(img);
#endif
}


////////////////////////////////////////////////////////////////////
d82 1
a82 1
} //end p
d89 2
a90 2
IA_ComplexDiscreteImage
permute_rows(IA_ComplexDiscreteImage img, int whichDimen) {
d93 2
a94 2
    IA_ComplexDiscreteImage tmp(img.domain(),0.0);
    IA_IPSIter psi(img.domain());
d96 1
a96 1
    IA_IntPoint ip, pip;
d108 1
a108 1
//     It handles the odd indices.
d110 2
a111 2
IA_ComplexDiscreteImage
f (int i, IA_ComplexDiscreteImage img, int d) {
d113 2
a114 2
    IA_ComplexDiscreteImage tmp( img );
    IA_IPSIter psi(img.domain());
d116 1
a116 1
    IA_IntPoint ip, fip;
d135 2
a136 2
IA_ComplexDiscreteImage
g(int i,IA_ComplexDiscreteImage img, int d) {
d138 2
a139 2
    IA_ComplexDiscreteImage tmp(img);
    IA_IPSIter psi(img.domain());
d141 1
a141 1
    IA_IntPoint ip, gip;
d161 1
a161 1
class w_image: public IA_ClosureI<IA_IntPoint,complex> {
d168 3
a170 3
    w_image( IA_IntPointSet ps, int i, int whichDimen )
    :IA_ClosureI<IA_IntPoint,complex>(ps), 
     w_i(i), u(ipow2(i-1)), d(whichDimen) {
d174 2
a175 2
    complex operator() (const IA_IntPoint& xy) const {
        int evenodd( idiv2(xy[d],w_i-1) % 2 );
d185 1
a185 1
    IA_ClosureI<IA_IntPoint,complex> *clone_self() const {
d193 1
a193 1
//   This version utilizes iterative versions of the f() and g()
d196 2
a197 2
IA_ComplexDiscreteImage
fft(const IA_ComplexDiscreteImage & image) {
d199 1
a199 1
    IA_ComplexDiscreteImage a(image);  // Make a local copy
d204 1
a204 1
    IA_IntPointSet  aps = a.domain();
d206 5
d219 1
a219 1
	    w_image w(a.domain(), i, d);
d231 3
a233 3
IA_ComplexDiscreteImage
fft_inv(const IA_ComplexDiscreteImage & a) {
    IA_IntPoint a_dim_width( a.domain().sup() - a.domain().inf() + 1 );
d235 1
a235 1
    double nm = coordProduct( a_dim_width );  // area of point set
d237 1
a237 1
    return conjugate( fft( conjugate(a) )/nm );
@
