635 lines
19 KiB
C++
635 lines
19 KiB
C++
/*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 "precomp.hpp"
|
|
|
|
//*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. */
|