matrixUtil.h

#ifndef MATRIXUTIL_HEADER_
#define MATRIXUTIL_HEADER_

struct Block {
    short x0, y0, xend, yend;

    int width()  const  { return xend-x0; }
    int height() const  { return yend-y0; }
    int size()   const  { return width()*height(); }

    bool contains(short x,short y) const
        { return x0<=x && x<xend && y0<=y && y<yend; }

    Block() {}
    Block(short x0_,short y0_,short xend_,short yend_)
    : x0(x0_), y0(y0_), xend(xend_), yend(yend_) {}
};


template<class T,class I=PtrInt> struct MatrixSlice {
    typedef MatrixSlice<const T,I> Const; 

    T *start;   
    I colSkip;  

    MatrixSlice(): start(0) {}

    static MatrixSlice makeRaw(T *start,I colSkip) {
        MatrixSlice result;
        result.start= start;
        result.colSkip= colSkip;
        return result;
    }

    operator Const() const {
        Const result;
        result.start= start;
        result.colSkip= colSkip;
        return result;
    }

    T* operator[](I column) {
        ASSERT( isValid() );
        return start+column*colSkip;
    }
    const T* operator[](I column) const {
        return constCast(*this)[column];
    }

    void allocate( I width, I height, T *memory=0 ) {
        ASSERT( width>0 && height>0 );
        free();
        start= memory ? memory : new T[width*height];
        colSkip= height;
    }
    void free() {
        delete[] start;
        start= 0;
    }
    bool isValid() const {
        return start;
    }

    void fillSubMatrix(const Block &block,T value) {
        ASSERT( isValid() );
    //  compute begin and end column starts
        T *begin= start+block.y0+colSkip*block.x0
        , *end= start+block.y0+colSkip*block.xend;
    //  fill the columns
        for (T *it= begin; it!=end; it+= colSkip)
            std::fill( it, it+block.height(), value );
    }
    MatrixSlice& shiftMatrix(I x0,I y0) {
        ASSERT( isValid() );
        start+= x0*colSkip+y0;
        return *this;
    }
    void getPosition(const T *elem,int &x,int &y) const {
        ASSERT( isValid() );
        PtrInt diff= elem-start;
        if (diff>=0) {
            x= diff/colSkip;
            y= diff%colSkip;
        } else {
            x= -((-diff-1)/colSkip) -1;
            y= colSkip - (-diff-1)%colSkip -1;
        }
    }
}; // MatrixSlice class template


template<class T,class I=PtrInt> struct MatrixSummer {
    typedef T Result;

    MatrixSlice<T,I> sums; 

    bool isValid()  const   { return sums.isValid(); }
    void free()             { sums.free(); };

    Result getSum(I x0,I y0,I xend,I yend) const {
        ASSERT( sums.isValid() );
        return sums[xend][yend] -(sums[x0][yend]+sums[xend][y0]) +sums[x0][y0];
    }
    Result getSum(const Block &b) const 
        { return getSum( b.x0, b.y0, b.xend, b.yend ); }

    template<class Input> void fill(Input inp,I width,I height) {
        if ( !sums.isValid() )
            sums.allocate(width+1,height+1);

    //  fill the edges with zeroes
        for (I i=0; i<=width; ++i)
            sums[i][0]= 0;
        for (I j=1; j<=height; ++j)
            sums[0][j]= 0;
    //  acummulate in the y-growing direction
        for (I i=1; i<=width; ++i)
            for (I j=1; j<=height; ++j)
                sums[i][j]= sums[i][j-1] + Result(inp[i-1][j-1]);
    //  acummulate in the x-growing direction
        for (I i=2; i<=width; ++i)
            for (I j=1; j<=height; ++j)
                sums[i][j]+= sums[i-1][j];
    }
}; // MatrixSummer class template

template<class Num> struct DoubleNum {
    Num value, square;
    
    DoubleNum() 
        { DEBUG_ONLY( value= square= std::numeric_limits<Num>::quiet_NaN(); ) }
    
    DoubleNum(Num val)
    : value(val), square(sqr(val)) {}
    
    DoubleNum(const DoubleNum &other)
    : value(other.value), square(other.square) {}
    
    void unpack(Num &val,Num &sq) const { val= value; sq= square; }
    
    DoubleNum& operator+=(const DoubleNum &other) {
        value+= other.value;
        square+= other.square;
        return *this;
    }
    DoubleNum& operator-=(const DoubleNum &other) {
        value-= other.value;
        square-= other.square;
        return *this;
    }
    friend DoubleNum operator+(const DoubleNum &a,const DoubleNum &b) 
        { return DoubleNum(a)+= b; }
    friend DoubleNum operator-(const DoubleNum &a,const DoubleNum &b) 
        { return DoubleNum(a)-= b; }
}; // DoubleNum template struct


template< class SumT, class PixT, class I=PtrInt >
struct SummedMatrix {
    typedef DoubleNum<SumT> BSumRes;
    typedef MatrixSummer<BSumRes,I> BSummer;
    
    I width                     
    , height;                   
    MatrixSlice<PixT> pixels;   
    BSummer summer;             
    bool sumsValid;             
    
    void setSize( I width_, I height_ ) {
        free();
        width= width_;
        height= height_;
        pixels.allocate(width,height);
    }
    void free(bool freePixels=true) {
        if (freePixels)
            pixels.free();
        else
            pixels.start= 0;
        summer.free();
        sumsValid= false;
    }
    
    void summers_makeValid() const {
        ASSERT(pixels.isValid());
        if (!sumsValid) {
            constCast(summer).fill(pixels,width,height);
            constCast(sumsValid)= true;
        }
    }
    void summers_invalidate() 
        { sumsValid= false; }
    BSumRes getSums(const Block &block) const 
        { return getSums( block.x0, block.y0, block.xend, block.yend ); }   
    BSumRes getSums( I x0, I y0, I xend, I yend ) const {
        ASSERT( sumsValid && x0>=0 && y0>=0 && xend>x0 && yend>y0 
            && xend<=width && yend<=height );
        return summer.getSum(x0,y0,xend,yend);
    }
    SumT getValueSum( I x0, I y0, I xend, I yend ) const {
        // \todo quicker version
        return getSums(x0,y0,xend,yend).value;
    }
}; // SummedPixels template struct



namespace MatrixWalkers {

    template < class Checked, class Unchecked, class Operator >
    Operator walkOperate( Checked checked, Unchecked unchecked, Operator oper ) {
    //  outer cycle start - to be always run at least once
        ASSERT( checked.outerCond() );
        do {
        //  inner initialization
            checked.innerInit();
            unchecked.innerInit();
        // inner cycle start - to be always run at least once
            ASSERT( checked.innerCond() );
            do {
            //  perform the operation and do the inner step for both iterators
                oper( checked.get(), unchecked.get() );
                checked.innerStep();
                unchecked.innerStep();
            } while ( checked.innerCond() );

        //  signal the end of inner cycle to the operator and do the outer step for both iterators
            oper.innerEnd();
            checked.outerStep();
            unchecked.outerStep();

        } while ( checked.outerCond() );

        return oper;
    }


    template<class T,class I> struct RotBase {
    public:
        typedef MatrixSlice<T,I> TMatrix;
    protected:
        TMatrix current;    
        T *lastStart;       

    public:
        RotBase( TMatrix matrix, int x0, int y0 )
        : current( matrix.shiftMatrix(x0,y0) ), lastStart(current.start) { 
            DEBUG_ONLY( current.start= 0; )
            ASSERT( matrix.isValid() );
        }

        void innerInit()    { current.start= lastStart; }
        T& get()            { return *current.start; }
    }; // RotBase class template

    #define ROTBASE_INHERIT \
        typedef typename RotBase<T,I>::TMatrix TMatrix; \
        using RotBase<T,I>::current; \
        using RotBase<T,I>::lastStart;

    template<class T,class I=PtrInt> struct Rotation_0: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_0( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.x0, block.y0 ) {}
        
        explicit Rotation_0( TMatrix matrix, I x0=0, I y0=0 )
        : RotBase<T,I>( matrix, x0, y0 ) {}

        void outerStep() { lastStart+= current.colSkip; }
        void innerStep() { ++current.start; }
    };

    template<class T,class I> struct Rotation_1_T: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_1_T( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.xend-1, block.y0 ) {}

        void outerStep() { lastStart-= current.colSkip; }
        void innerStep() { ++current.start; }
    };

    template<class T,class I> struct Rotation_2: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_2( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.xend-1, block.yend-1 ) {}

        void outerStep() { lastStart-= current.colSkip; }
        void innerStep() { --current.start; }
    };

    template<class T,class I> struct Rotation_3_T: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_3_T( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.x0, block.yend-1 ) {}

        void outerStep() { lastStart+= current.colSkip; }
        void innerStep() { --current.start; }
    };

    template<class T,class I> struct Rotation_0_T: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_0_T( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.x0, block.y0 ) {}
        
        explicit Rotation_0_T( TMatrix matrix, I x0=0, I y0=0 )
        : RotBase<T,I>( matrix, x0, y0 ) {}

        void outerStep() { ++lastStart; }
        void innerStep() { current.start+= current.colSkip; }
    };

    template<class T,class I> struct Rotation_1: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_1( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.xend-1, block.y0 ) {}

        void outerStep() { ++lastStart; }
        void innerStep() { current.start-= current.colSkip; }
    };

    template<class T,class I> struct Rotation_2_T: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_2_T( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.xend-1, block.yend-1 ) {}

        void outerStep() { --lastStart; }
        void innerStep() { current.start-= current.colSkip; }
    };

    template<class T,class I> struct Rotation_3: public RotBase<T,I> { ROTBASE_INHERIT
        Rotation_3( TMatrix matrix, const Block &block )
        : RotBase<T,I>( matrix, block.x0, block.yend-1 ) {}

        void outerStep() { --lastStart; }
        void innerStep() { current.start+= current.colSkip; }
    };

    
    template<class Check,class U,class Operator> 
    inline Operator walkOperateCheckRotate( Check checked, Operator oper
    , MatrixSlice<U> pixels2, const Block &block2, char rotation) {
        typedef PtrInt I;
        switch (rotation) {
            case 0: return walkOperate( checked, Rotation_0  <U,I>(pixels2,block2) , oper );
            case 1: return walkOperate( checked, Rotation_0_T<U,I>(pixels2,block2) , oper );
            case 2: return walkOperate( checked, Rotation_1  <U,I>(pixels2,block2) , oper );
            case 3: return walkOperate( checked, Rotation_1_T<U,I>(pixels2,block2) , oper );
            case 4: return walkOperate( checked, Rotation_2  <U,I>(pixels2,block2) , oper );
            case 5: return walkOperate( checked, Rotation_2_T<U,I>(pixels2,block2) , oper );
            case 6: return walkOperate( checked, Rotation_3  <U,I>(pixels2,block2) , oper );
            case 7: return walkOperate( checked, Rotation_3_T<U,I>(pixels2,block2) , oper );
            default: ASSERT(false); return oper;
        }
    }
    
    
    struct Rotation {
        static void check(int DEBUG_ONLY(r)) {
            ASSERT( 0<=r && r<8 );
        }
        static int invert(int r) {
            check(r);
            return (4-r/2)%4 *2 +r%2;
        }
        static int compose(int r1,int r2) {
            check(r1); check(r2);
            if (r2%2)
                r1= invert(r1);
            return (r1/2 + r2/2) %4 *2+ ( (r1%2)^(r2%2) );
        }
    }; // Rotation struct
    
    template<class T,class I=PtrInt> struct Checked: public Rotation_0<T,I> {
        typedef MatrixSlice<T,I> TMatrix;
        typedef Rotation_0<T,I> Base;
        using Base::current;
        using Base::lastStart;
        
        T *colEnd           
        , *colsEndStart;    

        Checked( TMatrix pixels, const Block &block )
        : Base( pixels, block ), colEnd( pixels.start+pixels.colSkip*block.x0+block.yend )
        , colsEndStart( pixels.start+pixels.colSkip*block.xend+block.y0 ) {
            ASSERT( block.xend>block.x0 && block.yend>block.y0 );
        }

        bool outerCond() {
            ASSERT(lastStart<=colsEndStart); 
            return lastStart!=colsEndStart; 
        }
        void outerStep() {
            colEnd+= current.colSkip; 
            Base::outerStep();
        }
        bool innerCond() {
            ASSERT(current.start<=colEnd);
            return current.start!=colEnd; 
        }
    }; // Checked class template

    template<class T,class U> struct CheckedImage {
        typedef T QImage;
        typedef U QRgb;
    protected:
        QImage &img;    
        int lineIndex   
        , width         
        , height;       
        QRgb *line      
        , *lineEnd;     

    public:
        CheckedImage(QImage &image)
        : img(image), lineIndex(0), width(image.width()), height(image.height())
            { DEBUG_ONLY(line=lineEnd=0;) }

        bool outerCond()    { return lineIndex<height; }
        void outerStep()    { ++lineIndex; }

        void innerInit()    { line= (QRgb*)img.scanLine(lineIndex); lineEnd= line+width; }
        bool innerCond()    { return line!=lineEnd; }
        void innerStep()    { ++line; }

        QRgb& get()         { return *line; }
    }; // CheckedImage class template
    

    struct OperatorBase {
        void innerEnd() {}
    };

    template<class TOut,class TIn> struct RDSummer: public OperatorBase {
        TOut totalSum, lineSum;

        RDSummer()
        : totalSum(0), lineSum(0) {}
        void operator()(const TIn &num1,const TIn& num2)
            { lineSum+= TOut(num1) * TOut(num2); }
        void innerEnd()
            { totalSum+= lineSum; lineSum= 0; }
        TOut result()   
            { ASSERT(!lineSum); return totalSum; }
    };

    template<class T> struct AddMulCopy: public OperatorBase {
        const T toAdd, toMul;

        AddMulCopy(T add,T mul)
        : toAdd(add), toMul(mul) {}

        template<class R1,class R2> void operator()(R1 &res,R2 f) const
            { res= (f+toAdd)*toMul; }
    };

    template<class T> struct MulAddCopyChecked: public OperatorBase {
        const T toMul, toAdd, min, max;

        MulAddCopyChecked(T mul,T add,T minVal,T maxVal)
        : toMul(mul), toAdd(add), min(minVal), max(maxVal) {}
        
        template<class R1,class R2> void operator()(R1 &res,R2 f) const
            { res= checkBoundsFunc( min, f*toMul+toAdd, max ); }
    };
} // MatrixWalkers namespace

#endif // MATRIXUTIL_HEADER_

---