Index: trunk/doc/pslib/psLibSDRS.tex =================================================================== --- trunk/doc/pslib/psLibSDRS.tex (revision 11598) +++ trunk/doc/pslib/psLibSDRS.tex (revision 11722) @@ -1,3 +1,3 @@ -%%% $Id: psLibSDRS.tex,v 1.440 2007-02-03 00:13:08 eugene Exp $ +%%% $Id: psLibSDRS.tex,v 1.441 2007-02-09 01:19:43 price Exp $ \documentclass[panstarrs,spec]{panstarrs} @@ -6435,76 +6435,69 @@ We require the ability to calculate the (fast) Fourier transforms of -floating-point one-dimensional vectors and two-dimensional images. We -expect that these will be implemented through wrapping an external -library. We define the following APIs to support FFT operations on vectors: - -\begin{prototype} -psVector *psVectorFFT(psVector *out, const psVector *in, psFFTFlags direction); -psVector *psVectorReal(psVector *out, const psVector *in); -psVector *psVectorImaginary(psVector *out, const psVector *in); -psVector *psVectorComplex(psVector *out, const psVector *real, const psVector *imag); -psVector *psVectorConjugate(psVector *out, const psVector *in); +floating-point one-dimensional vectors and two-dimensional images. +These will be implemented through wrapping an external library. We +expect that we shall only want to deal with purely real vectors and +images, rather than complex; this can save operations when performing +the FFT (roughly a factor of 2). We define the following APIs to +support FFT operations on vectors: + +\begin{prototype} +bool psVectorForwardFFT(psVector **real psVector **imag, const psVector *in); +bool psVectorBackwardFFT(psVector **out, const psVector *real, const psVector *imag, int origNum); psVector *psVectorPowerSpectrum(psVector *out, const psVector *in); -\end{prototype} - -The forward and reverse FFT is calculated using \code{psVectorFFT}, -which takes as input the vector of interest (\code{in}) and the -direction (\code{direction}), which is specified by an enumerated type -defined below. The input vector may be of type \code{psF32} or -\code{psC32}, the result is always \code{psC32}. If the input vector -is \code{psF32}, the direction must be forward. Neither the forward -or inverse transforms must multiply by $1/N$ (or $1/N^{1/2}$), and so -it falls to the responsibility of the user to multiply a vector that -has been forward- and reverse-transformed by $1/N$. - -Conversions between complex and real vectors requires the functions -\code{psVectorReal}, which returns the real part (\code{psF32}) of the -complex vector \code{in}, \code{psVectorImaginary}, which returns the -the magnitude of the imaginary part (\code{psF32}) of the complex -vector \code{in}, and \code{psVectorComplex}, which constructs a -complex vector (\code{psC32}) from the real and imaginary components -\code{real} and \code{imag}, both of type \code{psF32}. - -We also specify the additional utility functions -\code{psVectorConjugate} and \code{psVectorPowerSpectrum} which -construct the complex conjugate (\code{psC32}) and the magnitude -(\code{psF32}) vectors of the input complex vector (\code{psC32}). - -The direction of an FFT performed by \code{psVectorFFT} is specified -by an enumerated type, \code{psFFTFlags}: - -\begin{datatype} -/** Specify direction of FFT, and if the result is real or not */ -typedef enum { - PS_FFT_FORWARD = 1, ///< psImageFFT/psVectorFFT should perform a forward FFT. - PS_FFT_REVERSE = 2, ///< psImageFFT/psVectorFFT should perform a reverse FFT. - PS_FFT_REAL_RESULT = 4 ///< psImageFFT/psVectorFFT should return a real image. This is valid - ///< for only reverse FFT, i.e., the psImageFFT/psVectorFFT flag - ///< parameter should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT. -} psFFTFlags; -\end{datatype} - -The entry \code{PS_FFT_REAL_RESULT} means that the output of an -inverse FFT is to be purely real. An example of its use is: -\begin{verbatim} -out = psImageFFT(in, PS_FFT_REVERSE | PS_FFT_REAL_RESULT); -\end{verbatim} - -The output from a FFT must be of the same size as the input, and the -size of the power spectrum equal to $N/2 + 1$, where $N$ is the number -of elements in the input. In the future, if this adversely affects -performance, this will be revised so that the redundant information -will be neglected. +bool psVectorComplexMultiply(psVector **outReal, psVector **outImag, + const psVector *in1Real, const psVector *in1Imag, + const psVector *in2Real, const psVector *in2Imag); +\end{prototype} + +We might have chosen to use the \code{complex} types that C99 +provides, but unfortunately this isn't implemented on all the systems +that we would have desired. Instead, we use separate vectors +to carry the real and imaginary parts. + +The forward FFT (exponent $-1$) is calculated using +\code{psVectorForwardFFT}, which takes as input the vector of interest +(\code{in}) of type F32, returning the \code{real} and +\code{imaginary} parts in separate vectors, which are allocated if +required. + +The backward FFT (exponent $+1$) is calculated using +\code{psVectorBackwardFFT}, which takes as input the \code{real} and +\code{imaginary} parts, each of type F32, and returns the purely real +\code{out} vector, which is allocated if required. Due to the +ambiguity in the size of the original vector when the purely real +optimisation is made, this function also requires the number of +elements, \code{origNum} in the original vector before FFT-ing with +\code{psVectorForwardFFT}. + +Note that the FFTs are not normalised, and so it falls to the +responsibility of the user to multiply a vector that has been forward- +and reverse-transformed by $1/N$. + +\code{psVectorPowerSpectrum} calculates the power spectrum (magnitude +of the fourier transform) from the \code{in} vector (of type F32. + +Finally, for convenience, the function \code{psVectorComplexMultiply} +is provided to multiply two vectors of complex numbers which are +specified by their real and imaginary parts, and returning the real +and imaginary parts separately. In analogy with the vector FFT operations, we also define matching image operations as follows: \begin{prototype} -psImage *psImageFFT(psImage *out, const psImage *image, psFFTFlags direction); -psImage *psImageReal(psImage *out, const psImage *in); -psImage *psImageImaginary(psImage *out, const psImage *in); -psImage *psImageComplex(psImage *out, const psImage *real, const psImage *imag); -psImage *psImageConjugate(psImage *out, const psImage *in); +bool psImageForwardFFT(psImage **real psImage **imag, const psImage *in); +bool psImageBackwardFFT(psImage **out, const psImage *real, const psImage *imag, int origCols); psImage *psImagePowerSpectrum(psImage *out, const psImage *in); -\end{prototype} +bool psImageComplexMultiply(psImage **outReal, psImage **outImag, + const psImage *in1Real, const psImage *in1Imag, + const psImage *in2Real, const psImage *in2Imag); +\end{prototype} + +The only subtle difference is in \code{psImageBackwardFFT}, where the +number of columns in the original image (before applying +\code{psImageForwardFFT}) is required due to the ambiguity arising +from the optimisation made for purely real data. Note that the number +of rows need not be provided, because the FFT-ed image is shrunk in +the columns but not the rows. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -6612,21 +6605,25 @@ \subsubsection{Convolve an image with the kernel} -Given an input image and the convolution kernel, -\code{psImageConvolve} shall convolve the input image, -\code{in}, with the kernel, \code{kernel} and return the convolved -image, \code{out}. The API shall be the following: -\begin{prototype} -psImage *psImageConvolve(psImage *out, const psImage *in, - const psKernel *kernel, bool direct); -\end{prototype} - -Two methods shall be available for the convolution: if \code{direct} -is \code{true}, then the convolution shall be performed in real space -(appropriate for small kernels); otherwise, the convolution shall be -performed using Fast Fourier Transforms (FFTs; appropriate for larger -kernels). The latter option involves padding the input image, copying -the kernel into an image of the same size as the padded input image, -performing an FFT on each, multiplying the FFTs, and performing an -inverse FFT before trimming the image back to the original size. +Two methods are available for the convolution: direct and FFT. Given +an input image and the convolution kernel, +\code{psImageConvolveDirect} and \code{psImageConvolveFFT} shall +convolve the input image, \code{in}, with the kernel, \code{kernel} +and return the convolved image, \code{out}. The APIs shall be the +following: +\begin{prototype} +psImage *psImageConvolveDirect(psImage *out, const psImage *in, + const psKernel *kernel); +psImage *psImageConvolveFFT(psImage *out, const psImage *in, + const psKernel *kernel); +\end{prototype} + +\code{psImageConvolveDirect} shall perform the convolution in real +space (appropriate for small kernels); \code{psImageConvolveFFT} shall +perform the convolution using Fast Fourier Transforms (FFTs; +appropriate for larger kernels). The latter option involves padding +the input image, copying the kernel into an image of the same size as +the padded input image, performing an FFT on each, multiplying the +FFTs, and performing an inverse FFT before trimming the image back to +the original size. In the event that \code{out} is \code{NULL}, a new \code{psImage}