amathutils.cpp

微软的基于HMM的人脸识别原代码, 非常经典的说

                {
                    code = -1;
                    goto test_exit;
                }
            }

            err = fabs( fm - radius )/(radius + 1e-15);
            if( err > max_mag_err )
            {
                max_mag_err = err;
                mm_err_x = x;
                mm_err_y = y;
                if( max_mag_err > success_mag_error_level )
                {
                    code = -2;
                    goto test_exit;
                }
            }
            */
            data_y[k] = (float)y;
            data_x[k] = (float)x;
            std_phase[k] = (float)angle1;
            std_mag[k] = (float)radius;
        }
    }

    /* vector function test */

    /* shuffle data */
    atsRandInit( &state, 0, total, seed );

    for( i = 0; i < total/2; i++ )
    {
        int i0 = atsRand32s( &state );
        int i1 = atsRand32s( &state );

        float temp;

        ATS_SWAP( data_y[i0], data_y[i1], temp );
        ATS_SWAP( data_x[i0], data_x[i1], temp );
        ATS_SWAP( std_phase[i0], std_phase[i1], temp );
        ATS_SWAP( std_mag[i0], std_mag[i1], temp );
    }

    cvbCartToPolar( data_x, data_y, fast_mag, fast_phase, total );

    for( i = 0; i < total; i++ )
    {
        double err = atsCompareAngles( std_phase[i], fast_phase[i] );
        if( err > max_phase_berr )
        {
            max_phase_berr = err;
            mp_err_by = data_y[i];
            mp_err_bx = data_x[i];

            if( max_phase_berr > success_phase_error_level )
            {
                code = -4;
                goto test_exit;
            }
        }
        err = fabs( fast_mag[i] - std_mag[i] )/(std_mag[i] + 1e-15);
        if( err > max_mag_berr )
        {
            max_mag_berr = err;
            mm_err_by = data_y[i];
            mm_err_bx = data_x[i];
            if( max_mag_berr > success_mag_error_level )
            {
                code = -5;
                goto test_exit;
            }
        }
    }

test_exit:

    icvFree( &data_y );
    icvFree( &data_x );
    icvFree( &std_mag );
    icvFree( &fast_mag );
    icvFree( &std_phase );
    icvFree( &fast_phase );

    if( code != -3 )
    {
        trsWrite( ATS_LST, "phase:\nmax scalar err: %g at y %g, x %g\n"
                           "max vector err: %g at y %g, x %g, seed %08x\n",
                           max_phase_err, mp_err_y, mp_err_x,
                           max_phase_berr, mp_err_by, mp_err_bx, seed );

        trsWrite( ATS_LST, "magnitude:\nmax scalar err: %g at y %g, x %g\n"
                           "max vector err: %g at y %g, x %g, seed %08x",
                           max_mag_err, mm_err_y, mm_err_x,
                           max_mag_berr, mm_err_by, mm_err_bx, seed );

        return max_phase_err <= success_phase_error_level &&
               max_mag_err <= success_mag_error_level &&
               max_phase_berr <= success_phase_error_level &&
               max_mag_berr <= success_mag_error_level ?
            trsResult( TRS_OK, "No errors" ) :
               max_mag_err > success_mag_error_level ||
               max_phase_err > success_phase_error_level ?
            trsResult( TRS_FAIL, "Bad scalar accuracy" ) :
            trsResult( TRS_FAIL, "Bad vector accuracy" ) ;
    }
    else
    {
        return trsResult( TRS_FAIL, "Vector function returns error code" );
    }
}


static int test_recip( void )
{
    double success_error_level = 1e-6;
    const int mbits = 23;
    const int exp_max = 255;

    /* position where the maximum error occured */
    double merr_bx = 0;

    unsigned seed = atsGetSeed();

    /* test parameters */
    int     exp = 0, mant = 0;
    double  max_berr = 0.; /* maximum vector error */
    int     code = 0;
    int     k;
    int     i = 0, m = 0, e = 0, total = 0;
    float*  data_x   = 0;
    float*  std      = 0;
    float*  result   = 0;
    AtsRandState state;

    read_sqrt_params();

    for( mant = 0; mant < (1<<mbits); mant += mant_step, m++ );
    for( exp = 1; exp < exp_max; exp += exp_step, e++ );
    total = m*e;

    data_x  = (float*)icvAlloc( total*sizeof(data_x[0]) );
    std     = (float*)icvAlloc( total*sizeof(std[0]) );
    result  = (float*)icvAlloc( total*sizeof(result[0]) );

    /* scalar function test */
    for( exp = 1, k = 0; exp < exp_max; exp += exp_step )
    {
        for( mant = 0; mant < (1<<mbits); mant += mant_step, k++ )
        {
            float x;

            *((unsigned*)&x) = (exp << mbits) | mant;

            data_x[k] = x;
            std[k] = 1.f/x;
        }
    }

    /* vector function test */

    /* shuffle data */
    atsRandInit( &state, 0, total, seed );

    for( i = 0; i < total/2; i++ )
    {
        int i0 = atsRand32s( &state );
        int i1 = atsRand32s( &state );

        float temp;

        ATS_SWAP( data_x[i0], data_x[i1], temp );
        ATS_SWAP( std[i0], std[i1], temp );
    }

    cvbReciprocal( data_x, result, total );

    for( i = 0; i < total; i++ )
    {
        double err = fabs(result[i] - std[i])/(std[i] + DBL_EPSILON);
        if( err > max_berr )
        {
            max_berr = err;
            merr_bx = data_x[i];

            if( max_berr > success_error_level )
            {
                code = -3;
                goto test_exit;
            }
        }
    }

test_exit:

    icvFree( &data_x );
    icvFree( &std );
    icvFree( &result );

    if( code != -2 )
    {
        trsWrite( ATS_LST, "max vector err: %g at x = %g, seed %08x",
                           max_berr, merr_bx, seed );

        return max_berr <= success_error_level ?
            trsResult( TRS_OK, "No errors" ) :
            trsResult( TRS_FAIL, "Bad accuracy" ) ;
    }
    else
    {
        return trsResult( TRS_FAIL, "Vector function returns error code" );
    }
}


static int test_log( void )
{
    const double success_error_level = 1e-6;
    const int mbits = 52;
    const int exp_max = 2047;

    /* position where the maximum error occured */
    double  merr_bx = 0;

    unsigned seed = atsGetSeed();

    /* test parameters */
    int64   mant_step;
    int     exp_step;
    int     exp = 0;
    int64   mant = 0;
    double  max_berr = 0.; /* maximum vector error */
    int     code = 0;
    int     k;
    int     i = 0, m = 0, e = 0, total = 0;
    double* data_x  = 0;
    float*  std     = 0;
    float*  result  = 0;
    AtsRandState state;

    read_log_params();

    mant_step = log_mant_step;
    exp_step = log_exp_step;

    for( mant = 0; mant < ((int64)1<<mbits); mant += mant_step, m++ );
    for( exp = 1; exp < exp_max; exp += exp_step, e++ );
    total = m*e;

    data_x  = (double*)icvAlloc( total*sizeof(data_x[0]) );
    std     = (float*)icvAlloc( total*sizeof(std[0]) );
    result  = (float*)icvAlloc( total*sizeof(result[0]) );

    /* scalar function test */
    for( exp = 1, k = 0; exp < exp_max; exp += exp_step )
    {
        for( mant = 0; mant < ((int64)1<<mbits); mant += mant_step, k++ )
        {
            double x;

            *((int64*)&x) = ((int64)exp << mbits) | mant;

            data_x[k] = x;
            std[k] = (float)log(x);
        }
    }

    /* vector function test */

    /* shuffle data */
    atsRandInit( &state, 0, total, seed );

    for( i = 0; i < total/2; i++ )
    {
        int i0 = atsRand32s( &state );
        int i1 = atsRand32s( &state );

        float temp;
        double dtemp;

        ATS_SWAP( data_x[i0], data_x[i1], dtemp );
        ATS_SWAP( std[i0], std[i1], temp );
    }

    cvbFastLog( data_x, result, total );

    for( i = 0; i < total; i++ )
    {
        double err = fabs(result[i] - std[i])/(std[i] + DBL_EPSILON);
        if( err > max_berr )
        {
            max_berr = err;
            merr_bx = data_x[i];

            if( max_berr > success_error_level )
            {
                code = -3;
                goto test_exit;
            }
        }
    }

test_exit:

    icvFree( &data_x );
    icvFree( &std );
    icvFree( &result );

    if( code != -2 )
    {
        trsWrite( ATS_LST, "max vector err: %g at x = %g, seed %08x",
                           max_berr, merr_bx, seed );

        return max_berr <= success_error_level ?
            trsResult( TRS_OK, "No errors" ) :
            trsResult( TRS_FAIL, "Bad accuracy" ) ;
    }
    else
    {
        return trsResult( TRS_FAIL, "Vector function returns error code" );
    }
}



static int test_exp( void )
{
    const double success_error_level = 1e-5;
    const int mbits = 23;

    /* position where the maximum error occured */
    double  merr_bx = 0;

    unsigned seed = atsGetSeed();

    /* test parameters */
    int     mant_step = 0;
    int     exp_max = 0;
    int     _exp = 0;
    int     mant = 0;
    double  max_berr = 0.; /* maximum vector error */
    int     code = 0;
    int     k;
    int     i = 0, m = 0, e = 0, total = 0;
    float*  data_x  = 0;
    double* std     = 0;
    double* result  = 0;
    AtsRandState state;

    read_exp_params();

    mant_step = exp_mant_step;
    exp_max = exp_exp_max;

    for( mant = 0; mant < ((int64)1<<mbits); mant += mant_step, m++ );
    for( _exp = 127 - exp_max; _exp <= 127 + exp_max; _exp++, e++ );
    total = m*e;

    data_x  = (float*)icvAlloc( total*sizeof(data_x[0]) );
    std     = (double*)icvAlloc( total*sizeof(std[0]) );
    result  = (double*)icvAlloc( total*sizeof(result[0]) );

    /* scalar function test */
    for( _exp = 127 - exp_max, k = 0; _exp <= 127 + exp_max; _exp++ )
    {
        for( mant = 0; mant < ((int64)1<<mbits); mant += mant_step, k++ )
        {
            float x;

            *((int*)&x) = (_exp << mbits) | mant;

            data_x[k] = x;
            std[k] = exp(x);
        }
    }

    /* vector function test */

    /* shuffle data */
    atsRandInit( &state, 0, total, seed );

    for( i = 0; i < total/2; i++ )
    {
        int i0 = atsRand32s( &state );
        int i1 = atsRand32s( &state );

        float temp;
        double dtemp;

        ATS_SWAP( data_x[i0], data_x[i1], temp );
        ATS_SWAP( std[i0], std[i1], dtemp );
    }

    cvbFastExp( data_x, result, total );

    for( i = 0; i < total; i++ )
    {
        double err = fabs(result[i] - std[i])/(std[i] + DBL_EPSILON);
        if( err > max_berr )
        {
            max_berr = err;
            merr_bx = data_x[i];

            if( max_berr > success_error_level )
            {
                code = -3;
                goto test_exit;
            }
        }
    }

test_exit:

    icvFree( &data_x );
    icvFree( &std );
    icvFree( &result );

    if( code != -2 )
    {
        trsWrite( ATS_LST, "max vector err: %g at x = %g, seed %08x",
                           max_berr, merr_bx, seed );

        return max_berr <= success_error_level ?
            trsResult( TRS_OK, "No errors" ) :
            trsResult( TRS_FAIL, "Bad accuracy" ) ;
    }
    else
    {
        return trsResult( TRS_FAIL, "Vector function returns error code" );
    }
}


#define SQRT_32F      0
#define INV_SQRT_32F  1

void InitAMathUtils( void )
{
    /* Registering test functions */
    trsReg( func_names[0], test_desc, atsAlgoClass, test_arctan );
    trsRegArg( func_names[1], test_desc, atsAlgoClass, test_sqrt, SQRT_32F );
    trsRegArg( func_names[2], test_desc, atsAlgoClass, test_sqrt, INV_SQRT_32F );
    trsReg( func_names[3], test_desc, atsAlgoClass, test_cart_to_polar );
    trsReg( func_names[4], test_desc, atsAlgoClass, test_recip );
    trsReg( func_names[5], test_desc, atsAlgoClass, test_log );
    trsReg( func_names[6], test_desc, atsAlgoClass, test_exp );

} /* InitAMathUtils */

/* End of file. */