001/*
002 * ***** BEGIN LICENSE BLOCK ***** Version: MPL 1.1/GPL 2.0/LGPL 2.1
003 * 
004 * The contents of this file are subject to the Mozilla Public License Version
005 * 1.1 (the "License"); you may not use this file except in compliance with the
006 * License. You may obtain a copy of the License at http://www.mozilla.org/MPL/
007 * 
008 * Software distributed under the License is distributed on an "AS IS" basis,
009 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License for
010 * the specific language governing rights and limitations under the License.
011 * 
012 * The Original Code is JTransforms.
013 * 
014 * The Initial Developer of the Original Code is Piotr Wendykier, Emory
015 * University. Portions created by the Initial Developer are Copyright (C)
016 * 2007-2009 the Initial Developer. All Rights Reserved.
017 * 
018 * Alternatively, the contents of this file may be used under the terms of
019 * either the GNU General Public License Version 2 or later (the "GPL"), or the
020 * GNU Lesser General Public License Version 2.1 or later (the "LGPL"), in which
021 * case the provisions of the GPL or the LGPL are applicable instead of those
022 * above. If you wish to allow use of your version of this file only under the
023 * terms of either the GPL or the LGPL, and not to allow others to use your
024 * version of this file under the terms of the MPL, indicate your decision by
025 * deleting the provisions above and replace them with the notice and other
026 * provisions required by the GPL or the LGPL. If you do not delete the
027 * provisions above, a recipient may use your version of this file under the
028 * terms of any one of the MPL, the GPL or the LGPL.
029 * 
030 * ***** END LICENSE BLOCK *****
031 */
032
033package edu.emory.mathcs.jtransforms.fft;
034
035// @formatter:off
036/**
037 * <p>
038 * This is a set of utility methods for R/W access to data resulting from a call
039 * to the Fourier transform of <em>real</em> data. Memory optimized methods,
040 * namely
041 * <ul>
042 *   <li>{@link DoubleFFT_3D#realForward(double[])}</li>
043 *   <li>{@link DoubleFFT_3D#realForward(double[][][])}</li>
044 *   <li>{@link FloatFFT_3D#realForward(float[])}</li>
045 *   <li>{@link FloatFFT_3D#realForward(float[][][])}</li>
046 * </ul>
047 * are implemented to handle this case specifically. However, packing of the
048 * data in the data array is somewhat obscure. This class provides methods for
049 * direct access to the data, without the burden of all necessary tests.
050 * </p>
051 * <h3>Example for Fourier Transform of real, double precision 1d data</h3>
052 * <p>
053 * <pre>
054 *   DoubleFFT_3D fft = new DoubleFFT_2D(slices, rows, columns);
055 *   double[] data = new double[2 * slices * rows * columns];
056 *   ...
057 *   fft.realForwardFull(data);
058 *   data[(s1 * rows + r1) * 2 * columns + c1] = val1;
059 *   val2 = data[(s2 * rows + r2) * 2 * columns + c2];
060 * </pre>
061 * is equivalent to
062 * <pre>
063 *   DoubleFFT_3D fft = new DoubleFFT_3D(slices, rows, columns);
064 *   RealFFTUtils_3D unpacker = new RealFFTUtils_3D(slices, rows, columns);
065 *   double[] data = new double[slices * rows * columns];
066 *   ...
067 *   fft.realForward(data);
068 *   unpacker.pack(val1, s1, r1, c1, data);
069 *   val2 = unpacker.unpack(s2, r2, c2, data, 0);
070 * </pre>
071 * Even (resp. odd) values of <code>c</code> correspond to the real (resp.
072 * imaginary) part of the Fourier mode.
073 * </p>
074 * <h3>Example for Fourier Transform of real, double precision 3d data</h3>
075 * <p>
076 * <pre>
077 *   DoubleFFT_3D fft = new DoubleFFT_3D(slices, rows, columns);
078 *   double[][][] data = new double[slices][rows][2 * columns];
079 *   ...
080 *   fft.realForwardFull(data);
081 *   data[s1][r1][c1] = val1;
082 *   val2 = data[s2][r2][c2];
083 * </pre>
084 * is equivalent to
085 * <pre>
086 *   DoubleFFT_3D fft = new DoubleFFT_3D(slices, rows, columns);
087 *   RealFFTUtils_3D unpacker = new RealFFTUtils_3D(slices, rows, columns);
088 *   double[][][] data = new double[slices][rows][columns];
089 *   ...
090 *   fft.realForward(data);
091 *   unpacker.pack(val1, s1, r1, c1, data);
092 *   val2 = unpacker.unpack(s2, r2, c2, data, 0);
093 * </pre>
094 * Even (resp. odd) values of <code>c</code> correspond to the real (resp.
095 * imaginary) part of the Fourier mode.
096 * </p>
097 * 
098 * @author S&eacute;bastien Brisard
099 * 
100 */
101// @formatter:on
102public class RealFFTUtils_3D {
103    /** The constant <code>int</code> value of 1. */
104    private static final int ONE = 1;
105
106    /** The constant <code>int</code> value of 2. */
107    private static final int TWO = 2;
108
109    /** The constant <code>int</code> value of 0. */
110    private static final int ZERO = 0;
111
112    /** The size of the data in the third direction. */
113    private final int columns;
114
115    /** The size of the data in the second direction. */
116    private final int rows;
117
118    /** The constant value of <code>2 * columns</code>. */
119    private final int rowStride;
120
121    /** The size of the data in the first direction. */
122    private final int slices;
123
124    /** The constant value of <code>2 * rows * columns</code>. */
125    private final int sliceStride;
126
127    /**
128     * Creates a new instance of this class. The size of the underlying
129     * {@link DoubleFFT_3D} or {@link FloatFFT_3D} must be specified.
130     * 
131     * @param slices
132     *            number of slices
133     * @param rows
134     *            number of rows
135     * @param columns
136     *            number of columns
137     */
138    public RealFFTUtils_3D(final int slices, final int rows, final int columns) {
139        this.slices = slices;
140        this.rows = rows;
141        this.columns = columns;
142        this.rowStride = columns;
143        this.sliceStride = rows * this.rowStride;
144    }
145
146    /**
147     * <p>
148     * Returns the 1d index of the specified 3d Fourier mode. In other words, if
149     * <code>packed</code> contains the transformed data following a call to
150     * {@link DoubleFFT_3D#realForwardFull(double[])} or
151     * {@link FloatFFT_3D#realForward(float[])}, then the returned value
152     * <code>index</code> gives access to the <code>[s][r][c]</code> Fourier
153     * mode
154     * <ul>
155     * <li>if <code>index == {@link Integer#MIN_VALUE}</code>, then the Fourier
156     * mode is zero,</li>
157     * <li>if <code>index >= 0</code>, then the Fourier mode is
158     * <code>packed[index]</code>,</li>
159     * <li>if <code>index < 0</code>, then the Fourier mode is
160     * <code>-packed[-index]</code>,</li>
161     * </ul>
162     * </p>
163     * 
164     * @param s
165     *            the slice index
166     * @param r
167     *            the row index
168     * @param c
169     *            the column index
170     * @return the value of <code>index</code>
171     */
172    public int getIndex(final int s, final int r, final int c) {
173        final int cmod2 = c & ONE;
174        final int rmul2 = r << ONE;
175        final int smul2 = s << ONE;
176        final int ss = s == ZERO ? ZERO : slices - s;
177        final int rr = r == ZERO ? ZERO : rows - r;
178        if (c <= ONE) {
179            if (r == ZERO) {
180                if (s == ZERO) {
181                    return c == ZERO ? ZERO : Integer.MIN_VALUE;
182                } else if (smul2 < slices) {
183                    return s * sliceStride + c;
184                } else if (smul2 > slices) {
185                    final int index = ss * sliceStride;
186                    return cmod2 == ZERO ? index : -(index + ONE);
187                } else {
188                    return cmod2 == ZERO ? s * sliceStride : Integer.MIN_VALUE;
189                }
190            } else if (rmul2 < rows) {
191                return s * sliceStride + r * rowStride + c;
192            } else if (rmul2 > rows) {
193                final int index = ss * sliceStride + rr * rowStride;
194                return cmod2 == ZERO ? index : -(index + ONE);
195            } else {
196                if (s == ZERO) {
197                    return cmod2 == ZERO ? r * rowStride : Integer.MIN_VALUE;
198                } else if (smul2 < slices) {
199                    return s * sliceStride + r * rowStride + c;
200                } else if (smul2 > slices) {
201                    final int index = ss * sliceStride + r * rowStride;
202                    return cmod2 == ZERO ? index : -(index + ONE);
203                } else {
204                    final int index = s * sliceStride + r * rowStride;
205                    return cmod2 == ZERO ? index : Integer.MIN_VALUE;
206                }
207            }
208        } else if (c < columns) {
209            return s * sliceStride + r * rowStride + c;
210        } else if (c > columns + ONE) {
211            final int cc = (columns << ONE) - c;
212            final int index = ss * sliceStride + rr * rowStride + cc;
213            return cmod2 == ZERO ? index : -(index + TWO);
214        } else {
215            if (r == ZERO) {
216                if (s == ZERO) {
217                    return cmod2 == ZERO ? ONE : Integer.MIN_VALUE;
218                } else if (smul2 < slices) {
219                    final int index = ss * sliceStride;
220                    return cmod2 == ZERO ? index + ONE : -index;
221                } else if (smul2 > slices) {
222                    final int index = s * sliceStride;
223                    return cmod2 == ZERO ? index + ONE : index;
224                } else {
225                    final int index = s * sliceStride;
226                    return cmod2 == ZERO ? index + ONE : Integer.MIN_VALUE;
227                }
228            } else if (rmul2 < rows) {
229                final int index = ss * sliceStride + rr * rowStride;
230                return cmod2 == ZERO ? index + ONE : -index;
231            } else if (rmul2 > rows) {
232                final int index = s * sliceStride + r * rowStride;
233                return cmod2 == ZERO ? index + ONE : index;
234            } else {
235                if (s == ZERO) {
236                    final int index = r * rowStride + ONE;
237                    return cmod2 == ZERO ? index : Integer.MIN_VALUE;
238                } else if (smul2 < slices) {
239                    final int index = ss * sliceStride + r * rowStride;
240                    return cmod2 == ZERO ? index + ONE : -index;
241                } else if (smul2 > slices) {
242                    final int index = s * sliceStride + r * rowStride;
243                    return cmod2 == ZERO ? index + ONE : index;
244                } else {
245                    final int index = s * sliceStride + r * rowStride;
246                    return cmod2 == ZERO ? index + ONE : Integer.MIN_VALUE;
247                }
248            }
249        }
250    }
251
252    /**
253     * Sets the specified Fourier mode of the transformed data. The data array
254     * results from a call to {@link DoubleFFT_3D#realForward(double[])}.
255     * 
256     * @param val
257     *            the new value of the <code>[s][r][c]</code> Fourier mode
258     * @param s
259     *            the slice index
260     * @param r
261     *            the row index
262     * @param c
263     *            the column index
264     * @param packed
265     *            the transformed data
266     * @param pos
267     *            index of the first element in array <code>packed</code>
268     */
269    public void pack(final double val, final int s, final int r, final int c,
270            final double[] packed, final int pos) {
271        final int i = getIndex(s, r, c);
272        if (i >= 0) {
273            packed[pos + i] = val;
274        } else if (i > Integer.MIN_VALUE) {
275            packed[pos - i] = -val;
276        } else {
277            throw new IllegalArgumentException(String.format(
278                    "[%d][%d][%d] component cannot be modified (always zero)",
279                    s, r, c));
280        }
281    }
282
283    /**
284     * Sets the specified Fourier mode of the transformed data. The data array
285     * results from a call to {@link DoubleFFT_3D#realForward(double[][][])}.
286     * 
287     * @param val
288     *            the new value of the <code>[s][r][c]</code> Fourier mode
289     * @param s
290     *            the slice index
291     * @param r
292     *            the row index
293     * @param c
294     *            the column index
295     * @param packed
296     *            the transformed data
297     */
298    public void pack(final double val, final int s, final int r, final int c,
299            final double[][][] packed) {
300        final int i = getIndex(s, r, c);
301        final int ii = Math.abs(i);
302        final int ss = ii / sliceStride;
303        final int remainder = ii % sliceStride;
304        final int rr = remainder / rowStride;
305        final int cc = remainder % rowStride;
306        if (i >= 0) {
307            packed[ss][rr][cc] = val;
308        } else if (i > Integer.MIN_VALUE) {
309            packed[ss][rr][cc] = -val;
310        } else {
311            throw new IllegalArgumentException(
312                    String.format(
313                            "[%d][%d] component cannot be modified (always zero)",
314                            r, c));
315        }
316    }
317
318    /**
319     * Sets the specified Fourier mode of the transformed data. The data array
320     * results from a call to {@link FloatFFT_3D#realForward(float[])}.
321     * 
322     * @param val
323     *            the new value of the <code>[s][r][c]</code> Fourier mode
324     * @param s
325     *            the slice index
326     * @param r
327     *            the row index
328     * @param c
329     *            the column index
330     * @param packed
331     *            the transformed data
332     * @param pos
333     *            index of the first element in array <code>packed</code>
334     */
335    public void pack(final float val, final int s, final int r, final int c,
336            final float[] packed, final int pos) {
337        final int i = getIndex(s, r, c);
338        if (i >= 0) {
339            packed[pos + i] = val;
340        } else if (i > Integer.MIN_VALUE) {
341            packed[pos - i] = -val;
342        } else {
343            throw new IllegalArgumentException(String.format(
344                    "[%d][%d][%d] component cannot be modified (always zero)",
345                    s, r, c));
346        }
347    }
348
349    /**
350     * Sets the specified Fourier mode of the transformed data. The data array
351     * results from a call to {@link FloatFFT_3D#realForward(float[][][])}.
352     * 
353     * @param val
354     *            the new value of the <code>[s][r][c]</code> Fourier mode
355     * @param s
356     *            the slice index
357     * @param r
358     *            the row index
359     * @param c
360     *            the column index
361     * @param packed
362     *            the transformed data
363     */
364    public void pack(final float val, final int s, final int r, final int c,
365            final float[][][] packed) {
366        final int i = getIndex(s, r, c);
367        final int ii = Math.abs(i);
368        final int ss = ii / sliceStride;
369        final int remainder = ii % sliceStride;
370        final int rr = remainder / rowStride;
371        final int cc = remainder % rowStride;
372        if (i >= 0) {
373            packed[ss][rr][cc] = val;
374        } else if (i > Integer.MIN_VALUE) {
375            packed[ss][rr][cc] = -val;
376        } else {
377            throw new IllegalArgumentException(String.format(
378                    "[%d][%d][%d] component cannot be modified (always zero)",
379                    s, r, c));
380        }
381    }
382
383    /**
384     * Returns the specified Fourier mode of the transformed data. The data
385     * array results from a call to {@link DoubleFFT_3D#realForward(double[])}.
386     * 
387     * @param s
388     *            the slice index
389     * @param r
390     *            the row index
391     * @param c
392     *            the column index
393     * @param packed
394     *            the transformed data
395     * @param pos
396     *            index of the first element in array <code>packed</code>
397     * @return the value of the <code>[s][r][c]</code> Fourier mode
398     */
399    public double unpack(final int s, final int r, final int c,
400            final double[] packed, final int pos) {
401        final int i = getIndex(s, r, c);
402        if (i >= 0) {
403            return packed[pos + i];
404        } else if (i > Integer.MIN_VALUE) {
405            return -packed[pos - i];
406        } else {
407            return ZERO;
408        }
409    }
410
411    /**
412     * Returns the specified Fourier mode of the transformed data. The data
413     * array results from a call to
414     * {@link DoubleFFT_3D#realForward(double[][][])} .
415     * 
416     * @param s
417     *            the slice index
418     * @param r
419     *            the row index
420     * @param c
421     *            the column index
422     * @param packed
423     *            the transformed data
424     * @return the value of the <code>[s][r][c]</code> Fourier mode
425     */
426    public double unpack(final int s, final int r, final int c,
427            final double[][][] packed) {
428        final int i = getIndex(s, r, c);
429        final int ii = Math.abs(i);
430        final int ss = ii / sliceStride;
431        final int remainder = ii % sliceStride;
432        final int rr = remainder / rowStride;
433        final int cc = remainder % rowStride;
434        if (i >= 0) {
435            return packed[ss][rr][cc];
436        } else if (i > Integer.MIN_VALUE) {
437            return -packed[ss][rr][cc];
438        } else {
439            return ZERO;
440        }
441    }
442
443    /**
444     * Returns the specified Fourier mode of the transformed data. The data
445     * array results from a call to {@link FloatFFT_3D#realForward(float[])} .
446     * 
447     * @param s
448     *            the slice index
449     * @param r
450     *            the row index
451     * @param c
452     *            the column index
453     * @param packed
454     *            the transformed data
455     * @param pos
456     *            index of the first element in array <code>packed</code>
457     * @return the value of the <code>[s][r][c]</code> Fourier mode
458     */
459    public float unpack(final int s, final int r, final int c,
460            final float[] packed, final int pos) {
461        final int i = getIndex(s, r, c);
462        if (i >= 0) {
463            return packed[pos + i];
464        } else if (i > Integer.MIN_VALUE) {
465            return -packed[pos - i];
466        } else {
467            return ZERO;
468        }
469    }
470
471    /**
472     * Returns the specified Fourier mode of the transformed data. The data
473     * array results from a call to {@link FloatFFT_3D#realForward(float[][][])}
474     * .
475     * 
476     * @param s
477     *            the slice index
478     * @param r
479     *            the row index
480     * @param c
481     *            the column index
482     * @param packed
483     *            the transformed data
484     * @return the value of the <code>[s][r][c]</code> Fourier mode
485     */
486    public float unpack(final int s, final int r, final int c,
487            final float[][][] packed) {
488        final int i = getIndex(s, r, c);
489        final int ii = Math.abs(i);
490        final int ss = ii / sliceStride;
491        final int remainder = ii % sliceStride;
492        final int rr = remainder / rowStride;
493        final int cc = remainder % rowStride;
494        if (i >= 0) {
495            return packed[ss][rr][cc];
496        } else if (i > Integer.MIN_VALUE) {
497            return -packed[ss][rr][cc];
498        } else {
499            return ZERO;
500        }
501    }
502}