C++ PicArray类代码示例

void dirac::CWMFilterComponent( PicArray& pic_data, const int strength )
{
    // Do centre-weighted median denoising

    PicArray pic_copy( pic_data );

    const int width( 3 );
    const int offset( (width-1)/2 );
    const int centre_weight = std::max(1, (width*width+1)-strength );
    const int list_length = centre_weight+(width*width)-1;

    ValueType* val_list = new ValueType[list_length];

    for (int j=offset; j<pic_data.LengthY()-offset; ++j){
        for (int i=offset; i<pic_data.LastX()-offset; ++i){
            // Make the value list
            int pos=0;
            for (; pos<centre_weight-1; ++pos)
                val_list[pos] = pic_copy[j][i];

            for (int s=-offset; s<=offset; ++s){
                for (int r=-offset; r<=offset; ++r){
                    val_list[pos]=pic_copy[j+s][i+r];
                    pos++;
                }// r
            }// s

            pic_data[j][i] = Median( val_list, list_length );
        }// i
    }// j

    delete[] val_list;
}

void WaveletTransformTest::initPicData( PicArray& pic_data )
{
    for (int j=pic_data.FirstY() ; j<=pic_data.LastY() ; ++j)
    {
       for (int i=pic_data.FirstX() ; i<=pic_data.LastX() ; ++i)
       {
           pic_data[j][i] = (((i-j) % 13)*1024)/13;       
       }// i
    }// j
}

ValueType DiagFilterBchkD( const PicArray& pic,
                           const int xpos, const int ypos,
                           const TwoDArray<int>& filter,
                           const int shift)
{
    // Half the filter length
    const int len2 = 6;

    const int height = pic.LengthY();
    const int width = pic.LengthX();

    int uplus, uneg, vplus, vneg;
    int val = (1<<(shift-1));

    // Do 0 position horizontally
    val += filter[0][0]*pic[ypos][xpos];

    for (int i=1; i<=len2;++i){
            uplus = xpos + i;
            uplus = (uplus>=width ? width-1 : uplus);
            uneg = xpos - i;
            uneg = (uneg<0 ? 0 : uneg );
            val += filter[0][i]*(pic[ypos][uplus]+pic[ypos][uneg] );
    }

    // Do other positions vertically//
    //////////////////////////////////

    for (int j=1; j<=len2;++j){
        vplus = ypos + j;
        vplus = ( vplus>=height ? height-1 : vplus);
        vneg = ypos - j;
        vneg = (vneg<0 ? 0 : vneg );

        // Do 0 position horizontally
        val += filter[j][0]*(pic[vneg][xpos]+pic[vplus][xpos]);
        for (int i=1; i<=len2;++i){
            uplus = xpos + i;
            uplus = (uplus>=width ? width-1 : uplus);
            uneg = xpos - i;
            uneg = (uneg<0 ? 0 : uneg );
            val += filter[j][i]*(pic[vneg][uplus]+pic[vneg][uneg]+
                                 pic[vplus][uplus]+pic[vplus][uneg] );
        }
    }

    val >>= shift;

    return ValueType(val);
}

void VFilter( PicArray& pic_data, const OneDArray<int>& filter, const int bits )
{
    PicArray tmp_data( pic_data );
    const int offset = (1<<(bits-1));

    int sum;

    // Do the first bit
    for (int j=0; j<filter.Last(); ++j)
    {

        for (int i=0; i<pic_data.LengthX(); ++i)
        {
            sum = offset;
            for (int k=filter.Last(); k>=filter.First(); --k)
                sum += filter[k]*pic_data[std::max(j-k,0)][i];
            sum >>= bits;
            sum = std::min( 127, std::max( -128, sum) );
            tmp_data[j][i] = ValueType( sum );
        }// i

    }// j

    // Do the middle bit
    for (int j=filter.Last(); j<=pic_data.LastY()+filter.First(); ++j)
    {

        for (int i=0; i<pic_data.LengthX(); ++i)
        {
            sum = offset;
            for (int k=filter.Last(); k>=filter.First(); --k)
                sum += filter[k]*pic_data[j-k][i];
            sum >>= bits;
            sum = std::min( 127, std::max( -128, sum) );
            tmp_data[j][i] = ValueType( sum );
        }// i

    }// j

    // Do the last bit
    for (int j=pic_data.LastY()+filter.First()+1; j<pic_data.LengthY(); ++j)
    {

        for (int i=0; i<pic_data.LengthX(); ++i)
        {
            sum = offset;
            for (int k=filter.Last(); k>=filter.First(); --k)
                sum += filter[k]*pic_data[std::min(j-k,pic_data.LastY())][i];
            sum >>= bits;
            sum = std::min( 127, std::max( -128, sum) );
            tmp_data[j][i] = ValueType( sum );
        }// i

    }// j

    // Copy data across
    pic_data = tmp_data;

}

// The loop over the columns is the same every time so lends itself to isolation
// as an individual function.
void DownConverter::RowLoop( const int colpos , const PicArray& old_data , PicArray& new_data)
{

     //Calculation variables
    int sum;
    const int xlen = 2*new_data.LengthX();
    int linepos=0;    

     // Leading Column Edge
     // Similar loops to the x case in ByHalf_opto, for explanation look there.
     // Note the factor of two difference as we only want to fill in every other
     // line as the others have already been created by the line loops.

    for( int x=0; x<(2*Stage_I_Size) ; x+=2 , linepos++ )
    {
        sum =  (m_row_buffer[((x)>=0)?(x):0]     + m_row_buffer[x+1])*StageI_I;
        sum += (m_row_buffer[((x-1)>=0)?(x-1):0] + m_row_buffer[x+2])*StageI_II;
        sum += (m_row_buffer[((x-2)>=0)?(x-2):0] + m_row_buffer[x+3])*StageI_III;
        sum += (m_row_buffer[((x-3)>=0)?(x-3):0] + m_row_buffer[x+4])*StageI_IV;
        sum += (m_row_buffer[((x-4)>=0)?(x-4):0] + m_row_buffer[x+5])*StageI_V;
        sum += (m_row_buffer[((x-5)>=0)?(x-5):0] + m_row_buffer[x+6])*StageI_VI;
        sum += 1<<(StageI_Shift-1);//do rounding right

        new_data[colpos][linepos] = sum >> StageI_Shift;

    }
     //Middle of column
    for( int x=(2*Stage_I_Size) ; x<xlen-(2*Stage_I_Size) ; x+=2 , linepos++) 
    {
        sum =  (m_row_buffer[x]   + m_row_buffer[x+1])*StageI_I;
        sum += (m_row_buffer[x-1] + m_row_buffer[x+2])*StageI_II;
        sum += (m_row_buffer[x-2] + m_row_buffer[x+3])*StageI_III;
        sum += (m_row_buffer[x-3] + m_row_buffer[x+4])*StageI_IV;
        sum += (m_row_buffer[x-4] + m_row_buffer[x+5])*StageI_V;
        sum += (m_row_buffer[x-5] + m_row_buffer[x+6])*StageI_VI;
        sum += 1<<(StageI_Shift-1);//do rounding right

        new_data[colpos][linepos] = sum >> StageI_Shift;

    }
     //Trailing column edge
    for( int x=xlen-(2*Stage_I_Size) ; x< xlen-1 ; x+=2 , linepos++ )
    {
        sum =  (m_row_buffer[x]   + m_row_buffer[((x+1)<xlen)?(x+1):(xlen-1)])*StageI_I;
        sum += (m_row_buffer[x-1] + m_row_buffer[((x+2)<xlen)?(x+2):(xlen-1)])*StageI_II;
        sum += (m_row_buffer[x-2] + m_row_buffer[((x+3)<xlen)?(x+3):(xlen-1)])*StageI_III;
        sum += (m_row_buffer[x-3] + m_row_buffer[((x+4)<xlen)?(x+4):(xlen-1)])*StageI_IV;
        sum += (m_row_buffer[x-4] + m_row_buffer[((x+5)<xlen)?(x+5):(xlen-1)])*StageI_V;
        sum += (m_row_buffer[x-5] + m_row_buffer[((x+6)<xlen)?(x+6):(xlen-1)])*StageI_VI;
        sum += 1<<(StageI_Shift-1);//do rounding right

        new_data[colpos][linepos] = sum >> StageI_Shift;

    }

}

BlockMatcher::BlockMatcher( const PicArray& pic_data , 
                            const PicArray& ref_data , 
                            const OLBParams& bparams ,
                            const int precision , 
                            const MvArray& mv_array , 
                            const TwoDArray< MvCostData >& cost_array):
    m_pic_data(pic_data), 
    m_ref_data(ref_data),
    m_mv_array(mv_array),
    m_cost_array(cost_array),
    m_peldiff( ref_data , pic_data ), //NB: ORDER!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    m_subpeldiff( 3 ),
    m_bparams( bparams ),
    m_var_max( (pic_data.LengthX()+pic_data.LengthY() )/216 ),
    m_var_max_up( (pic_data.LengthX()+pic_data.LengthY() )/27 ),
    m_precision( precision )
{
    m_subpeldiff[0] = new BlockDiffHalfPel( ref_data, pic_data ); 
    m_subpeldiff[1] = new BlockDiffQuarterPel( ref_data, pic_data );
    m_subpeldiff[2] = new BlockDiffEighthPel( ref_data, pic_data );
}

void HFilter( PicArray& pic_data, const OneDArray<int>& filter, const int bits )
{
    ValueType* line_data = new ValueType[pic_data.LengthX()];
    const int offset = (1<<(bits-1));

    int sum;

    for (int j=0; j<pic_data.LengthY(); ++j)
    {
        // Do the first bit
        for (int i=0; i<filter.Last(); ++i)
        {
            sum = offset;
            for (int k=filter.Last(); k>=filter.First(); --k)
                sum += filter[k]*pic_data[j][std::max(i-k,0)];
            sum >>= bits;
            sum = std::min( 127, std::max( -128, sum) );
            line_data[i] = ValueType( sum );
        }// i
        // Do the middle bit
        for (int i=filter.Last(); i<=pic_data.LastX()+filter.First(); ++i)
        {
            sum = offset;
            for (int k=filter.Last(); k>=filter.First(); --k)
                sum += filter[k]*pic_data[j][i-k];
            sum >>= bits;
            sum = std::min( 127, std::max( -128, sum) );
            line_data[i] = ValueType( sum );
        }// i
        // Do the last bit
        for (int i=pic_data.LastX()+filter.First()+1; i<pic_data.LengthX(); ++i)
        {
            sum = offset;
            for (int k=filter.Last(); k>=filter.First(); --k)
                sum += filter[k]*pic_data[j][std::min(i-k,pic_data.LastX())];
            sum >>= bits;
            sum = std::min( 127, std::max( -128, sum) );
            line_data[i] = ValueType( sum );
        }// i

        // Copy data back

        for (int i=0; i<pic_data.LengthX(); ++i )
            pic_data[j][i] = line_data[i];

    }// j

    delete[] line_data;
}

void WaveletTransform::Transform(const Direction d, PicArray& pic_data)
{
    int xl,yl; 

    if (d == FORWARD)
    {
        //do work
        xl=pic_data.LengthX(); 
        yl=pic_data.LengthY(); 
        
        for (int l = 1; l <= depth; ++l)
        {
            VHSplit(0,0,xl,yl,pic_data); 
            xl /= 2; 
            yl /= 2; 
        }

        band_list.Init( depth , pic_data.LengthX() , pic_data.LengthY() ); 
    }
    else
    {
        //do work
        xl = pic_data.LengthX()/(1<<(depth-1)); 
        yl = pic_data.LengthY()/(1<<(depth-1)); 
        
        for (int l = 1; l <= depth; ++l)
        {
            VHSynth(0,0,xl,yl,pic_data); 
            xl *= 2; 
            yl *= 2; 
        }
        
        //band list now inaccurate, so clear        
        band_list.Clear();     
    }
}

//Up-convert by a factor of two.
void UpConverter::DoUpConverter(const PicArray& pic_data, PicArray& up_data)
{

    m_width_old = std::min (pic_data.LengthX(), m_orig_xl);
    m_height_old = std::min (pic_data.LengthY(), m_orig_yl);
    m_width_new = std::min(2*m_width_old, up_data.LengthX());
    m_height_new = std::min(2*m_height_old, up_data.LengthY());

    // Filter params
    const int filter_size = 4;
    const int filter_shift = 5;
    const short taps[4] = {21,-7,3,-1}; 


    //Variables that will be used by the filter calculations
    ValueType sum;
    int ypos(0);
    

    //There are three y loops to cope with the leading edge, middle 
    //and trailing edge of each column.

    for(int y = 0 ; y < filter_size; ++y , ypos += 2)
    {

        //We are filtering each column but doing it bit by bit.
        //This means our main loop is in the x direction and
        //there is a much greater chance the data we need will
        //be in the cache.
        for(int x = 0 , xpos = 0; x < m_width_old; x++ , xpos+=2 )
        {
            // Copy a Pixel from the original image in each even position
            up_data[ypos][xpos] = pic_data[y][x];

            //Work out the next pixel from filtered values.
            //Excuse the complicated ternary stuff but it sorts out the edge
            sum  = 1 << (filter_shift-1);
            sum += (pic_data[y][x]              + pic_data[y+1][x]) * taps[0];
            sum += (pic_data[(y>=1)?(y-1):0][x] + pic_data[y+2][x]) * taps[1];
            sum += (pic_data[(y>=2)?(y-2):0][x] + pic_data[y+3][x]) * taps[2];
            sum += (pic_data[(y>=3)?(y-3):0][x] + pic_data[y+4][x]) * taps[3];

            sum >>= filter_shift;
            up_data[ypos+1][xpos] = CLIP(sum, m_min_val, m_max_val);
        }// x, xpos

        // The row loop.
        RowLoop( up_data , ypos, filter_size, filter_shift, taps );
    }// y, ypos
    // This loop is like the last one but it deals with the centre
    // section of the image and so the ternary operations are dropped
    // from the filter section.
    for(int y = filter_size; y < m_height_old - filter_size; ++y , ypos += 2)
    {
        for(int x = 0 , xpos=0; x < m_width_old; x++ , xpos+=2 )
        {
            up_data[ypos][xpos] = pic_data[y][x];

            sum  = 1 << (filter_shift-1);
            
            for (int t=0; t<filter_size; ++t)
                sum += (pic_data[y-t][x]   + pic_data[y+1+t][x]) * taps[t];

            sum >>= filter_shift;
            up_data[ypos+1][xpos] = CLIP(sum, m_min_val, m_max_val);
        }// x,xpos
        RowLoop( up_data , ypos, filter_size, filter_shift, taps );

    }// y, ypos 
    // Another similar loop! - this time we are dealing with
    // the trailing edge so the ternary stuff is back in the
    // filter calcs but in the second parameter.    
    for(int y = m_height_old - filter_size; y < m_height_old; ++y , ypos+=2)
    {
        for(int x = 0 , xpos=0 ; x < m_width_old; x++ , xpos+=2)
        {
            up_data[ypos][xpos]=pic_data[y][x];

            sum  = 1 << (filter_shift-1);
            sum += (pic_data[y][x]     + pic_data[((y+1)<m_height_old)?(y+1):(m_height_old-1)][x]) * taps[0];
            sum += (pic_data[y - 1][x] + pic_data[((y+2)<m_height_old)?(y+2):(m_height_old-1)][x]) * taps[1];
            sum += (pic_data[y - 2][x] + pic_data[((y+3)<m_height_old)?(y+3):(m_height_old-1)][x]) * taps[2];
            sum += (pic_data[y - 3][x] + pic_data[((y+4)<m_height_old)?(y+4):(m_height_old-1)][x]) * taps[3];

            sum >>= filter_shift;
            up_data[ypos+1][xpos] = CLIP(sum, m_min_val, m_max_val);
        }//x,xpos
        RowLoop( up_data , ypos, filter_size, filter_shift, taps );

    }//y,ypos
}

bool StreamFieldInput::ReadFieldComponent(PicArray& pic_data1,
                                   PicArray& pic_data2,
                                   const CompSort& cs)
{
    if (! *m_ip_pic_ptr)
        return false;

    //initially set up for 8-bit file input expanded to 10 bits for array output

    int xl,yl;
    if (cs == Y_COMP){
        xl = m_sparams.Xl();
        yl = m_sparams.Yl();
    }
    else{
        if (m_sparams.CFormat()==format420)
        {
            xl = m_sparams.Xl()/2;
            yl = m_sparams.Yl()/2;
        }
        else if (m_sparams.CFormat() == format422)
        {
            xl = m_sparams.Xl()/2;
            yl = m_sparams.Yl();
        }
        else{
            xl = m_sparams.Xl();
            yl = m_sparams.Yl();
        }
    }

    unsigned char * temp = new unsigned char[xl];//array big enough for one line
    ValueType *pic;

    for (int j=0 ; j<yl ; j++)
    {
        m_ip_pic_ptr->read((char*) temp, xl);
        if (j % 2 == 0)
        {
            pic = m_sparams.TopFieldFirst() ?
                    &pic_data1[j/2][0] : &pic_data2[j/2][0];
        }
        else
        {
            pic = m_sparams.TopFieldFirst() ?
                    &pic_data2[j/2][0] : &pic_data1[j/2][0];
        }
        for (int i=0 ; i<xl ; ++i)
        {
            pic[i] = (ValueType) temp[i];
        }//I
        for (int i=0 ; i<xl ; ++i)
        {
            pic[i] -= 128;
        }//I


        //pad the columns on the rhs using the edge value
        for (int i=xl ; i<pic_data1.LengthX() ; ++i ){
            pic[i] = pic[xl-1];
        }//I

    }//J

    delete [] temp;

    //now do the padded lines, using the last true line
    for (int j=yl/2 ; j<pic_data1.LengthY() ; ++j )
    {
        for (int i=0 ; i<pic_data1.LengthX() ; ++i )
        {
            pic_data1[j][i] = pic_data1[yl/2-1][i];
            pic_data2[j][i] = pic_data2[yl/2-1][i];
        }//I
    }//J

    return true;
}

bool StreamFrameInput::ReadFrameComponent(PicArray& pic_data, const CompSort& cs)
{

    if (! *m_ip_pic_ptr)
        return false;

    int xl,yl;
    if (cs == Y_COMP){
        xl = m_sparams.Xl();
        yl = m_sparams.Yl();
    }
    else{
        if (m_sparams.CFormat()==format420)
        {
            xl = m_sparams.Xl()/2;
            yl = m_sparams.Yl()/2;
        }
        else if (m_sparams.CFormat() == format422)
        {
            xl = m_sparams.Xl()/2;
            yl = m_sparams.Yl();
        }
        else{
            xl = m_sparams.Xl();
            yl = m_sparams.Yl();
        }
    }

    unsigned char * temp = new unsigned char[xl];//array big enough for one line

    for (int j=0 ; j<yl ; ++j)
    {
        m_ip_pic_ptr->read((char*) temp, xl);

        for (int i=0 ; i<xl ; ++i)
        {
            pic_data[j][i] = (ValueType) temp[i];
        }//I
        for (int i=0 ; i<xl ; ++i)
        {
            pic_data[j][i] -= 128;
        }//I


        //pad the columns on the rhs using the edge value
        for (int i=xl ; i<pic_data.LengthX() ; ++i ){
            pic_data[j][i] = pic_data[j][xl-1];
        }//I

    }//J

    delete [] temp;

    //now do the padded lines, using the last true line
    for (int j=yl ; j<pic_data.LengthY() ; ++j )
    {
        for (int i=0 ; i<pic_data.LengthX() ; ++i )
        {
            pic_data[j][i] = pic_data[yl-1][i];
        }//I
    }//J

    return true;
}

//General function - does some admin and calls the correct function
void DownConverter::DoDownConvert(const PicArray& old_data, PicArray& new_data)
{
    //Down-convert by a factor of two.
    m_row_buffer= new ValueType[old_data.LengthX()];    
    //Variables that will be used by the filter calculations
    int sum;
    int colpos;

    // The area of the picture that will be downconverted
    const int xlen = 2*new_data.LengthX();    
    const int ylen = 2*new_data.LengthY();    


    //There are three y loops to cope with the leading edge, middle 
    //and trailing edge of each column.
    colpos=0;
    for( int y=0; y<Stage_I_Size*2 ; y+=2 , colpos++ )
    {
        // We are filtering each column but doing it bit by bit.
        // This means our main loop is in the x direction and
        // there is a much greater chance the data we need will
        // be in the cache.

        for( int x=0 ; x<xlen ; x++ )
        {            
            // In down conversion we interpolate every pixel
            // so there is no copying.
            // Excuse the complicated ternary stuff but it sorts out the edge
            sum =  (old_data[y][x] + old_data[y+1][x])*StageI_I;
            sum += (old_data[((y-1)>=0)?(y-1):0][x] + old_data[y+2][x])*StageI_II;
            sum += (old_data[((y-2)>=0)?(y-2):0][x] + old_data[y+3][x])*StageI_III;
            sum += (old_data[((y-3)>=0)?(y-3):0][x] + old_data[y+4][x])*StageI_IV;
            sum += (old_data[((y-4)>=0)?(y-4):0][x] + old_data[y+5][x])*StageI_V;
            sum += (old_data[((y-5)>=0)?(y-5):0][x] + old_data[y+6][x])*StageI_VI;
            sum += 1<<(StageI_Shift-1);//do rounding right
            m_row_buffer[x] = sum >> StageI_Shift;
        }// x
        //Speaking of which - the row loop.

        RowLoop(colpos,old_data,new_data);
    }// y 

    // This loop is like the last one but it deals with the center
    // section of the image and so the ternary operations are dropped
    // from the filter section.
    for( int y=Stage_I_Size*2 ; y<ylen-Stage_I_Size*2 ; y+=2 , colpos++ )
    {
        for( int x=0 ; x<xlen ; x++ )
        {

            sum =  (old_data[y][x]   + old_data[y+1][x])*StageI_I;
            sum += (old_data[y-1][x] + old_data[y+2][x])*StageI_II;
            sum += (old_data[y-2][x] + old_data[y+3][x])*StageI_III;
            sum += (old_data[y-3][x] + old_data[y+4][x])*StageI_IV;
            sum += (old_data[y-4][x] + old_data[y+5][x])*StageI_V;
            sum += (old_data[y-5][x] + old_data[y+6][x])*StageI_VI;
            sum += 1<<(StageI_Shift-1);//do rounding right
            m_row_buffer[x] = sum >> StageI_Shift;
        }// x

        RowLoop( colpos , old_data , new_data );
    }// y

    // Another similar loop! - this time we are dealing with
    // the trailing edge so the ternary stuff is back in the
    // filter calcs but in the second parameter.

    for( int y=ylen-(Stage_I_Size*2) ; y<ylen-1 ; y+=2 , colpos++ )
    {
        for( int x=0; x<xlen ; x++ )
        {

            sum =  (old_data[y][x]   + old_data[((y+1)<ylen)?(y+1):(ylen-1)][x])*StageI_I;
            sum += (old_data[y-1][x] + old_data[((y+2)<ylen)?(y+2):(ylen-1)][x])*StageI_II;
            sum += (old_data[y-2][x] + old_data[((y+3)<ylen)?(y+3):(ylen-1)][x])*StageI_III;
            sum += (old_data[y-3][x] + old_data[((y+4)<ylen)?(y+4):(ylen-1)][x])*StageI_IV;
            sum += (old_data[y-4][x] + old_data[((y+5)<ylen)?(y+5):(ylen-1)][x])*StageI_V;
            sum += (old_data[y-5][x] + old_data[((y+6)<ylen)?(y+6):(ylen-1)][x])*StageI_VI;

            // Do rounding right
            sum += 1<<(StageI_Shift-1);
            m_row_buffer[x] = sum >> StageI_Shift;

        }// x

        RowLoop( colpos , old_data , new_data );

    }//  y

    // Tidy up the data
    delete[] m_row_buffer;

}

    /* 
    * NOTE: we are not doing any bounds checks here. This function must
    * be invoked only when the reference images start and stop fall
    * withing bounds
    */
    float simple_block_diff_up_mmx_4(
            const PicArray& pic_data, const PicArray& ref_data, 
            const ImageCoords& start_pos, const ImageCoords& end_pos, 
            const ImageCoords& ref_start, const ImageCoords& ref_stop,
            const MVector& rmdr, float cost_so_far, 
            float best_total_cost_so_far)
    {
        ValueType *pic_curr = &pic_data[start_pos.y][start_pos.x];
        ValueType *ref_curr = &ref_data[ref_start.y][ref_start.x];

        const int width = end_pos.x - start_pos.x;
        int height = end_pos.y - start_pos.y;
        const int ref_stride = ref_data.LengthX();

        // go down a row and back up
        const int pic_next = pic_data.LengthX() - width;
        // go down 2 rows and back up
        const int ref_next = ref_data.LengthX()*2 - width*2;

        REPORTM (ref_start.x>=0 && ref_stop.x < ref_data.LengthX() &&
                ref_start.y>=0 && ref_stop.y < ref_data.LengthY(), 
                "Reference image coordinates within bounds");

        CalcValueType sum = 0;
        CalcValueType mop_sum(0);
        int stopX = (width>>2)<<2;
        __m64 m_sum = _mm_set_pi16(0, 0, 0, 0);
        u_mmx_val u_sum;
        if (rmdr.x == 0 && rmdr.y == 0 )
        {
            //std::cerr << "Inmmx routine rmdr.x = rmdr.y = 0" << std::endl;
#if 1
            for( int y=0; y < height; y++, pic_curr+=pic_next, ref_curr+=ref_next )
            {
                m_sum = _mm_xor_si64 (m_sum, m_sum);
                mop_sum= 0;
                for( int x=0; x < stopX; x+=4, pic_curr+=4, ref_curr+=8 )
                {
                    __m64 pic = *(__m64 *)pic_curr;
                    __m64 ref = _mm_unpacklo_pi16 (*(__m64 *)ref_curr, *(__m64 *)(ref_curr+4));
                    __m64 ref2 = _mm_unpackhi_pi16 (*(__m64 *)ref_curr, *(__m64 *)(ref_curr+4));
                    ref = _mm_unpacklo_pi16 ( ref, ref2);
                    // ref - pic
                    pic = _mm_sub_pi16 (pic, ref);
                    // abs (ref - pic)
                    ref = _mm_srai_pi16(pic, 15);
                    pic = _mm_xor_si64(pic, ref);
                    pic = _mm_sub_pi16 (pic, ref);
                    // sum += abs(ref -pic)
                    /** 
                    * Since we are re-initialising m_sum with every loop
                    * maybe we don't need the following since overflow may
                    * not occur
                    ref = _mm_xor_si64(ref, ref);
                    ref = _mm_unpackhi_pi16(pic, ref);
                    pic = _mm_unpacklo_pi16(pic, pic);
                    pic = _mm_srai_pi32 (pic, 16);
                    pic = _mm_add_pi32 (pic, ref);
                    m_sum = _mm_add_pi32 (m_sum, pic);
                    **/
                    m_sum = _mm_add_pi16 (m_sum, pic);
                }
                // mopup;
                for (int x = stopX; x < width; ++x, ++pic_curr,ref_curr+=2)
                {
                    mop_sum += std::abs (*ref_curr - *pic_curr);
                }
                u_sum.m = m_sum;
                //sum += (u_sum.i[0] + u_sum.i[1] + mop_sum);
                sum += (u_sum.h[0] + u_sum.h[1] + u_sum.h[2] + u_sum.h[3] + mop_sum);
                _mm_empty();
                if ((sum + cost_so_far )>= best_total_cost_so_far)
                {
                    return best_total_cost_so_far;
                }
            }
            _mm_empty();
            return sum + cost_so_far;
#else
            float sum = cost_so_far;
            for( int y=0; y < height; ++y, pic_curr+=pic_next, ref_curr+=ref_next )
            {
                for( int x=0; x < width; ++x, ++pic_curr, ref_curr+=2 )
                {
                    sum += std::abs( *ref_curr - *pic_curr );
                }// x
                
                if ( sum>= best_total_cost_so_far)
                    return best_total_cost_so_far;

            }// y
            return sum;
#endif

        }
//.........这里部分代码省略.........

//General function - does some admin and calls the correct function
void DownConverter::DoDownConvert(const PicArray& old_data, PicArray& new_data)
{
    //Down-convert by a factor of two.
    m_row_buffer= new ValueType[old_data.LengthX()];    
    //Variables that will be used by the filter calculations
    int sum;
    int colpos;

    // The area of the picture that will be downconverted
    const int xlen = 2*new_data.LengthX();    
    const int ylen = 2*new_data.LengthY();    


    //There are three y loops to cope with the leading edge, middle 
    //and trailing edge of each column.
    colpos=0;

    static __m64 zero = _mm_set_pi16(0, 0, 0, 0);
    static __m64 tap0 = _mm_set_pi16 (0, StageI_I, 0, StageI_I);
    static __m64 tap1 = _mm_set_pi16 (0, StageI_II, 0, StageI_II);
    static __m64 tap2 = _mm_set_pi16 (0, StageI_III, 0, StageI_III);
    static __m64 tap3 = _mm_set_pi16 (0, StageI_IV, 0, StageI_IV);
    static __m64 tap4 = _mm_set_pi16 (0, StageI_V, 0, StageI_V);
    static __m64 tap5 = _mm_set_pi16 (0, StageI_VI, 0, StageI_VI);
    static __m64 round = _mm_set_pi32 ( 1<<(StageI_Shift-1), 1<<(StageI_Shift-1));

    u_sum sum1, sum2;
    __m64 tmp, m1, m2;

    int stopX = (xlen >> 2)<<2;
    for( int y=0; y<Stage_I_Size*2 ; y+=2 , colpos++ )
    {
        // We are filtering each column but doing it bit by bit.
        // This means our main loop is in the x direction and
        // there is a much greater chance the data we need will
        // be in the cache.

        for( int x=0 ; x<stopX ; x+=4 )
        {            
            // In down conversion we interpolate every pixel
            // so there is no copying.
            // Excuse the complicated ternary stuff but it sorts out the edge
            sum1.m = _mm_set_pi32 (0, 0);
            sum2.m = _mm_set_pi32 (0, 0);

            mmx_add (&old_data[y][x], &old_data[y+1][x], tap0, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[((y-1)>=0)?(y-1):0][x] , &old_data[y+2][x], tap1, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[((y-2)>=0)?(y-2):0][x] , &old_data[y+3][x], tap2, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[((y-3)>=0)?(y-3):0][x] , &old_data[y+4][x], tap3, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[((y-4)>=0)?(y-4):0][x] , &old_data[y+5][x], tap4, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[((y-5)>=0)?(y-5):0][x] , &old_data[y+6][x], tap5, zero, &sum1.m, &sum2.m);

            sum1.m = _mm_add_pi32 (sum1.m, round);
            sum2.m = _mm_add_pi32 (sum2.m, round);
            sum1.m = _mm_srai_pi32 (sum1.m, StageI_Shift);
            sum2.m = _mm_srai_pi32 (sum2.m, StageI_Shift);
            m_row_buffer[x] = sum1.i[0];
            m_row_buffer[x+1] = sum1.i[1];
            m_row_buffer[x+2] = sum2.i[0];
            m_row_buffer[x+3] = sum2.i[1];
        }// x
        _mm_empty();

        for( int x=stopX ; x<xlen ; x++ )
        {            
            // In down conversion we interpolate every pixel
            // so there is no copying.
            // Excuse the complicated ternary stuff but it sorts out the edge
            sum =  (old_data[y][x] + old_data[y+1][x])*StageI_I;
            sum += (old_data[((y-1)>=0)?(y-1):0][x] + old_data[y+2][x])*StageI_II;
            sum += (old_data[((y-2)>=0)?(y-2):0][x] + old_data[y+3][x])*StageI_III;
            sum += (old_data[((y-3)>=0)?(y-3):0][x] + old_data[y+4][x])*StageI_IV;
            sum += (old_data[((y-4)>=0)?(y-4):0][x] + old_data[y+5][x])*StageI_V;
            sum += (old_data[((y-5)>=0)?(y-5):0][x] + old_data[y+6][x])*StageI_VI;
            sum += 1<<(StageI_Shift-1);//do rounding right
            m_row_buffer[x] = sum >> StageI_Shift;
        }// x
        //Speaking of which - the row loop.

        RowLoop(colpos,new_data);
    }// y 

    // This loop is like the last one but it deals with the center
    // section of the image and so the ternary operations are dropped
    // from the filter section.
    for( int y=Stage_I_Size*2 ; y<ylen-Stage_I_Size*2 ; y+=2 , colpos++ )
    {
        for( int x=0 ; x<stopX ; x+=4 )
        {            
            // In down conversion we interpolate every pixel
            // so there is no copying.
            // Excuse the complicated ternary stuff but it sorts out the edge
            sum1.m = _mm_set_pi32 (0, 0);
            sum2.m = _mm_set_pi32 (0, 0);

            mmx_add (&old_data[y][x], &old_data[y+1][x], tap0, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[y-1][x] , &old_data[y+2][x], tap1, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[y-2][x] , &old_data[y+3][x], tap2, zero, &sum1.m, &sum2.m);
            mmx_add(&old_data[y-3][x] , &old_data[y+4][x], tap3, zero, &sum1.m, &sum2.m);
//.........这里部分代码省略.........

    CalcValueType simple_intra_block_diff_mmx_4 ( 
            const BlockDiffParams& dparams, 
            const PicArray& pic_data, ValueType &dc_val)
    {
        __m64 tmp = _mm_set_pi16(0, 0, 0, 0);
        u_mmx_val u_sum;
        u_sum.i[0] = u_sum.i[1] = 0;

        ValueType *src = &(pic_data[dparams.Yp()][dparams.Xp()]);

        int height = dparams.Yl();
        int width = dparams.Xl();
        int stopX = (width>>2)<<2;
        int pic_next = (pic_data.LengthX() - width);
        CalcValueType mop_sum = 0;
        for (int j = 0; j < height; j++)
        {
            for (int i = 0; i < stopX; i+=4) 
            {
                __m64 pic = *(__m64 *)src;
                // sum += (pic)
                tmp = _mm_xor_si64(tmp, tmp);
                tmp = _mm_unpackhi_pi16(pic, tmp);
                tmp = _mm_slli_pi32 (tmp, 16);
                tmp = _mm_srai_pi32 (tmp, 16);
                pic = _mm_unpacklo_pi16(pic, pic);
                pic = _mm_srai_pi32 (pic, 16);
                pic = _mm_add_pi32 (pic, tmp);
                u_sum.m = _mm_add_pi32 (u_sum.m, pic);
                src += 4;
            }
            // Mop up
            for (int i = stopX; i < width; ++i)
            {
                mop_sum += *src;
                src++;
            }
            src += pic_next;
        }

        CalcValueType int_dc =  (u_sum.i[0] + u_sum.i[1] + mop_sum)/(width*height);

        dc_val = static_cast<ValueType>( int_dc );

        // Now compute the resulting SAD
        __m64 dc = _mm_set_pi16 ( dc_val, dc_val , dc_val , dc_val);
        u_sum.m = _mm_xor_si64(u_sum.m, u_sum.m); // initialise sum to 0
        mop_sum = 0;
        
        src = &(pic_data[dparams.Yp()][dparams.Xp()]);
        for (int j = 0; j < height; ++j)
        {
            for (int i = 0; i < stopX; i+=4)
            {
                __m64 pic = *(__m64 *)src;
                // pic - dc
                pic = _mm_sub_pi16 (pic, dc);
                // abs (pic - dc)
                tmp = _mm_srai_pi16(pic, 15);
                pic = _mm_xor_si64(pic, tmp);
                pic = _mm_sub_pi16 (pic, tmp);
                // sum += abs(pic -dc)
                tmp = _mm_xor_si64(tmp, tmp);
                tmp = _mm_unpackhi_pi16(pic, tmp);
                pic = _mm_unpacklo_pi16(pic, pic);
                pic = _mm_srai_pi32 (pic, 16);
                pic = _mm_add_pi32 (pic, tmp);
                u_sum.m = _mm_add_pi32 (u_sum.m, pic);
                src += 4;
            }
            // Mop up
            for (int i = stopX; i < width; ++i)
            {
                mop_sum += std::abs(*src - dc_val);
                src++;
            }
            src += pic_next;
        }
        CalcValueType intra_cost = u_sum.i[0] + u_sum.i[1] + mop_sum;
        _mm_empty();

        return intra_cost;

    }