/* This file is part of jpegpixi, a program to interpolate pixels in JFIF image files. Copyright (C) 2003 Martin Dickopp Jpegpixi is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. Jpegpixi is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with jpegpixi; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. This file is in part derived from the JPEG library of the Independent JPEG Group, see the file `README.jpeglib' for copyright and license information. */ #ifdef HAVE_CONFIG_H # include #endif #include "util.h" #include "jpegpixi.h" /* Scaling factors. */ static const double scale [DCTSIZE] = { 1.0, 1.3870398453221474618, 1.3065629648763765279, 1.1758756024193587170, 1.0, 0.7856949583871021813, 0.5411961001461969844, 0.2758993792829430123 }; static JCOEF coeff (double x, UINT16 quantval) gcc_attr_const; /* Perform forward discrete cosine transformation and quantization. */ void fdct (JCOEF *freqs, const double *values, const UINT16 *quant) { int i; double tmp [DCTSIZE2], *ptr; ptr = tmp; for (i = 0; i < DCTSIZE; ++i) { double z2, z3, z4, z5, z11, z13; const double tmp0 = (values [0] + values [7]) / 64.0; const double tmp7 = (values [0] - values [7]) / 64.0; const double tmp1 = (values [1] + values [6]) / 64.0; const double tmp6 = (values [1] - values [6]) / 64.0; const double tmp2 = (values [2] + values [5]) / 64.0; const double tmp5 = (values [2] - values [5]) / 64.0; const double tmp3 = (values [3] + values [4]) / 64.0; const double tmp4 = (values [3] - values [4]) / 64.0; double tmp10 = tmp0 + tmp3; double tmp13 = tmp0 - tmp3; double tmp11 = tmp1 + tmp2; double tmp12 = tmp1 - tmp2; const double z1 = (tmp12 + tmp13) * 0.7071067811865475244; ptr [0] = tmp10 + tmp11; ptr [4] = tmp10 - tmp11; ptr [2] = tmp13 + z1; ptr [6] = tmp13 - z1; tmp10 = tmp4 + tmp5; tmp11 = tmp5 + tmp6; tmp12 = tmp6 + tmp7; z5 = (tmp10 - tmp12) * 0.3826834323650897717; z2 = 0.5411961001461969844 * tmp10 + z5; z4 = 1.3065629648763765279 * tmp12 + z5; z3 = 0.7071067811865475244 * tmp11; z11 = tmp7 + z3; z13 = tmp7 - z3; ptr [5] = z13 + z2; ptr [3] = z13 - z2; ptr [1] = z11 + z4; ptr [7] = z11 - z4; values += DCTSIZE; ptr += DCTSIZE; } ptr = tmp; for (i = 0; i < DCTSIZE; ++i) { double z2, z3, z4, z5, z11, z13; const double tmp0 = ptr [DCTSIZE * 0] + ptr [DCTSIZE * 7]; const double tmp7 = ptr [DCTSIZE * 0] - ptr [DCTSIZE * 7]; const double tmp1 = ptr [DCTSIZE * 1] + ptr [DCTSIZE * 6]; const double tmp6 = ptr [DCTSIZE * 1] - ptr [DCTSIZE * 6]; const double tmp2 = ptr [DCTSIZE * 2] + ptr [DCTSIZE * 5]; const double tmp5 = ptr [DCTSIZE * 2] - ptr [DCTSIZE * 5]; const double tmp3 = ptr [DCTSIZE * 3] + ptr [DCTSIZE * 4]; const double tmp4 = ptr [DCTSIZE * 3] - ptr [DCTSIZE * 4]; double tmp10 = tmp0 + tmp3; double tmp13 = tmp0 - tmp3; double tmp11 = tmp1 + tmp2; double tmp12 = tmp1 - tmp2; const double z1 = (tmp12 + tmp13) * 0.7071067811865475244; freqs [DCTSIZE * 0] = coeff ((tmp10 + tmp11) / (scale [i]), quant [DCTSIZE * 0]); freqs [DCTSIZE * 4] = coeff ((tmp10 - tmp11) / (scale [i]), quant [DCTSIZE * 4]); freqs [DCTSIZE * 2] = coeff ((tmp13 + z1) / (scale [i] * scale [2]), quant [DCTSIZE * 2]); freqs [DCTSIZE * 6] = coeff ((tmp13 - z1) / (scale [i] * scale [6]), quant [DCTSIZE * 6]); tmp10 = tmp4 + tmp5; tmp11 = tmp5 + tmp6; tmp12 = tmp6 + tmp7; z5 = (tmp10 - tmp12) * 0.3826834323650897717; z2 = 0.5411961001461969844 * tmp10 + z5; z4 = 1.3065629648763765279 * tmp12 + z5; z3 = 0.7071067811865475244 * tmp11; z11 = tmp7 + z3; z13 = tmp7 - z3; freqs [DCTSIZE * 5] = coeff ((z13 + z2) / (scale [i] * scale [5]), quant [DCTSIZE * 5]); freqs [DCTSIZE * 3] = coeff ((z13 - z2) / (scale [i] * scale [3]), quant [DCTSIZE * 3]); freqs [DCTSIZE * 1] = coeff ((z11 + z4) / (scale [i] * scale [1]), quant [DCTSIZE * 1]); freqs [DCTSIZE * 7] = coeff ((z11 - z4) / (scale [i] * scale [7]), quant [DCTSIZE * 7]); ++ptr; ++freqs; ++quant; } } /* Perform inverse discrete cosine transformation and dequantization. */ void idct (double *values, const JCOEF *freqs, const UINT16 *quant) { int i; double tmp [DCTSIZE2], *ptr; ptr = tmp; for (i = 0; i < DCTSIZE; ++i) { if (freqs [DCTSIZE * 1] == 0 && freqs [DCTSIZE * 2] && freqs [DCTSIZE * 3] && freqs [DCTSIZE * 4] && freqs [DCTSIZE * 5] && freqs [DCTSIZE * 6] && freqs [DCTSIZE * 6]) { const double dc_value = (double)freqs [DCTSIZE * 0] * (double)quant [DCTSIZE * 0] * scale [i]; ptr [DCTSIZE * 0] = dc_value; ptr [DCTSIZE * 1] = dc_value; ptr [DCTSIZE * 2] = dc_value; ptr [DCTSIZE * 3] = dc_value; ptr [DCTSIZE * 4] = dc_value; ptr [DCTSIZE * 5] = dc_value; ptr [DCTSIZE * 6] = dc_value; ptr [DCTSIZE * 7] = dc_value; } else { double tmp0 = (double)freqs [DCTSIZE * 0] * (double)quant [DCTSIZE * 0] * scale [i]; double tmp4 = (double)freqs [DCTSIZE * 1] * (double)quant [DCTSIZE * 1] * scale [i] * scale [1]; double tmp1 = (double)freqs [DCTSIZE * 2] * (double)quant [DCTSIZE * 2] * scale [i] * scale [2]; double tmp5 = (double)freqs [DCTSIZE * 3] * (double)quant [DCTSIZE * 3] * scale [i] * scale [3]; double tmp2 = (double)freqs [DCTSIZE * 4] * (double)quant [DCTSIZE * 4] * scale [i]; double tmp6 = (double)freqs [DCTSIZE * 5] * (double)quant [DCTSIZE * 5] * scale [i] * scale [5]; double tmp3 = (double)freqs [DCTSIZE * 6] * (double)quant [DCTSIZE * 6] * scale [i] * scale [6]; double tmp7 = (double)freqs [DCTSIZE * 7] * (double)quant [DCTSIZE * 7] * scale [i] * scale [7]; double tmp10 = tmp0 + tmp2; double tmp11 = tmp0 - tmp2; double tmp13 = tmp1 + tmp3; double tmp12 = (tmp1 - tmp3) * 1.4142135623730950488 - tmp13; const double z13 = tmp6 + tmp5; const double z10 = tmp6 - tmp5; const double z11 = tmp4 + tmp7; const double z12 = tmp4 - tmp7; const double z5 = (z10 + z12) * 1.8477590650225735123; tmp0 = tmp10 + tmp13; tmp3 = tmp10 - tmp13; tmp1 = tmp11 + tmp12; tmp2 = tmp11 - tmp12; tmp7 = z11 + z13; tmp11 = (z11 - z13) * 1.4142135623730950488; tmp10 = 1.0823922002923939688 * z12 - z5; tmp12 = -2.6131259297527530557 * z10 + z5; tmp6 = tmp12 - tmp7; tmp5 = tmp11 - tmp6; tmp4 = tmp10 + tmp5; ptr [DCTSIZE * 0] = tmp0 + tmp7; ptr [DCTSIZE * 1] = tmp1 + tmp6; ptr [DCTSIZE * 2] = tmp2 + tmp5; ptr [DCTSIZE * 3] = tmp3 - tmp4; ptr [DCTSIZE * 4] = tmp3 + tmp4; ptr [DCTSIZE * 5] = tmp2 - tmp5; ptr [DCTSIZE * 6] = tmp1 - tmp6; ptr [DCTSIZE * 7] = tmp0 - tmp7; } ++freqs; ++ptr; ++quant; } ptr = tmp; for (i = 0; i < DCTSIZE; ++i) { double tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; double tmp10 = ptr [0] + ptr [4]; double tmp11 = ptr [0] - ptr [4]; double tmp13 = ptr [2] + ptr [6]; double tmp12 = (ptr [2] - ptr [6]) * 1.4142135623730950488 - tmp13; const double z13 = ptr [5] + ptr [3]; const double z10 = ptr [5] - ptr [3]; const double z11 = ptr [1] + ptr [7]; const double z12 = ptr [1] - ptr [7]; const double z5 = (z10 + z12) * 1.8477590650225735123; tmp0 = tmp10 + tmp13; tmp3 = tmp10 - tmp13; tmp1 = tmp11 + tmp12; tmp2 = tmp11 - tmp12; tmp7 = z11 + z13; tmp11 = (z11 - z13) * 1.4142135623730950488; tmp10 = 1.0823922002923939688 * z12 - z5; tmp12 = -2.6131259297527530557 * z10 + z5; tmp6 = tmp12 - tmp7; tmp5 = tmp11 - tmp6; tmp4 = tmp10 + tmp5; values [0] = (tmp0 + tmp7); values [1] = (tmp1 + tmp6); values [2] = (tmp2 + tmp5); values [3] = (tmp3 - tmp4); values [4] = (tmp3 + tmp4); values [5] = (tmp2 - tmp5); values [6] = (tmp1 - tmp6); values [7] = (tmp0 - tmp7); ptr += DCTSIZE; values += DCTSIZE; } } /* Convert floating-point value to JCOEF. */ static JCOEF coeff (const double x, const UINT16 quantval) { int result; /* Round up if the fractional part of the quantized values is equal to or greater than 0.72 (1 - 0.28). This empirically determined threshold gives the best results. */ if (x >= 0.0) result = (int)(x / (double)quantval + 0.28); else result = -(int)(-x / (double)quantval + 0.28); if (result < MIN_DCT_COEFF) return (JCOEF)MIN_DCT_COEFF; else if (result > MAX_DCT_COEFF) return (JCOEF)MAX_DCT_COEFF; else return (JCOEF)result; }