matrixUtil.h

   1 #ifndef MATRIXUTIL_HEADER_
   2 #define MATRIXUTIL_HEADER_
   3
   4 #include "debug.h"
   5
   6
   7 /** Structure representing a rectangle */
   8 struct Block {
   9         short x0, y0, xend, yend;
  10
  11         int width()      const  { return xend-x0; }
  12         int height() const      { return yend-y0; }
  13         int size()       const  { return width()*height(); }
  14
  15         bool contains(short x,short y) const
  16                 { return x0<=x && x<xend && y0<=y && y<yend; }
  17
  18         Block() {}
  19         Block(short x0_,short y0_,short xend_,short yend_)
  20         : x0(x0_), y0(y0_), xend(xend_), yend(yend_) {}
  21 };
  22
  23 ////    Matrix templates
  24
  25 /** A simple generic template for matrices of fixed size, uses shallow copying
  26  *      and manual memory management */
  27 template<class T,class I=size_t> class Matrix {
  28 public:
  29         friend class Matrix< typename NonConstType<T>::Result, I >; ///< because of conversion
  30         typedef Matrix<const T,I> ConstMatrix; ///< the class is convertible to ConstMatrix
  31
  32 protected:
  33         T *start;       ///< pointer to the top-left pixel
  34         I colSkip;      ///< how many pixels to skip to get in the next column
  35         DEBUG_ONLY(I width;);
  36 public:
  37         /** Initializes an empty matrix */
  38         Matrix()
  39         : start(0) {}
  40
  41         /** Allocates a matrix of given size via #allocate (see for details) */
  42         Matrix( I width_, I height_, T *memory=0 )
  43         : start(0)
  44                 { allocate(width_,height_,memory); }
  45
  46         /** Creates a shallow copy */
  47         Matrix(const Matrix &m)
  48         : start(m.start), colSkip(m.colSkip)
  49                 { DEBUG_ONLY( width= m.width; ); }
  50
  51         /** Converts to a matrix of constant objects (shallow copy) */
  52         operator ConstMatrix() const {
  53                 ConstMatrix result;
  54                 result.start= start;
  55                 result.colSkip= colSkip;
  56                 DEBUG_ONLY( result.width= width; )
  57                 return result;
  58         }
  59
  60         /** Indexing operator - returns pointer to a column */
  61         T* operator[](I column) {
  62                 ASSERT( start && /*column>=0 &&*/ column<width );
  63                 return start+column*colSkip;
  64         }
  65         /** Const version of indexing operator - doesn't allow to change the elements */
  66         const T* operator[](I column) const {
  67                 return constCast(*this)[column];
  68         }
  69
  70         /** Reallocates the matrix for a new size. If \p memory parameter is given,
  71          *      it is used for storage (the user is responsible that the matrix fits in it, etc.) */
  72         void allocate( I width_, I height_, T *memory=0 ) {
  73                 free();
  74                 start= memory ? memory : new T[width_*height_];
  75                 colSkip= height_;
  76                 DEBUG_ONLY( width= width_; )
  77         }
  78         /** Releases the memory */
  79         void free() {
  80                 delete[] start;
  81                 start= 0;
  82         }
  83         /** Returns whether the matrix is allocated (and thus usable for indexing) */
  84         bool isValid() const {
  85                 return start;
  86         }
  87
  88         /** Fills a submatrix of a valid matrix with a value */
  89         void fillSubMatrix(const Block &block,T value) {
  90                 ASSERT( isValid() && T(block.xend)<=width );
  91         //      compute begin and end column starts
  92                 T *begin= start+block.y0+colSkip*block.x0
  93                 , *end= start+block.y0+colSkip*block.xend;
  94         //      fill the columns
  95             for (T *it= begin; it!=end; it+= colSkip)
  96                         std::fill( it, it+block.height(), value );
  97         }
  98 /** \name Low-level functions
  99  *      mainly for use by image iterators, etc.
 100  *      @{ */
 101         /** Shifts the indexing of this matrix - dangerous.
 102          *      After calling this, only addressing or more shifts can be done (not checked).
 103          *      Also shifts out of the allocated matrix aren't detected */
 104         void shiftMatrix(I x0,I y0) {
 105                 ASSERT( isValid() );
 106                 start+= x0+y0*colSkip;
 107                 DEBUG_ONLY( width-= x0; );
 108                 ASSERT(width>0);
 109         }
 110         /** Returns internal pointer #start */
 111         T* getStart()                   { return start; }
 112         /** Returns internal value #colSkip */
 113         I  getColSkip() const   { return colSkip; }
 114 ///     @}
 115 }; // Matrix class template
 116
 117
 118
 119 /** MatrixSummer objects store partial sums of a matrix, allowing quick computation
 120  *      of sum of any rectangle in the matrix. It's parametrized by "type of the result",
 121  *      "type of the input" and "indexing type" (defaults to size_t) */
 122 template<class T,class U,class I=size_t> class MatrixSummer {
 123 public:
 124         typedef T Result;
 125         typedef Matrix<const U,I> InpMatrix;
 126         /** The type of filling the summer - normal sums, sums of squares */
 127         enum SumType { Values=0, Squares=1 };
 128
 129 private:
 130         Matrix<T,I> sums; ///< Internal matrix containing precomputed partial sums
 131
 132 public:
 133         /** Creates an empty summer */
 134         MatrixSummer() {}
 135         /** Creates a summer filled from a matrix */
 136         MatrixSummer( InpMatrix m, SumType sumType, I width, I height )
 137                 { fill(m,sumType,width,height); }
 138         /** Frees the resources */
 139         ~MatrixSummer()
 140                 { sums.free(); }
 141
 142         /** Only empty objects are allowed to be copied (assertion) */
 143         MatrixSummer( const MatrixSummer &DEBUG_ONLY(other) )
 144                 { ASSERT( !other.isValid() ); }
 145         /** Only empty objects are allowed to be assigned (assertion) */
 146         MatrixSummer& operator=( const MatrixSummer &DEBUG_ONLY(other) )
 147                 { return ASSERT( !other.isValid() && !isValid() ), *this; }
 148
 149         /** Returns whether the object is filled with data */
 150         bool isValid() const
 151                 { return sums.isValid(); }
 152         /** Clears the object */
 153         void invalidate()
 154                 { sums.free(); };
 155
 156         /** Prepares object to make sums for a matrix. If the summer has already been used before,
 157          *      the method assumes it was for a matrix of the same size */
 158         void fill(InpMatrix m,SumType sumType,I width,I height,I x0=0,I y0=0) {
 159                 ASSERT( m.isValid() && (sumType==Values || sumType==Squares) );
 160
 161                 m.shiftMatrix(x0,y0);
 162
 163                 if ( !sums.isValid() )
 164                         sums.allocate(width+1,height+1);
 165
 166                 I yEnd= height+y0;
 167         //      fill the edges with zeroes
 168                 for (I i=0; i<=width; ++i)
 169                         sums[i][0]= 0;
 170                 for (I j=1; j<=height; ++j)
 171                         sums[0][j]= 0;
 172         //      acummulate in the y-growing direction
 173                 if (sumType==Values)
 174                         for (I i=1; i<=width; ++i)
 175                                 for (I j=1; j<=yEnd; ++j)
 176                                         sums[i][j]= sums[i][j-1]+m[i-1][j-1];
 177                 else // sumType==Squares
 178                         for (I i=1; i<=width; ++i)
 179                                 for (I j=1+y0; j<=yEnd; ++j)
 180                                         sums[i][j]= sums[i][j-1]+sqr<Result>(m[i-1][j-1]);
 181         //      acummulate in the x-growing direction
 182                 for (I i=2; i<=width; ++i)
 183                         for (I j=1+y0; j<=yEnd; ++j)
 184                                 sums[i][j]+= sums[i-1][j];
 185         }
 186
 187         /** Computes the sum of a rectangle (in constant time) */
 188         Result getSum(I x0,I y0,I xend,I yend) const {
 189                 ASSERT( sums.isValid() && xend>x0 && yend>y0 ); /*x0>=0 && y0>=0 && xend>=0 && yend>=0 &&*/
 190                 return sums[xend][yend] -sums[x0][yend] -sums[xend][y0] +sums[x0][y0];
 191         }
 192 }; // MatrixSummer class template
 193
 194 /** Contains various iterators for matrices (see Matrix)
 195  *      to be used in walkOperate() and walkOperateCheckRotate() */
 196 namespace MatrixWalkers {
 197
 198         /** Performs manipulations with rotation codes 0-7 (dihedral group of order eight) */
 199         struct Rotation {
 200                 /** Asserts the parameter is within 0-7 */
 201                 static void check(int DEBUG_ONLY(r)) {
 202                         ASSERT( 0<=r && r<8 );
 203                 }
 204                 /** Returns inverted rotation (the one that takes this one back to identity) */
 205                 static int invert(int r) {
 206                         check(r);
 207                         return (4-r/2)%4 *2 +r%2;
 208                 }
 209                 /** Computes rotation equal to projecting at first with \p r1 and the result with \p r2 */
 210                 static int compose(int r1,int r2) {
 211                         check(r1); check(r2);
 212                         if (r2%2)
 213                                 r1= invert(r1);
 214                         return (r1/2 + r2/2) %4 *2+ ( (r1%2)^(r2%2) );
 215                 }
 216         };
 217
 218         /** Checked_ iterator for a rectangle in a matrix, no rotation */
 219         template<class T,class I=size_t> class Checked {
 220         public:
 221                 typedef Matrix<T,I> TMatrix;
 222         protected:
 223                 I colSkip                       ///  offset of the "next row"
 224                 , height;                       ///< height of the iterated rectangle
 225                 T *col                          ///  pointer the current column
 226                 , *colEnd                       ///  pointer to the end column
 227                 , *elem                         ///  pointer to the current element
 228                 , *elemsEnd;            ///< pointer to the end element of the current column
 229
 230         public:
 231                 /** Initializes a new iterator for a \p block of \p pixels */
 232                 Checked( TMatrix pixels, const Block &block )
 233                 : colSkip       ( pixels.getColSkip() )
 234                 , height        ( block.height() )
 235                 , col           ( pixels.getStart()+block.y0+colSkip*block.x0 )
 236                 , colEnd        ( pixels.getStart()+block.y0+colSkip*block.xend ) {
 237                         DEBUG_ONLY( elem= elemsEnd= 0; )
 238                         ASSERT( block.xend>block.x0 && block.yend>block.y0 );
 239                 }
 240
 241                 bool outerCond()        { return col!=colEnd; }
 242                 void outerStep()        { col+= colSkip; }
 243
 244                 void innerInit()        { elem= col; elemsEnd= col+height; }
 245                 bool innerCond()        { return elem!=elemsEnd; }
 246                 void innerStep()        { ++elem; }
 247
 248                 T& get()                        { return *elem; }
 249         }; // Checked class template
 250
 251         /** Checked_ iterator for a whole QImage; no rotation, but always transposed */
 252         template<class T,class U> struct CheckedImage {
 253                 typedef T QImage;
 254                 typedef U QRgb;
 255         protected:
 256                 QImage &img;    ///< reference to the image
 257                 int lineIndex;  ///< index of the current line
 258                 QRgb *line              ///  pointer to the current pixel
 259                 , *lineEnd;             ///< pointer to the end of the line
 260
 261         public:
 262                 /** Initializes a new iterator for an instance of QImage (Qt class) */
 263                 CheckedImage(QImage &image)
 264                 : img(image), lineIndex(0)
 265                         { DEBUG_ONLY(line=lineEnd=0;) }
 266
 267                 bool outerCond()        { return lineIndex<img.height(); }
 268                 void outerStep()        { ++lineIndex; }
 269
 270                 void innerInit()        { line= (QRgb*)img.scanLine(lineIndex); lineEnd= line+img.width(); }
 271                 bool innerCond()        { return line!=lineEnd; }
 272                 void innerStep()        { ++line; }
 273
 274                 QRgb& get()                     { return *line; }
 275         }; // CheckedImage class template
 276
 277         /** Base structure for walkers changing 'x' in the outer and 'y' in the inner loop */
 278         template<class T,class I=size_t> struct RotBase {
 279         public:
 280                 typedef Matrix<T,I> MType;
 281         protected:
 282                 const I colSkip;        ///< the number of elements to skip to get on the next column
 283                 T *list                         ///  points to the beginning of the current column/row
 284                 , *elem;                        ///< points to the current element
 285
 286         public:
 287                 RotBase( MType matrix, I x0, I y0 )
 288                 : colSkip( matrix.getColSkip() ), list( matrix.getStart()+y0+colSkip*x0 )
 289                         { DEBUG_ONLY( elem= 0; ) }
 290
 291                 void innerInit()        { elem= list; }
 292                 T& get()                        { return *elem; }
 293         }; // RotBase class template
 294
 295         #define ROTBASE_INHERIT \
 296                 using RotBase<T,I>::colSkip; \
 297                 using RotBase<T,I>::list; \
 298                 using RotBase<T,I>::elem; \
 299                 typedef Matrix<T,I> MType;
 300
 301
 302         /** No rotation: x->, y-> */
 303         template<class T,class I> struct Rotation_0: public RotBase<T,I> { ROTBASE_INHERIT
 304                 Rotation_0( MType matrix, const Block &block )
 305                 : RotBase<T,I>( matrix, block.x0, block.y0 ) {}
 306
 307                 void outerStep() { list+= colSkip; }
 308                 void innerStep() { ++elem; }
 309         };
 310
 311         /** Rotated 90deg\. cw\., transposed: x<-, y-> */
 312         template<class T,class I> struct Rotation_1_T: public RotBase<T,I> { ROTBASE_INHERIT
 313                 Rotation_1_T( MType matrix, const Block &block )
 314                 : RotBase<T,I>( matrix, block.xend-1, block.y0 ) {}
 315
 316                 void outerStep() { list-= colSkip; }
 317                 void innerStep() { ++elem; }
 318         };
 319
 320         /** Rotated 180deg\. cw\.: x<-, y<- */
 321         template<class T,class I> struct Rotation_2: public RotBase<T,I> { ROTBASE_INHERIT
 322                 Rotation_2( MType matrix, const Block &block )
 323                 : RotBase<T,I>( matrix, block.xend-1, block.yend-1 ) {}
 324
 325                 void outerStep() { list-= colSkip; }
 326                 void innerStep() { --elem; }
 327         };
 328
 329         /** Rotated 270deg\. cw\., transposed: x->, y<- */
 330         template<class T,class I> struct Rotation_3_T: public RotBase<T,I> { ROTBASE_INHERIT
 331                 Rotation_3_T( MType matrix, const Block &block )
 332                 : RotBase<T,I>( matrix, block.x0, block.yend-1 ) {}
 333
 334                 void outerStep() { list+= colSkip; }
 335                 void innerStep() { --elem; }
 336         };
 337
 338         /** No rotation, transposed: y->, x-> */
 339         template<class T,class I> struct Rotation_0_T: public RotBase<T,I> { ROTBASE_INHERIT
 340                 Rotation_0_T( MType matrix, const Block &block )
 341                 : RotBase<T,I>( matrix, block.x0, block.y0 ) {}
 342
 343                 void outerStep() { ++list; }
 344                 void innerStep() { elem+= colSkip; }
 345         };
 346
 347         /** Rotated 90deg\. cw\.: y->, x<- */
 348         template<class T,class I> struct Rotation_1: public RotBase<T,I> { ROTBASE_INHERIT
 349                 Rotation_1( MType matrix, const Block &block )
 350                 : RotBase<T,I>( matrix, block.xend-1, block.y0 ) {}
 351
 352                 void outerStep() { ++list; }
 353                 void innerStep() { elem-= colSkip; }
 354         };
 355
 356         /** Rotated 180deg\. cw\., transposed: y<-, x<- */
 357         template<class T,class I> struct Rotation_2_T: public RotBase<T,I> { ROTBASE_INHERIT
 358                 Rotation_2_T( MType matrix, const Block &block )
 359                 : RotBase<T,I>( matrix, block.xend-1, block.yend-1 ) {}
 360
 361                 void outerStep() { --list; }
 362                 void innerStep() { elem-= colSkip; }
 363         };
 364
 365         /** Rotated 270deg\. cw\.: y<-, x-> */
 366         template<class T,class I> struct Rotation_3: public RotBase<T,I> { ROTBASE_INHERIT
 367                 Rotation_3( MType matrix, const Block &block )
 368                 : RotBase<T,I>( matrix, block.x0, block.yend-1 ) {}
 369
 370                 void outerStep() { --list; }
 371                 void innerStep() { elem+= colSkip; }
 372         };
 373
 374         /** Iterates two matrix iterators and performs an action.
 375          *      The loop is controled by the first iterator (\p checked)
 376          *      and on every corresponding pair (a,b) \p oper(a,b) is invoked. Returns \p oper. */
 377         template < class Check, class Unchecked, class Operator >
 378         Operator walkOperate( Check checked, Unchecked unchecked, Operator oper ) {
 379         //      outer cycle start - to be always run at least once
 380                 ASSERT( checked.outerCond() );
 381                 do {
 382                 //      inner initialization
 383                         checked.innerInit();
 384                         unchecked.innerInit();
 385                 // inner cycle start - to be always run at least once
 386                         ASSERT( checked.innerCond() );
 387                         do {
 388                         //      perform the operation and do the inner step for both iterators
 389                                 oper( checked.get(), unchecked.get() );
 390                                 checked.innerStep();
 391                                 unchecked.innerStep();
 392                         } while ( checked.innerCond() );
 393
 394                 //      signal the end of inner cycle to the operator and do the outer step for both iterators
 395                         oper.innerEnd();
 396                         checked.outerStep();
 397                         unchecked.outerStep();
 398
 399                 } while ( checked.outerCond() );
 400
 401                 return oper;
 402         }
 403
 404         /** A flavour of walkOperate() choosing the right Rotation_* iterator based on \p rotation */
 405         template<class Check,class U,class Operator> inline Operator walkOperateCheckRotate
 406         ( Check checked, Operator oper, Matrix<U> pixels2, const Block &block2, char rotation) {
 407                 typedef size_t I;
 408                 switch (rotation) {
 409                         case 0: return walkOperate( checked, Rotation_0  <U,I>(pixels2,block2) , oper );
 410                         case 1: return walkOperate( checked, Rotation_0_T<U,I>(pixels2,block2) , oper );
 411                         case 2: return walkOperate( checked, Rotation_1  <U,I>(pixels2,block2) , oper );
 412                         case 3: return walkOperate( checked, Rotation_1_T<U,I>(pixels2,block2) , oper );
 413                         case 4: return walkOperate( checked, Rotation_2  <U,I>(pixels2,block2) , oper );
 414                         case 5: return walkOperate( checked, Rotation_2_T<U,I>(pixels2,block2) , oper );
 415                         case 6: return walkOperate( checked, Rotation_3  <U,I>(pixels2,block2) , oper );
 416                         case 7: return walkOperate( checked, Rotation_3_T<U,I>(pixels2,block2) , oper );
 417                         default: ASSERT(false); return oper;
 418                 }
 419         }
 420
 421         /** A convenience base type for operators to use with walkOperate() */
 422         struct OperatorBase {
 423                 void innerEnd() {}
 424         };
 425
 426         /** An operator computing the sum of products */
 427         template<class TOut,class TIn> struct RDSummer: public OperatorBase {
 428                 TOut totalSum, lineSum;
 429
 430                 RDSummer()
 431                 : totalSum(0), lineSum(0) {}
 432                 void operator()(const TIn &num1,const TIn& num2)
 433                         { lineSum+= TOut(num1) * TOut(num2); }
 434                 void innerEnd()
 435                         { totalSum+= lineSum; lineSum= 0; }
 436                 TOut result()   ///< returns the result
 437                         { ASSERT(!lineSum); return totalSum; }
 438         };
 439
 440         /** An operator performing a= (b+#toAdd)*#toMul */
 441         template<class T> struct AddMulCopy: public OperatorBase {
 442                 const T toAdd, toMul;
 443
 444                 AddMulCopy(T add,T mul)
 445                 : toAdd(add), toMul(mul) {}
 446
 447                 template<class R1,class R2> void operator()(R1 &res,R2 f) const
 448                         { res= (f+toAdd)*toMul; }
 449         };
 450
 451         /** An operator performing a= b*#toMul+#toAdd and moving the result into [#min,#max] bounds */
 452         template<class T> struct MulAddCopyChecked: public OperatorBase {
 453                 const T toMul, toAdd, min, max;
 454
 455                 MulAddCopyChecked(T mul,T add,T minVal,T maxVal)
 456                 : toMul(mul), toAdd(add), min(minVal), max(maxVal) {}
 457                 template<class R1,class R2> void operator()(R1 &res,R2 f) const
 458                         { res= checkBoundsFunc( min, f*toMul+toAdd, max ); }
 459         };
 460 }
 461
 462 #endif // MATRIXUTIL_HEADER_