Click here to Skip to main content
15,896,496 members
Articles / Programming Languages / C

Linear Rank Deficient Filters in Computer Vision Applications

Rate me:
Please Sign up or sign in to vote.
4.84/5 (17 votes)
1 Jun 2008GPL34 min read 40.6K   1.4K   33  
The article describes the optimization of an image convolution operation with the application of linear rank deficient filters in computer vision problems.
#ifndef VEC2D_H
#define VEC2D_H


class vec2D       //16 byte aligned
{
public:
        enum CONV { SAME, VALID };

        vec2D(const wchar_t* file);
        vec2D(unsigned int ysize = 1, unsigned int xsize = 1,
              int yoffset = 0, int xoffset = 0, const float* data = 0);           // NxM vector,  data copied in [row1][row2] ... [rowN] form
        vec2D(const vec2D& v);
        ~vec2D();

// Operators
        inline const vec2D& operator=(const vec2D& v);                             //this = v        refs
        inline const vec2D& operator=(const float* pf);                            //this = pf
        inline float& operator()(int y, int x);                                    //operator v(y,x) = v.m_data[y][x]
        inline float operator()(int y, int x) const;                               //const operator v(y,x) = v.m_data[y][x]
        inline float& operator[](unsigned int i);                                  //operator v[i] = v.m_data[y][x]  0-offset operator MATLAB order like
        inline float operator[](unsigned int i) const;                             //const operator v[i] = v.m_data[y][x]  0-offset operator MATLAB order like

// Operations
        void print(const wchar_t* file = 0) const;                                 //dump contents
        inline float get(int y, int x) const;                                      //get periodic boundary conditions
        void set(float scalar);                                                    //set all array to a scalar
        void set(float scalar, RECT& r);                                           //set selected rect in array to a scalar [ left, right ) sse 0-offset
        void setrand();                                                            //set all array to rand values
        bool copy(const vec2D& v, int left = 0, int top = 0);                      //copy v's (this_size)[hxw] region from dx,dy offset to this 0-offset  [v>=this in size]
        void minmax(float& min, float& max) const;                                  //get min max values
        float maxval() const;                                                       //get max value
        float minval() const;                                                       //get min value
        void maxval(float& max, int& x, int& y,                                     //get max value with [x,y] coord
                    int sizex = 1, int sizey = 1, int dx = 0, int dy = 0) const;    //from offset dx,dy in rect [sizex,sizey]
        void fliplr();                                                              //flip matrix along vertical axis
        void flipud();                                                              //flip matrix along horizontal axis
        void rotate180();                                                           //rotate matrix by 180'
        vec2D* repmat(unsigned int v, unsigned int h);                              //replicate matrix v time vertically and h times horizontally


        //0-offset array mat operation//xfisrt,yfirst=0 xlast,ylast=width-1,height-1//////////////////////////////////////////////////////
        inline const vec2D& operator*(const vec2D& b) const;                    //new c = this*b      refs
        void add(const vec2D& a, const vec2D& b);                               //this = a.+b   sse optimized
        void sub(const vec2D& a, const vec2D& b);                               //this = a.-b   sse optimized
        void mul(const vec2D& a, const vec2D& b);                               //this = a*b
        void mult(const vec2D& a, const vec2D& b);                              //this = a*b'   sse optimized
        void div(const vec2D& a, const vec2D& b);                               //this = a./b   sse optimized
        void mule(const vec2D& a, const vec2D& b);                              //this = a.*b   sse optimized
        void add(float scalar);                                                 //this = this.+scalar   sse optimized
        void sub(float scalar);                                                 //this = this.-scalar   sse optimized
        void mul(float scalar);                                                 //this = this.*scalar   sse optimized
        void div(float scalar);                                                 //this = this./scalar   sse optimized
        float prod() const;                                                     //f = x0*x1*x2*...xN  sse optimized        
        void inter2(const vec2D& src, vec2D& dst_grdx, vec2D& dst_grdy);        //biliniar 2d interpolation, grd temp buffers width = this->size height=2
        void conv2(const vec2D& a, const vec2D& filter, enum CONV type = SAME); //this = a*filter  convolution with filter
        void conv2(const vec2D& a, const vec2D& re, const vec2D& im);           //this = convolution with complex filter [re i*im]        
        void conv2(const vec2D& a,                                              //this = sum(s(i) conv2( conv2(a, u[i]), v[i]) ))
                   const std::vector<vec2D>& u,                                 //col vectors                
                   const std::vector<vec2D>& v,                                 //row vectors
                   const vec2D& s,                                              //row vector of singular values
                   vec2D& tmp1,                                                 //VALID tmp1 = [(a.height - u.height + 1) x a.width]
                   vec2D& tmp2,                                                 //VALID tmp2 = this = [(a.height - u.height + 1) x a.width - a,width + 1]
                   enum CONV type = SAME);     
        float dot(const float* a, const float* b, unsigned int size) const;     //a0*b0+a1*b1+a2*b2+ ... An*Bn  sse 16-bit aligned
        //0-offset array mat operation//xfisrt,yfirst=0 xlast,ylast=width-1,height-1//////////////////////////////////////////////////////

        void normalize(float a, float b);                                       //normalize to a...b range  x1 = x-min
        void histeq(vec2D& hist);                                               //histogram normalization vector in 0...255 range -> to 0...1.0 range        


// Access
        inline int xfirst() const;                                                 //return first col index into array
        inline int yfirst() const;                                                 //return first row index into array
        inline int xlast() const;                                                  //return last col index into array
        inline int ylast() const;                                                  //return last row index into array
        inline unsigned int width() const;                                         //width of array
        inline unsigned int height() const;                                        //height size of array
        inline unsigned int length() const;                                        //total size of array width*height
        inline const float* data(int y, int x) const;


private:
        int m_xoffset, m_yoffset;
        int m_xlast, m_ylast;
        unsigned int m_width;
        unsigned int m_height;
        float** m_data;               //offseted data


        void init(unsigned int ysize, unsigned int xsize,
                  int yoffset = 0, int xoffset = 0);                            //allocate m_data only
        void init(const vec2D &v);                                              //fill this->m_data with v.m_data[:][:]
        void close(void);                                                       //deallocate m_data


};



/*
  arrange    -3 -2 -1 0 +1 +2  indeces array
          -3
          -2
          -1
           0          *
          +1
          +2

    vec1D v(6,6,-3,-3);
     v(-3,-3) = val1;
     v(-2,-3) = val2;
     ...
     v(2,2) = valN

*/


inline int vec2D::xfirst() const
{
        return m_xoffset;
}

inline int vec2D::yfirst() const
{
        return m_yoffset;
}

inline int vec2D::xlast() const
{
        return m_xlast;
}

inline int vec2D::ylast() const
{
        return m_ylast;
}

inline unsigned int vec2D::width() const
{
        return m_width;
}

inline unsigned int vec2D::height() const
{
        return m_height;
}

inline unsigned int vec2D::length() const
{
        return m_width*m_height;
}

inline const vec2D& vec2D::operator=(const float * pf)
{
        vec2D& v = *this;
        for (int y = v.yfirst(); y <= v.ylast(); y++)
                for (int x = v.xfirst(); x <= v.xlast(); x++)
                        v(y, x) = *pf++;

        return *this;
}

inline const vec2D& vec2D::operator=(const vec2D & v)
{
        if (this == &v) {
                return *this;
        } else if (m_width == v.width() && m_height == v.height()) {  //equal size arrays?
                if (m_xoffset != v.xfirst() || m_yoffset != v.yfirst()) {   //equal offsets? correct them if not
                        m_data += m_yoffset;
                        for (unsigned int j = 0; j < m_height; j++)
                                m_data[j] = m_data[j] + m_xoffset - v.xfirst();
                        m_data -= v.yfirst();

                        m_xoffset = v.xfirst();
                        m_yoffset = v.yfirst();
                }
                init(v);
        } else {  //create a complete copy from v
                close();
                init(v.height(), v.width(), v.yfirst(), v.xfirst());
                init(v);
        }
        return *this;
}

inline float& vec2D::operator()(int y, int x)
{
        return m_data[y][x];
}

inline float vec2D::operator()(int y, int x) const
{
        return m_data[y][x];
}

inline float& vec2D::operator[](unsigned int i)
{
        return m_data[i%m_height][i/m_height];     //m_height - MATLAB order operator v[:]
}

inline float vec2D::operator[](unsigned int i) const
{
        return m_data[i%m_height][i/m_height];
}

inline float vec2D::get(int y, int x) const
{
        if (x < xfirst()) x = xfirst() + (xfirst() - x);
        else if (x > xlast()) x = xlast() - (x - xlast());
        if (y < yfirst()) y = yfirst() + (yfirst() - y);
        else if (y > ylast()) y = ylast() - (y - ylast());

        return m_data[y][x];
}

inline const vec2D& vec2D::operator*(const vec2D& b) const     // c = a*b -> C = this*B
{
        const vec2D& a = *this;
        if (a.xfirst() != 0 || a.yfirst() != 0)
                return a;

        vec2D* pc = new vec2D(a.height(), b.width());             //0-offseted
        vec2D& c = *pc;

        for (unsigned int y = 0; y < c.height(); y++) {
                for (unsigned int x = 0; x < c.width(); x++) {
                        for (unsigned int i = 0; i < a.width(); i++)                     //a and b must be with zero offsets
                                c(y, x) += a(y, i) * b(i, x);
                }
        }

        return c;
}

inline const float* vec2D::data(int y, int x) const
{
        return m_data[y] + x;
}

inline float vec2D::dot(const float* a, const float* b, unsigned int size) const
{
        float z = 0.0f, fres = 0.0f;
        __declspec(align(16)) float ftmp[4] = { 0.0f, 0.0f, 0.0f, 0.0f };
        __m128 mres;

        if ((size / 4) != 0) {
                mres = _mm_load_ss(&z);
                for (unsigned int i = 0; i < size / 4; i++)
                        mres = _mm_add_ps(mres, _mm_mul_ps(_mm_load_ps(&a[4*i]), _mm_load_ps(&b[4*i])));

                //mres = a,b,c,d
                __m128 mv1 = _mm_movelh_ps(mres, mres);     //a,b,a,b
                __m128 mv2 = _mm_movehl_ps(mres, mres);     //c,d,c,d
                mres = _mm_add_ps(mv1, mv2);                //res[0],res[1]

                _mm_store_ps(ftmp, mres);

                fres = ftmp[0] + ftmp[1];
        }

        if ((size % 4) != 0) {
                for (unsigned int i = size - size % 4; i < size; i++)
                        fres += a[i] * b[i];
        }

        return fres;
}


#endif

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

License

This article, along with any associated source code and files, is licensed under The GNU General Public License (GPLv3)


Written By
Engineer
Russian Federation Russian Federation
Highly skilled Engineer with 14 years of experience in academia, R&D and commercial product development supporting full software life-cycle from idea to implementation and further support. During my academic career I was able to succeed in MIT Computers in Cardiology 2006 international challenge, as a R&D and SW engineer gain CodeProject MVP, find algorithmic solutions to quickly resolve tough customer problems to pass product requirements in tight deadlines. My key areas of expertise involve Object-Oriented
Analysis and Design OOAD, OOP, machine learning, natural language processing, face recognition, computer vision and image processing, wavelet analysis, digital signal processing in cardiology.

Comments and Discussions