001/* ***** BEGIN LICENSE BLOCK *****
002 * 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
006 * the License. You may obtain a copy of the License at
007 * http://www.mozilla.org/MPL/
008 *
009 * Software distributed under the License is distributed on an "AS IS" basis,
010 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
011 * for the specific language governing rights and limitations under the
012 * License.
013 *
014 * The Original Code is JTransforms.
015 *
016 * The Initial Developer of the Original Code is
017 * Piotr Wendykier, Emory University.
018 * Portions created by the Initial Developer are Copyright (C) 2007-2009
019 * the Initial Developer. All Rights Reserved.
020 *
021 * Alternatively, the contents of this file may be used under the terms of
022 * either the GNU General Public License Version 2 or later (the "GPL"), or
023 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
024 * in which case the provisions of the GPL or the LGPL are applicable instead
025 * of those above. If you wish to allow use of your version of this file only
026 * under the terms of either the GPL or the LGPL, and not to allow others to
027 * use your version of this file under the terms of the MPL, indicate your
028 * decision by deleting the provisions above and replace them with the notice
029 * and other provisions required by the GPL or the LGPL. If you do not delete
030 * the provisions above, a recipient may use your version of this file under
031 * the terms of any one of the MPL, the GPL or the LGPL.
032 *
033 * ***** END LICENSE BLOCK ***** */
034
035package edu.emory.mathcs.jtransforms.dct;
036
037import edu.emory.mathcs.utils.IOUtils;
038
039/**
040 * Accuracy check of single precision DCT's
041 * 
042 * @author Piotr Wendykier (piotr.wendykier@gmail.com)
043 * 
044 */
045public class AccuracyCheckFloatDCT {
046
047    private static int[] sizes1D = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 120, 128, 256, 310, 512, 1024, 1056, 2048, 8192, 10158, 16384, 32768, 65536, 131072 };
048
049    private static int[] sizes2D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 120, 128, 256, 310, 511, 512, 1024 };
050
051    private static int[] sizes3D = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 32, 64, 100, 128 };
052
053    private static double eps = Math.pow(2, -23);
054
055    private AccuracyCheckFloatDCT() {
056
057    }
058
059    public static void checkAccuracyDCT_1D() {
060        System.out.println("Checking accuracy of 1D DCT...");
061        for (int i = 0; i < sizes1D.length; i++) {
062            FloatDCT_1D dct = new FloatDCT_1D(sizes1D[i]);
063            double err = 0;
064            float[] a = new float[sizes1D[i]];
065            IOUtils.fillMatrix_1D(sizes1D[i], a);
066            float[] b = new float[sizes1D[i]];
067            IOUtils.fillMatrix_1D(sizes1D[i], b);
068            dct.forward(a, true);
069            dct.inverse(a, true);
070            err = computeRMSE(a, b);
071            if (err > eps) {
072                System.err.println("\tsize = " + sizes1D[i] + ";\terror = " + err);
073            } else {
074                System.out.println("\tsize = " + sizes1D[i] + ";\terror = " + err);
075            }
076            a = null;
077            b = null;
078            dct = null;
079            System.gc();
080        }
081    }
082
083    public static void checkAccuracyDCT_2D() {
084        System.out.println("Checking accuracy of 2D DCT (float[] input)...");
085        for (int i = 0; i < sizes2D.length; i++) {
086            FloatDCT_2D dct2 = new FloatDCT_2D(sizes2D[i], sizes2D[i]);
087            double err = 0;
088            float[] a = new float[sizes2D[i] * sizes2D[i]];
089            IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a);
090            float[] b = new float[sizes2D[i] * sizes2D[i]];
091            IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], b);
092            dct2.forward(a, true);
093            dct2.inverse(a, true);
094            err = computeRMSE(a, b);
095            if (err > eps) {
096                System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err);
097            } else {
098                System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err);
099            }
100            a = null;
101            b = null;
102            dct2 = null;
103            System.gc();
104        }
105        System.out.println("Checking accuracy of 2D DCT (float[][] input)...");
106        for (int i = 0; i < sizes2D.length; i++) {
107            FloatDCT_2D dct2 = new FloatDCT_2D(sizes2D[i], sizes2D[i]);
108            double err = 0;
109            float[][] a = new float[sizes2D[i]][sizes2D[i]];
110            IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], a);
111            float[][] b = new float[sizes2D[i]][sizes2D[i]];
112            IOUtils.fillMatrix_2D(sizes2D[i], sizes2D[i], b);
113            dct2.forward(a, true);
114            dct2.inverse(a, true);
115            err = computeRMSE(a, b);
116            if (err > eps) {
117                System.err.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err);
118            } else {
119                System.out.println("\tsize = " + sizes2D[i] + " x " + sizes2D[i] + ";\terror = " + err);
120            }
121            a = null;
122            b = null;
123            dct2 = null;
124            System.gc();
125        }
126
127    }
128
129    public static void checkAccuracyDCT_3D() {
130        System.out.println("Checking accuracy of 3D DCT (float[] input)...");
131        for (int i = 0; i < sizes3D.length; i++) {
132            FloatDCT_3D dct3 = new FloatDCT_3D(sizes3D[i], sizes3D[i], sizes3D[i]);
133            double err = 0;
134            float[] a = new float[sizes3D[i] * sizes3D[i] * sizes3D[i]];
135            IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a);
136            float[] b = new float[sizes3D[i] * sizes3D[i] * sizes3D[i]];
137            IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], b);
138            dct3.forward(a, true);
139            dct3.inverse(a, true);
140            err = computeRMSE(a, b);
141            if (err > eps) {
142                System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err);
143            } else {
144                System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err);
145            }
146            a = null;
147            b = null;
148            dct3 = null;
149            System.gc();
150        }
151
152        System.out.println("Checking accuracy of 3D DCT (float[][][] input)...");
153        for (int i = 0; i < sizes3D.length; i++) {
154            FloatDCT_3D dct3 = new FloatDCT_3D(sizes3D[i], sizes3D[i], sizes3D[i]);
155            double err = 0;
156            float[][][] a = new float[sizes3D[i]][sizes3D[i]][sizes3D[i]];
157            IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], a);
158            float[][][] b = new float[sizes3D[i]][sizes3D[i]][sizes3D[i]];
159            IOUtils.fillMatrix_3D(sizes3D[i], sizes3D[i], sizes3D[i], b);
160            dct3.forward(a, true);
161            dct3.inverse(a, true);
162            err = computeRMSE(a, b);
163            if (err > eps) {
164                System.err.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err);
165            } else {
166                System.out.println("\tsize = " + sizes3D[i] + " x " + sizes3D[i] + " x " + sizes3D[i] + ";\t\terror = " + err);
167            }
168            a = null;
169            b = null;
170            dct3 = null;
171            System.gc();
172        }
173    }
174
175    private static double computeRMSE(float[] a, float[] b) {
176        if (a.length != b.length) {
177            throw new IllegalArgumentException("Arrays are not the same size");
178        }
179        double rms = 0;
180        double tmp;
181        for (int i = 0; i < a.length; i++) {
182            tmp = (a[i] - b[i]);
183            rms += tmp * tmp;
184        }
185        return Math.sqrt(rms / (float) a.length);
186    }
187
188    private static double computeRMSE(float[][] a, float[][] b) {
189        if (a.length != b.length || a[0].length != b[0].length) {
190            throw new IllegalArgumentException("Arrays are not the same size");
191        }
192        double rms = 0;
193        double tmp;
194        for (int r = 0; r < a.length; r++) {
195            for (int c = 0; c < a[0].length; c++) {
196                tmp = (a[r][c] - b[r][c]);
197                rms += tmp * tmp;
198            }
199        }
200        return Math.sqrt(rms / (a.length * a[0].length));
201    }
202
203    private static double computeRMSE(float[][][] a, float[][][] b) {
204        if (a.length != b.length || a[0].length != b[0].length || a[0][0].length != b[0][0].length) {
205            throw new IllegalArgumentException("Arrays are not the same size");
206        }
207        double rms = 0;
208        double tmp;
209        for (int s = 0; s < a.length; s++) {
210            for (int r = 0; r < a[0].length; r++) {
211                for (int c = 0; c < a[0][0].length; c++) {
212                    tmp = (a[s][r][c] - b[s][r][c]);
213                    rms += tmp * tmp;
214                }
215            }
216        }
217        return Math.sqrt(rms / (a.length * a[0].length * a[0][0].length));
218    }
219
220    public static void main(String[] args) {
221        checkAccuracyDCT_1D();
222        checkAccuracyDCT_2D();
223        checkAccuracyDCT_3D();
224        System.exit(0);
225    }
226}