| /*M/////////////////////////////////////////////////////////////////////////////////////// |
| // |
| // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. |
| // |
| // By downloading, copying, installing or using the software you agree to this license. |
| // If you do not agree to this license, do not download, install, |
| // copy or use the software. |
| // |
| // |
| // Intel License Agreement |
| // For Open Source Computer Vision Library |
| // |
| // Copyright (C) 2000, Intel Corporation, all rights reserved. |
| // Third party copyrights are property of their respective owners. |
| // |
| // Redistribution and use in source and binary forms, with or without modification, |
| // are permitted provided that the following conditions are met: |
| // |
| // * Redistribution's of source code must retain the above copyright notice, |
| // this list of conditions and the following disclaimer. |
| // |
| // * Redistribution's in binary form must reproduce the above copyright notice, |
| // this list of conditions and the following disclaimer in the documentation |
| // and/or other materials provided with the distribution. |
| // |
| // * The name of Intel Corporation may not be used to endorse or promote products |
| // derived from this software without specific prior written permission. |
| // |
| // This software is provided by the copyright holders and contributors "as is" and |
| // any express or implied warranties, including, but not limited to, the implied |
| // warranties of merchantability and fitness for a particular purpose are disclaimed. |
| // In no event shall the Intel Corporation or contributors be liable for any direct, |
| // indirect, incidental, special, exemplary, or consequential damages |
| // (including, but not limited to, procurement of substitute goods or services; |
| // loss of use, data, or profits; or business interruption) however caused |
| // and on any theory of liability, whether in contract, strict liability, |
| // or tort (including negligence or otherwise) arising in any way out of |
| // the use of this software, even if advised of the possibility of such damage. |
| // |
| //M*/ |
| #include "_cvaux.h" |
| |
| //*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: icvImgToObs_DCT_8u32f_C1R |
| // Purpose: The function takes as input an image and returns the sequnce of observations |
| // to be used with an embedded HMM; Each observation is top-left block of DCT |
| // coefficient matrix. |
| // Context: |
| // Parameters: img - pointer to the original image ROI |
| // imgStep - full row width of the image in bytes |
| // roi - width and height of ROI in pixels |
| // obs - pointer to resultant observation vectors |
| // dctSize - size of the block for which DCT is calculated |
| // obsSize - size of top-left block of DCT coeffs matrix, which is treated |
| // as observation. Each observation vector consists of |
| // obsSize.width * obsSize.height floats. |
| // The following conditions should be satisfied: |
| // 0 < objSize.width <= dctSize.width, |
| // 0 < objSize.height <= dctSize.height. |
| // delta - dctBlocks are overlapped and this parameter specifies horizontal |
| // and vertical shift. |
| // Returns: |
| // CV_NO_ERR or error code |
| // Notes: |
| // The algorithm is following: |
| // 1. First, number of observation vectors per row and per column are calculated: |
| // |
| // Nx = floor((roi.width - dctSize.width + delta.width)/delta.width); |
| // Ny = floor((roi.height - dctSize.height + delta.height)/delta.height); |
| // |
| // So, total number of observation vectors is Nx*Ny, and total size of |
| // array obs must be >= Nx*Ny*obsSize.width*obsSize.height*sizeof(float). |
| // 2. Observation vectors are calculated in the following loop |
| // ( actual implementation may be different ), where |
| // I[x1:x2,y1:y2] means block of pixels from source image with |
| // x1 <= x < x2, y1 <= y < y2, |
| // D[x1:x2,y1:y2] means sub matrix of DCT matrix D. |
| // O[x,y] means observation vector that corresponds to position |
| // (x*delta.width,y*delta.height) in the source image |
| // ( all indices are counted from 0 ). |
| // |
| // for( y = 0; y < Ny; y++ ) |
| // { |
| // for( x = 0; x < Nx; x++ ) |
| // { |
| // D = DCT(I[x*delta.width : x*delta.width + dctSize.width, |
| // y*delta.height : y*delta.height + dctSize.height]); |
| // O[x,y] = D[0:obsSize.width, 0:obsSize.height]; |
| // } |
| // } |
| //F*/ |
| |
| /*comment out the following line to make DCT be calculated in floating-point arithmetics*/ |
| //#define _CV_INT_DCT |
| |
| /* for integer DCT only */ |
| #define DCT_SCALE 15 |
| |
| #ifdef _CV_INT_DCT |
| typedef int work_t; |
| |
| #define DESCALE CV_DESCALE |
| #define SCALE(x) CV_FLT_TO_FIX((x),DCT_SCALE) |
| #else |
| typedef float work_t; |
| |
| #define DESCALE(x,n) (float)(x) |
| #define SCALE(x) (float)(x) |
| #endif |
| |
| /* calculate dct transform matrix */ |
| static void icvCalcDCTMatrix( work_t * cfs, int n ); |
| |
| #define MAX_DCT_SIZE 32 |
| |
| static CvStatus CV_STDCALL |
| icvImgToObs_DCT_8u32f_C1R( uchar * img, int imgStep, CvSize roi, |
| float *obs, CvSize dctSize, |
| CvSize obsSize, CvSize delta ) |
| { |
| /* dct transform matrices: horizontal and vertical */ |
| work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; |
| work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; |
| |
| /* temporary buffers for dct */ |
| work_t temp0[MAX_DCT_SIZE * 4]; |
| work_t temp1[MAX_DCT_SIZE * 4]; |
| work_t *buffer = 0; |
| work_t *buf_limit; |
| |
| double s; |
| |
| int y; |
| int Nx, Ny; |
| |
| int n1 = dctSize.height, m1 = n1 / 2; |
| int n2 = dctSize.width, m2 = n2 / 2; |
| |
| if( !img || !obs ) |
| return CV_NULLPTR_ERR; |
| |
| if( roi.width <= 0 || roi.height <= 0 ) |
| return CV_BADSIZE_ERR; |
| |
| if( delta.width <= 0 || delta.height <= 0 ) |
| return CV_BADRANGE_ERR; |
| |
| if( obsSize.width <= 0 || dctSize.width < obsSize.width || |
| obsSize.height <= 0 || dctSize.height < obsSize.height ) |
| return CV_BADRANGE_ERR; |
| |
| if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE ) |
| return CV_BADRANGE_ERR; |
| |
| Nx = (roi.width - dctSize.width + delta.width) / delta.width; |
| Ny = (roi.height - dctSize.height + delta.height) / delta.height; |
| |
| if( Nx <= 0 || Ny <= 0 ) |
| return CV_BADRANGE_ERR; |
| |
| buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] )); |
| if( !buffer ) |
| return CV_OUTOFMEM_ERR; |
| |
| icvCalcDCTMatrix( tab_x, dctSize.width ); |
| icvCalcDCTMatrix( tab_y, dctSize.height ); |
| |
| buf_limit = buffer + obsSize.height * roi.width; |
| |
| for( y = 0; y < Ny; y++, img += delta.height * imgStep ) |
| { |
| int x, i, j, k; |
| work_t k0 = 0; |
| |
| /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */ |
| for( x = 0; x < roi.width; x++ ) |
| { |
| float is = 0; |
| work_t *buf = buffer + x; |
| work_t *tab = tab_y + 2; |
| |
| if( n1 & 1 ) |
| { |
| is = img[x + m1 * imgStep]; |
| k0 = ((work_t) is) * tab[-1]; |
| } |
| |
| /* first coefficient */ |
| for( j = 0; j < m1; j++ ) |
| { |
| float t0 = img[x + j * imgStep]; |
| float t1 = img[x + (n1 - 1 - j) * imgStep]; |
| float t2 = t0 + t1; |
| |
| t0 -= t1; |
| temp0[j] = (work_t) t2; |
| is += t2; |
| temp1[j] = (work_t) t0; |
| } |
| |
| buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT ); |
| if( (buf += roi.width) >= buf_limit ) |
| continue; |
| |
| /* other coefficients */ |
| for( ;; ) |
| { |
| s = 0; |
| |
| for( k = 0; k < m1; k++ ) |
| s += temp1[k] * tab[k]; |
| |
| buf[0] = DESCALE( s, PASS1_SHIFT ); |
| if( (buf += roi.width) >= buf_limit ) |
| break; |
| |
| tab += m1; |
| s = 0; |
| |
| if( n1 & 1 ) |
| { |
| k0 = -k0; |
| s = k0; |
| } |
| for( k = 0; k < m1; k++ ) |
| s += temp0[k] * tab[k]; |
| |
| buf[0] = DESCALE( s, PASS1_SHIFT ); |
| tab += m1; |
| |
| if( (buf += roi.width) >= buf_limit ) |
| break; |
| } |
| } |
| |
| k0 = 0; |
| |
| /* do transforms for rows. */ |
| for( x = 0; x + dctSize.width <= roi.width; x += delta.width ) |
| { |
| for( i = 0; i < obsSize.height; i++ ) |
| { |
| work_t *buf = buffer + x + roi.width * i; |
| work_t *tab = tab_x + 2; |
| float *obs_limit = obs + obsSize.width; |
| |
| s = 0; |
| |
| if( n2 & 1 ) |
| { |
| s = buf[m2]; |
| k0 = (work_t) (s * tab[-1]); |
| } |
| |
| /* first coefficient */ |
| for( j = 0; j < m2; j++ ) |
| { |
| work_t t0 = buf[j]; |
| work_t t1 = buf[n2 - 1 - j]; |
| work_t t2 = t0 + t1; |
| |
| t0 -= t1; |
| temp0[j] = (work_t) t2; |
| s += t2; |
| temp1[j] = (work_t) t0; |
| } |
| |
| *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT ); |
| |
| if( obs == obs_limit ) |
| continue; |
| |
| /* other coefficients */ |
| for( ;; ) |
| { |
| s = 0; |
| |
| for( k = 0; k < m2; k++ ) |
| s += temp1[k] * tab[k]; |
| |
| obs[0] = (float) DESCALE( s, PASS2_SHIFT ); |
| if( ++obs == obs_limit ) |
| break; |
| |
| tab += m2; |
| |
| s = 0; |
| |
| if( n2 & 1 ) |
| { |
| k0 = -k0; |
| s = k0; |
| } |
| for( k = 0; k < m2; k++ ) |
| s += temp0[k] * tab[k]; |
| obs[0] = (float) DESCALE( s, PASS2_SHIFT ); |
| |
| tab += m2; |
| if( ++obs == obs_limit ) |
| break; |
| } |
| } |
| } |
| } |
| |
| cvFree( &buffer ); |
| return CV_NO_ERR; |
| } |
| |
| |
| static CvStatus CV_STDCALL |
| icvImgToObs_DCT_32f_C1R( float * img, int imgStep, CvSize roi, |
| float *obs, CvSize dctSize, |
| CvSize obsSize, CvSize delta ) |
| { |
| /* dct transform matrices: horizontal and vertical */ |
| work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; |
| work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; |
| |
| /* temporary buffers for dct */ |
| work_t temp0[MAX_DCT_SIZE * 4]; |
| work_t temp1[MAX_DCT_SIZE * 4]; |
| work_t *buffer = 0; |
| work_t *buf_limit; |
| |
| double s; |
| |
| int y; |
| int Nx, Ny; |
| |
| int n1 = dctSize.height, m1 = n1 / 2; |
| int n2 = dctSize.width, m2 = n2 / 2; |
| |
| if( !img || !obs ) |
| return CV_NULLPTR_ERR; |
| |
| if( roi.width <= 0 || roi.height <= 0 ) |
| return CV_BADSIZE_ERR; |
| |
| if( delta.width <= 0 || delta.height <= 0 ) |
| return CV_BADRANGE_ERR; |
| |
| if( obsSize.width <= 0 || dctSize.width < obsSize.width || |
| obsSize.height <= 0 || dctSize.height < obsSize.height ) |
| return CV_BADRANGE_ERR; |
| |
| if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE ) |
| return CV_BADRANGE_ERR; |
| |
| Nx = (roi.width - dctSize.width + delta.width) / delta.width; |
| Ny = (roi.height - dctSize.height + delta.height) / delta.height; |
| |
| if( Nx <= 0 || Ny <= 0 ) |
| return CV_BADRANGE_ERR; |
| |
| buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] )); |
| if( !buffer ) |
| return CV_OUTOFMEM_ERR; |
| |
| icvCalcDCTMatrix( tab_x, dctSize.width ); |
| icvCalcDCTMatrix( tab_y, dctSize.height ); |
| |
| buf_limit = buffer + obsSize.height * roi.width; |
| |
| imgStep /= sizeof(img[0]); |
| |
| for( y = 0; y < Ny; y++, img += delta.height * imgStep ) |
| { |
| int x, i, j, k; |
| work_t k0 = 0; |
| |
| /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */ |
| for( x = 0; x < roi.width; x++ ) |
| { |
| float is = 0; |
| work_t *buf = buffer + x; |
| work_t *tab = tab_y + 2; |
| |
| if( n1 & 1 ) |
| { |
| is = img[x + m1 * imgStep]; |
| k0 = ((work_t) is) * tab[-1]; |
| } |
| |
| /* first coefficient */ |
| for( j = 0; j < m1; j++ ) |
| { |
| float t0 = img[x + j * imgStep]; |
| float t1 = img[x + (n1 - 1 - j) * imgStep]; |
| float t2 = t0 + t1; |
| |
| t0 -= t1; |
| temp0[j] = (work_t) t2; |
| is += t2; |
| temp1[j] = (work_t) t0; |
| } |
| |
| buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT ); |
| if( (buf += roi.width) >= buf_limit ) |
| continue; |
| |
| /* other coefficients */ |
| for( ;; ) |
| { |
| s = 0; |
| |
| for( k = 0; k < m1; k++ ) |
| s += temp1[k] * tab[k]; |
| |
| buf[0] = DESCALE( s, PASS1_SHIFT ); |
| if( (buf += roi.width) >= buf_limit ) |
| break; |
| |
| tab += m1; |
| s = 0; |
| |
| if( n1 & 1 ) |
| { |
| k0 = -k0; |
| s = k0; |
| } |
| for( k = 0; k < m1; k++ ) |
| s += temp0[k] * tab[k]; |
| |
| buf[0] = DESCALE( s, PASS1_SHIFT ); |
| tab += m1; |
| |
| if( (buf += roi.width) >= buf_limit ) |
| break; |
| } |
| } |
| |
| k0 = 0; |
| |
| /* do transforms for rows. */ |
| for( x = 0; x + dctSize.width <= roi.width; x += delta.width ) |
| { |
| for( i = 0; i < obsSize.height; i++ ) |
| { |
| work_t *buf = buffer + x + roi.width * i; |
| work_t *tab = tab_x + 2; |
| float *obs_limit = obs + obsSize.width; |
| |
| s = 0; |
| |
| if( n2 & 1 ) |
| { |
| s = buf[m2]; |
| k0 = (work_t) (s * tab[-1]); |
| } |
| |
| /* first coefficient */ |
| for( j = 0; j < m2; j++ ) |
| { |
| work_t t0 = buf[j]; |
| work_t t1 = buf[n2 - 1 - j]; |
| work_t t2 = t0 + t1; |
| |
| t0 -= t1; |
| temp0[j] = (work_t) t2; |
| s += t2; |
| temp1[j] = (work_t) t0; |
| } |
| |
| *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT ); |
| |
| if( obs == obs_limit ) |
| continue; |
| |
| /* other coefficients */ |
| for( ;; ) |
| { |
| s = 0; |
| |
| for( k = 0; k < m2; k++ ) |
| s += temp1[k] * tab[k]; |
| |
| obs[0] = (float) DESCALE( s, PASS2_SHIFT ); |
| if( ++obs == obs_limit ) |
| break; |
| |
| tab += m2; |
| |
| s = 0; |
| |
| if( n2 & 1 ) |
| { |
| k0 = -k0; |
| s = k0; |
| } |
| for( k = 0; k < m2; k++ ) |
| s += temp0[k] * tab[k]; |
| obs[0] = (float) DESCALE( s, PASS2_SHIFT ); |
| |
| tab += m2; |
| if( ++obs == obs_limit ) |
| break; |
| } |
| } |
| } |
| } |
| |
| cvFree( &buffer ); |
| return CV_NO_ERR; |
| } |
| |
| |
| static void |
| icvCalcDCTMatrix( work_t * cfs, int n ) |
| { |
| static const double sqrt2 = 1.4142135623730950488016887242097; |
| static const double pi = 3.1415926535897932384626433832795; |
| |
| static const double sincos[16 * 2] = { |
| 1.00000000000000000, 0.00000000000000006, |
| 0.70710678118654746, 0.70710678118654757, |
| 0.49999999999999994, 0.86602540378443871, |
| 0.38268343236508978, 0.92387953251128674, |
| 0.30901699437494740, 0.95105651629515353, |
| 0.25881904510252074, 0.96592582628906831, |
| 0.22252093395631439, 0.97492791218182362, |
| 0.19509032201612825, 0.98078528040323043, |
| 0.17364817766693033, 0.98480775301220802, |
| 0.15643446504023087, 0.98768834059513777, |
| 0.14231483827328514, 0.98982144188093268, |
| 0.13052619222005157, 0.99144486137381038, |
| 0.12053668025532305, 0.99270887409805397, |
| 0.11196447610330786, 0.99371220989324260, |
| 0.10452846326765346, 0.99452189536827329, |
| 0.09801714032956060, 0.99518472667219693, |
| }; |
| |
| #define ROTATE( c, s, dc, ds ) \ |
| { \ |
| t = c*dc - s*ds; \ |
| s = c*ds + s*dc; \ |
| c = t; \ |
| } |
| |
| #define WRITE2( j, a, b ) \ |
| { \ |
| cfs[j] = SCALE(a); \ |
| cfs2[j] = SCALE(b); \ |
| } |
| |
| double t, scale = 1. / sqrt( (double)n ); |
| int i, j, m = n / 2; |
| |
| cfs[0] = SCALE( scale ); |
| scale *= sqrt2; |
| cfs[1] = SCALE( scale ); |
| cfs += 2 - m; |
| |
| if( n > 1 ) |
| { |
| double a0, b0; |
| double da0, db0; |
| work_t *cfs2 = cfs + m * n; |
| |
| if( n <= 16 ) |
| { |
| da0 = a0 = sincos[2 * n - 1]; |
| db0 = b0 = sincos[2 * n - 2]; |
| } |
| else |
| { |
| t = pi / (2 * n); |
| da0 = a0 = cos( t ); |
| db0 = b0 = sin( t ); |
| } |
| |
| /* other rows */ |
| for( i = 1; i <= m; i++ ) |
| { |
| double a = a0 * scale; |
| double b = b0 * scale; |
| double da = a0 * a0 - b0 * b0; |
| double db = a0 * b0 + a0 * b0; |
| |
| cfs += m; |
| cfs2 -= m; |
| |
| for( j = 0; j < m; j += 2 ) |
| { |
| WRITE2( j, a, b ); |
| ROTATE( a, b, da, db ); |
| if( j + 1 < m ) |
| { |
| WRITE2( j + 1, a, -b ); |
| ROTATE( a, b, da, db ); |
| } |
| } |
| |
| ROTATE( a0, b0, da0, db0 ); |
| } |
| } |
| #undef ROTATE |
| #undef WRITE2 |
| } |
| |
| |
| CV_IMPL void |
| cvImgToObs_DCT( const void* arr, float *obs, CvSize dctSize, |
| CvSize obsSize, CvSize delta ) |
| { |
| CV_FUNCNAME( "cvImgToObs_DCT" ); |
| |
| __BEGIN__; |
| |
| CvMat stub, *mat = (CvMat*)arr; |
| |
| CV_CALL( mat = cvGetMat( arr, &stub )); |
| |
| switch( CV_MAT_TYPE( mat->type )) |
| { |
| case CV_8UC1: |
| IPPI_CALL( icvImgToObs_DCT_8u32f_C1R( mat->data.ptr, mat->step, |
| cvGetMatSize(mat), obs, |
| dctSize, obsSize, delta )); |
| break; |
| case CV_32FC1: |
| IPPI_CALL( icvImgToObs_DCT_32f_C1R( mat->data.fl, mat->step, |
| cvGetMatSize(mat), obs, |
| dctSize, obsSize, delta )); |
| break; |
| default: |
| CV_ERROR( CV_StsUnsupportedFormat, "" ); |
| } |
| |
| __END__; |
| } |
| |
| |
| /* End of file. */ |
