cvhmmobs.cpp.svn-base

来自「非结构化路识别」· SVN-BASE 代码 · 共 641 行 · 第 1/2 页

    /* 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 *) icvAlloc( 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;
                }
            }
        }
    }

    icvFree( &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,
                                           icvGetMatSize(mat), obs,
                                           dctSize, obsSize, delta ));
        break;
    case CV_32FC1:
        IPPI_CALL( icvImgToObs_DCT_32f_C1R( mat->data.fl, mat->step,
                                           icvGetMatSize(mat), obs,
                                           dctSize, obsSize, delta ));
        break;
    default:
        CV_ERROR( CV_StsUnsupportedFormat, "" );
    }

    __END__;
}


/* End of file. */