uiplutils.cpp

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

                    dstData[y*dstStep + x*ach + c] = (uchar)val;
                }
        break;

    case IPL_DEPTH_32F:
        {
        int *fltSrc = (int*)srcData;
        int *fltDst = (int*)dstData;

        srcStep /= sizeof(int);
        dstStep /= sizeof(int);

        for( y = 0; y < sz.height; y++ )
            for( x = 0; x < sz.width; x++ )
                for( c = 0; c < ach; c++ )
                {
                    int  val = int_extr0;
                    int  x1, y1 = y - anchorY;
                    int  i, j;

                    for( i = 0; i < rows; i++, y1++ )
                        for( j = 0, x1 = x - anchorX; j < cols; j++, x1++ )
                        {
                            int x2 = x1, y2 = y1;
                            if( x1 < 0 ) x2 = 0;
                            if( x1 >= sz.width ) x2 = sz.width - 1;
                            if( y1 < 0 ) y2 = 0;
                            if( y1 >= sz.height ) y2 = sz.height - 1;

                            /* unsigned comparsion helps to check for positivity */
                            if( mask[i*cols + j] != 0 )
                            {
                                int temp = fltSrc[y2*srcStep + x2*ach + c];
                                temp = ATS_TOGGLE_FLT( temp );
                                val = ATS_MIN( val, temp );
                            }
                        }

                    fltDst[y*dstStep + x*ach + c] = ATS_TOGGLE_FLT(val);
                }
        }
        break;
    default: assert(0);
    }
}


/*
   The function applies min filter using specified structuring element.
   float numbers are processed using integer arithmetics (IEEE 754 format is assumed)
*/
void atsMaxFilterEx( IplImage* src, IplImage* dst,
                     IplConvKernel* B )
{
    const    int  int_extr0   = 0x80000000;

    uchar*   srcData;
    uchar*   dstData;

    int      srcStep, dstStep;
    int      depth;
    int      ach; /* how many non-alpha channels and all the channels */
    CvSize  sz;

    int      cols    = B->nCols;
    int      rows    = B->nRows;
    int      anchorX = B->anchorX;
    int      anchorY = B->anchorY;
    int*     mask = B->values;

    int      x, y, c;

    atsGetImageInfo( src, (void**)&srcData, &srcStep, &sz, &depth, &ach, 0 );
    atsGetImageInfo( dst, (void**)&dstData, &dstStep, 0, 0, 0, 0 );

    switch( depth )
    {
    case IPL_DEPTH_8U:
        for( y = 0; y < sz.height; y++ )
            for( x = 0; x < sz.width; x++ )
                for( c = 0; c < ach; c++ )
                {
                    int val = int_extr0;
                    int  x1, y1 = y - anchorY;
                    int  i, j;

                    for( i = 0; i < rows; i++, y1++ )
                        for( j = 0, x1 = x - anchorX; j < cols; j++, x1++ )
                        {
                            int x2 = x1, y2 = y1;
                            if( x1 < 0 ) x2 = 0;
                            if( x1 >= sz.width ) x2 = sz.width - 1;
                            if( y1 < 0 ) y2 = 0;
                            if( y1 >= sz.height ) y2 = sz.height - 1;

                            /* unsigned comparsion helps to check for positivity */
                            if( mask[i*cols + j] != 0 )
                            {
                                int temp = srcData[y2*srcStep + x2*ach + c];
                                val = ATS_MAX( val, temp );
                            }
                        }

                    dstData[y*dstStep + x*ach + c] = (uchar)val;
                }
        break;

    case IPL_DEPTH_32F:
        {
        int *fltSrc = (int*)srcData;
        int *fltDst = (int*)dstData;

        srcStep /= sizeof(int);
        dstStep /= sizeof(int);

        for( y = 0; y < sz.height; y++ )
            for( x = 0; x < sz.width; x++ )
                for( c = 0; c < ach; c++ )
                {
                    int  val = int_extr0;
                    int  x1, y1 = y - anchorY;
                    int  i, j;

                    for( i = 0; i < rows; i++, y1++ )
                        for( j = 0, x1 = x - anchorX; j < cols; j++, x1++ )
                        {
                            int x2 = x1, y2 = y1;
                            if( x1 < 0 ) x2 = 0;
                            if( x1 >= sz.width ) x2 = sz.width - 1;
                            if( y1 < 0 ) y2 = 0;
                            if( y1 >= sz.height ) y2 = sz.height - 1;

                            /* unsigned comparsion helps to check for positivity */
                            if( mask[i*cols + j] != 0 )
                            {
                                int temp = fltSrc[y2*srcStep + x2*ach + c];
                                temp = ATS_TOGGLE_FLT( temp );
                                val = ATS_MAX( val, temp );
                            }
                        }

                    fltDst[y*dstStep + x*ach + c] = ATS_TOGGLE_FLT(val);
                }
        }
        break;
    default: assert(0);
    }
}



/*
   Replicates left and right ROI borders dx times,
   top and bottom ROI borders dy times.
*/
void atsReplicateBorders( IplImage* img, int dx, int dy )
{
    int      step;
    uchar*   data;
    int      bt_pix;
    int      x, y;
    int      x0, x1, y0, y1;
    int      count;
    int      copy_width;
    CvSize  sz;

    iplSetBorderMode( img, IPL_BORDER_REPLICATE, IPL_SIDE_ALL, 0 );

    if( !img->roi ) return;

    atsGetImageInfo( img, (void**)&data, &step, &sz, 0, 0, &bt_pix );

    x0 = img->roi->xOffset;
    y0 = img->roi->yOffset;

    x1 = x0 + sz.width - 1;
    y1 = y0 + sz.height - 1;

    copy_width = sz.width * bt_pix;

    /* left border */
    count = (x0 - dx >= 0 ? dx : x0)*bt_pix;
    copy_width += count;

    if( count > 0 )
    {
        for( y = 0; y < sz.height; y++, data += step )
            for( x = count; x > 0; x -= bt_pix ) memcpy( data - x, data, bt_pix );
        data -= sz.height*step;
    }

    /* right border */
    count = (x1 + dx < img->width ? dx : img->width - x1 - 1)*bt_pix;
    copy_width += count;

    if( count > 0 )
    {
        data += (sz.width - 1)*bt_pix;
        for( y = 0; y < sz.height; y++, data += step )
            for( x = count; x > 0; x -= bt_pix ) memcpy( data + x, data, bt_pix );
        data -= sz.height*step + (sz.width - 1)*bt_pix;
    }

    /* top border */
    count = (y0 - dy >= 0 ? dy : y0)*step;
    if( count > 0 )
    {
        for( y = count; y > 0; y -= step ) memcpy( data - y, data, copy_width );
    }

    /* bottom border */
    count = (y1 + dy < img->height ? dy : img->height - y1 - 1)*step;
    if( count > 0 )
    {
        data += (sz.height - 1)*step;
        for( y = count; y > 0; y -= step ) memcpy( data + y, data, copy_width );
    }
}


/*
   The convolution function.
   Supports only 32fC1 images
*/
void atsConvolve( IplImage* src, IplImage* dst,
                  IplConvKernelFP* ker )
{
    float*   srcData;
    float*   dstData;

    int      srcStep, dstStep;
    CvSize  sz;

    int      cols     = ker->nCols;
    int      rows     = ker->nRows;
    int      anchorX  = ker->anchorX;
    int      anchorY  = ker->anchorY;
    float*   ker_data = ker->values;
    int      x, y;

    atsGetImageInfo( src, (void**)&srcData, &srcStep, &sz, 0, 0, 0 );
    atsGetImageInfo( dst, (void**)&dstData, &dstStep, 0, 0, 0, 0 );

    srcStep /= sizeof(int);
    dstStep /= sizeof(int);

    for( y = 0; y < sz.height; y++ )
        for( x = 0; x < sz.width; x++ )
        {
            double sum = 0;
            int   x1, y1 = y - anchorY;
            int  i, j;

            for( i = 0; i < rows; i++, y1++ )
                for( j = 0, x1 = x - anchorX; j < cols; j++, x1++ )
                {
                    /* unsigned comparsion helps to check for positivity */
                    if( (unsigned)x1 < (unsigned)sz.width  &&
                        (unsigned)y1 < (unsigned)sz.height )
                    {
                        sum += srcData[y1*srcStep + x1] * ker_data[i*rows + j];
                    }
                    else
                    {
                        int x2 = x1, y2 = y1;
                        if( x2 < 0 ) x2 = 0;
                        if( x2 >= sz.width ) x2 = sz.width - 1;
                        if( y2 < 0 ) y2 = 0;
                        if( y2 >= sz.height ) y2 = sz.height - 1;

                        sum += srcData[y2*srcStep + x2] * ker_data[i*rows + j];
                    }
                }

            dstData[y*dstStep + x] = (float)sum;
        }
}


/*
 Function caluclates miscellaneous statistics for the image:
 min & max with locations, number of non-zero pixels, sum of all pixels,
 mean value, standard deviation, c_norm( maximum of absolute values ),
 l1_norm ( sum of absolute values ), l2_norm (square root of sum of values squares)
*/
void  atsCalcImageStatistics(
        IplImage* img, IplImage* mask,
        double* _min_val, double* _max_val,
        CvPoint* _min_loc, CvPoint* _max_loc,
        int* _non_zero, double* _sum,
        double* _mean, double* _svd,
        double* _c_norm, double* _l1_norm, double* _l2_norm,
        int* _mask_pix )
{
    float*   img_data = 0;
    uchar*   mask_data = 0;
    int      img_step = 0, mask_step = 0;

    double   sum = 0, /* pixels sum */
             l1_norm = 0, /* sum of absolute values */
             l2_norm = 0; /* sum of squares (used in stddev and l2 norm) */
    int      non_zero = 0; /* number of non-zero pixels */
    int      mask_counter = 0; /* number of 1-pixels in mask */
    double   c_norm = 0, /* maximum of absolute value */
             min_val = DBL_MAX, /* minimum value */
             max_val = -DBL_MAX; /* maximum value */
    CvPoint min_loc = {0,0}, max_loc = {0,0}; /* current positions of minimum and maximum */
    double   inv_mask_counter = 0;
    int      i, j;
    CvSize  sz;
    int      coi = img->roi ? img->roi->coi : 0;
    int      depth = 0, ch = 0;

    atsGetImageInfo( img, (void**)&img_data, &img_step, &sz, &depth, &ch, 0 );
    assert( depth == IPL_DEPTH_32F && (ch == 1 || (ch == 3 && coi != 0)));

    if( mask )
    {
        atsGetImageInfo( mask, (void**)&mask_data, &mask_step, 0, 0, 0, 0 );
    }

    img_step /= 4;
    img_data += coi == 0 ? 0 : coi - 1;

    /* iterate through the image */
    for( i = 0; i < sz.height; i++, img_data += img_step, mask_data += mask_step )
    {
        for( j = 0; j < sz.width; j++ )
        {
            /* check mask pixel. if mask is absent - all pixels considered */
            if( !mask_data || mask_data[j] != 0 )
            {
                double val = img_data[j*ch];

                /* update minimum */
                if( val < min_val )
                {
                    min_val = val;
                    min_loc.x = j;
                    min_loc.y = i;
                }

                /* update maximum */
                if( val > max_val )
                {
                    max_val = val;
                    max_loc.x = j;
                    max_loc.y = i;
                }

                /* update non-zero */
                non_zero += val != 0;

                /* update sum */
                sum += val;
                val = fabs(val);

                /* update norms */
                if( val > c_norm ) c_norm = val;
                l1_norm += val;
                l2_norm += val*val;

                /* increase mask pixels counter */
                mask_counter++;
            }
        }
    }

    if( mask_counter != 0 )
    {
        inv_mask_counter = 1./mask_counter;
        if( _min_val ) *_min_val = min_val;
        if( _max_val ) *_max_val = max_val;
        if( _min_loc ) *_min_loc = min_loc;
        if( _max_loc ) *_max_loc = max_loc;
    }
    else
    {
        /* set to default if no pixel processed */
        min_loc.x = min_loc.y = 0;
        max_loc.x = max_loc.y = 0;
        min_val = max_val = 0.f;
    }

    if( _mask_pix ) *_mask_pix = mask_counter;
    if( _non_zero ) *_non_zero = non_zero;
    if( _c_norm  ) *_c_norm  = c_norm;
    if( _l1_norm ) *_l1_norm = l1_norm;
    if( _l2_norm ) *_l2_norm = sqrt(l2_norm);
    if( _sum )     *_sum = sum;
    sum *= inv_mask_counter;
    if( _mean )    *_mean = sum;
    if( _svd )     *_svd = sqrt(l2_norm*inv_mask_counter - sum*sum);
}


/*
  Function calculates spatial and central moments up to third order.
  <binary> mode means that pixels values treated as 1 if they are non zero and 0 if zero.
*/
void    atsCalcMoments( IplImage* img, AtsMomentState* state, int binary )
{
    int      x, y;
    uchar*   img_data;
    int      img_step;
    CvSize  sz;
    int      coi = img->roi ? img->roi->coi : 0;
    int      depth = 0, ch = 0;
    double   cx = 0, cy = 0;
    AtsMomentState s;

    memset( &s, 0, sizeof(s));

    atsGetImageInfo( img, (void**)&img_data, &img_step, &sz, &depth, &ch, 0 );
    assert( depth == IPL_DEPTH_32F || depth == IPL_DEPTH_8U || depth == IPL_DEPTH_8S );
    assert( ch == 1 || (ch == 3 && coi != 0));
    assert( sz.width <= 2048 && sz.height <= 2048 );

    img_data += (coi == 0 ? 0 : coi - 1)*((depth & 255) >> 3);

    if( depth == IPL_DEPTH_8U || depth == IPL_DEPTH_8S ||
        (depth == IPL_DEPTH_32F && binary)) /* integer-value spatial moments */
    {
        int64  m00 = 0,
                 m10 = 0, m01 = 0,
                 m20 = 0, m11 = 0, m02 = 0,
                 m30 = 0, m21 = 0, m12 = 0, m03 = 0;

        int x2, y2; /* x^2 & y^2 */
        uchar* data = (uchar*)img_data;

        if( binary != 0 && depth == IPL_DEPTH_8S ) depth = IPL_DEPTH_8U;

        /* calc spatial moments */
        for( y = 0, y2 = 0; y < sz.height; y2 += 2*y + 1, y++, data += img_step )
        {
            int tm0 = 0, tm1 = 0;
            int64 tm2 = 0, tm3 = 0;

            for( x = 0, x2 = 0; x < sz.width; x2 += 2*x + 1, x++ )
            {
                int p, xp;

                if( depth == IPL_DEPTH_8U )
                {
                    p = data[x*ch];
                    if( binary ) p = p != 0;
                }
                else if( depth == IPL_DEPTH_8S )