matrixUtil.h

   1 #ifndef MATRIXUTIL_HEADER_
   2 #define MATRIXUTIL_HEADER_
   3
   4 ////    Matrix templates
   5 /** Initializes column pointers so a linear field acts like a matrix */
   6 template<class T> inline void
   7 initMatrixPointers(int width,int height,T *linear,T **matrix) {
   8     for (T **mend=matrix+width; matrix!=mend; ++matrix,linear+=height)
   9         *matrix= linear;
  10 }
  11 /** Allocates a new matrix of type T with given width and height */
  12 template<class T> T** newMatrix(int width,int height) {
  13         assert( width>0 && height>0 );
  14     T **result= new T*[width];
  15     T *linear= new T[width*height];
  16     initMatrixPointers(width,height,linear,result);
  17     return result;
  18 }
  19 /** Releases a matrix of type T (zero is ignored) */
  20 template<class T> inline void delMatrix(T **m) {
  21     if (!m)
  22         return;
  23         delete[] m[0];
  24     delete[] m;
  25 }
  26 /** Fills a submatrix with a value */
  27 template<class T,class Block> void fillSubMatrix(T **m,const Block &block,T value) {
  28         T **mend= m+block.xend;
  29     for (m+=block.x0; m!=mend; ++m)
  30                 std::fill( *m+block.y0, *m+block.yend, value );
  31 }
  32
  33 /** MatrixSummer objects store partial sums of a matrix, allowing quick computation
  34  *      of sum of any rectangle in the matrix */
  35 template<class T,class U,class I> class MatrixSummer {
  36 public:
  37         typedef T result_type;
  38         /** The type of filling the summer - normal sums, sums of squares */
  39         enum SumType { Values=0, Squares=1 };
  40
  41 private:
  42         result_type **sums; ///< Internal matrix containing precomputed partial sums
  43
  44 public:
  45         /** Creates an empty summer */
  46         MatrixSummer()
  47         : sums(0) {}
  48         /** Creates a summer filled from a matrix */
  49         MatrixSummer( const U **m, SumType sumType, I width, I height )
  50         : sums(0)
  51                 { fill(m,sumType,width,height); }
  52         /** Frees the resources */
  53         ~MatrixSummer()
  54                 { free(); }
  55
  56         /** Only empty objects are allowed to be copied (assertion) */
  57         MatrixSummer( const MatrixSummer &DEBUG_ONLY(other) )
  58         : sums(0)
  59                 { assert(!other.sums); }
  60         /** Only empty objects are allowed to be assigned (assertion) */
  61         MatrixSummer& operator=( const MatrixSummer &DEBUG_ONLY(other) )
  62                 { return assert( !other.sums && !sums ), *this; }
  63
  64         /** Returns whether the object is filled with data */
  65         bool isValid() const
  66                 { return sums; }
  67         /** Clears the object */
  68         void invalidate()
  69                 { free(); };
  70
  71         /** Prepares object to make sums for a matrix. If the summer has already been used before,
  72          *      the method assumes it was for a matrix with the same size */
  73         void fill(const U **m,SumType sumType,I width,I height,I x0=0,I y0=0) {
  74                 assert( m && (sumType==Values || sumType==Squares) );
  75                 m+= x0;
  76                 if (!sums)
  77                         allocate(width,height);
  78                 I yEnd= height+y0;
  79         //      fill the edges with zeroes
  80                 for (I i=0; i<=width; ++i)
  81                         sums[i][0+y0]= 0;
  82                 for (I j=1+y0; j<=yEnd; ++j)
  83                         sums[0][j]= 0;
  84         //      acummulate in the y-growing direction
  85                 if (sumType==Values)
  86                         for (I i=1; i<=width; ++i)
  87                                 for (I j=1+y0; j<=yEnd; ++j)
  88                                         sums[i][j]= sums[i][j-1]+m[i-1][j-1];
  89                 else // sumType==Squares
  90                         for (I i=1; i<=width; ++i)
  91                                 for (I j=1+y0; j<=yEnd; ++j)
  92                                         sums[i][j]= sums[i][j-1]+sqr<result_type>(m[i-1][j-1]);
  93         //      acummulate in the x-growing direction
  94                 for (I i=2; i<=width; ++i)
  95                         for (I j=1+y0; j<=yEnd; ++j)
  96                                 sums[i][j]+= sums[i-1][j];
  97         }
  98         /** Just a shortcut for filling from a non-const matrix */
  99         void fill(U **m,SumType sumType,I width,I height,I x0=0,I y0=0)
 100                 { fill( bogoCast(m), sumType, width, height, x0, y0 ); }
 101
 102         /** Computes the sum of a rectangle (in constant time) */
 103         result_type getSum(I x0,I y0,I xend,I yend) const {
 104                 assert( sums && x0>=0 && y0>=0 && xend>=0 && yend>=0 && xend>x0 && yend>y0 );
 105                 return sums[xend][yend] -sums[x0][yend] -sums[xend][y0] +sums[x0][y0];
 106         }
 107 private:
 108         /** Allocates internal matrix */
 109         void allocate(I width,I height) {
 110                 assert( width>0 && height>0 );
 111                 delMatrix(sums);
 112                 sums= newMatrix<T>(width+1,height+1);
 113         }
 114         /** Frees all resources */
 115         void free() {
 116                 delMatrix(sums);
 117                 sums= 0;
 118         }
 119 };
 120
 121 namespace MatrixWalkers {
 122         using MTypes::Block;
 123
 124         /** Performs manipulations with rotation codes 0..7 (dihedral group of order eight) */
 125         struct Rotation {
 126                 /** Asserts the parameter is within 0..7 */
 127                 static void check(int DEBUG_ONLY(r)) {
 128                         assert( 0<=r && r<8 );
 129                 }
 130                 /** Returns inverted rotation (the one that takes this one back to identity) */
 131                 static int invert(int r) {
 132                         check(r);
 133                         return (4-r/2)%4 *2 +r%2;
 134                 }
 135                 /** Computes rotation equal to projecting at first with \p r1 and the result with \p r2 */
 136                 static int compose(int r1,int r2) {
 137                         check(r1); check(r2);
 138                         if (r2%2)
 139                                 r1= invert(r1);
 140                         return (r1/2 + r2/2) %4 *2+ ( (r1%2)^(r2%2) );
 141                 }
 142         };
 143
 144         /** Checked iterator for a rectangle in a matrix, no rotation */
 145         template<class T> class Checked {
 146         protected:
 147                 T **list                        ///  pointer the current column
 148                 , **listsEnd;           ///< pointer to the end column
 149                 size_t addElem          ///  offset of the first row from the top
 150                 , addElemsEnd;          ///< offset of the end row from the top
 151                 T *elem                         ///  pointer to the current element
 152                 , *elemsEnd;            ///< pointer to the end element of the current column
 153
 154         public:
 155                 Checked( T **pixels, const Block &block )
 156                 : list(pixels+block.x0), listsEnd(pixels+block.xend)
 157                 , addElem(block.y0), addElemsEnd(block.yend) {
 158                         DEBUG_ONLY( elem=elemsEnd=0; )
 159                         assert( block.xend>block.x0 && block.yend>block.y0 );
 160                 }
 161
 162                 bool outerCond()        { return list!=listsEnd; }
 163                 void outerStep()        { ++list; }
 164
 165                 void innerInit()        { elem= (*list)+addElem; elemsEnd= (*list)+addElemsEnd; }
 166                 bool innerCond()        { return elem!=elemsEnd; }
 167                 void innerStep()        { ++elem; }
 168
 169                 T& get()                        { return *elem; }
 170         };
 171
 172         /** Checked iterator for a whole QImage; no rotation, but always transposed */
 173         template<class T,class U> struct CheckedImage {
 174                 typedef T QImage;
 175                 typedef U QRgb;
 176         protected:
 177                 QImage &img;    ///< reference to the image
 178                 int lineIndex;  ///< index of the current line
 179                 QRgb *line              ///  pointer to the current pixel
 180                 , *lineEnd;             ///< pointer to the end of the line
 181
 182         public:
 183                 CheckedImage(QImage &image)
 184                 : img(image), lineIndex(0)
 185                         { DEBUG_ONLY(line=lineEnd=0;) }
 186
 187                 bool outerCond()        { return lineIndex<img.height(); }
 188                 void outerStep()        { ++lineIndex; }
 189
 190                 void innerInit()        { line= (QRgb*)img.scanLine(lineIndex); lineEnd= line+img.width(); }
 191                 bool innerCond()        { return line!=lineEnd; }
 192                 void innerStep()        { ++line; }
 193
 194                 QRgb& get()                     { return *line; }
 195         };
 196
 197         namespace NOSPACE {
 198                 /** Base structure for walkers changing 'x' in the outer and 'y' in the inner loop */
 199                 template<class T> class RBase_it {
 200                 protected:
 201                         T **list; ///< points to the beginning of the current column
 202                         const size_t addElem; ///< position of the first-to-be-walked elements from *#list[i]
 203                         T *elem; ///< points to the current element
 204
 205                 public:
 206                         RBase_it( T **listsBegin, size_t addElem_ )
 207                         : list(listsBegin), addElem(addElem_)
 208                                 { DEBUG_ONLY(elem=0;) }
 209
 210                         void innerInit()        { elem= *list+addElem; }
 211                         T& get()                        { return *elem; }
 212                 };
 213         }
 214
 215         /** No rotation: x->, y-> */
 216         template<class T> struct Rotation_0: public RBase_it<T> {
 217                 using RBase_it<T>::list; using RBase_it<T>::elem;
 218
 219                 Rotation_0( T **pixels, const Block &block )
 220                 : RBase_it<T>( pixels+block.x0, block.y0 ) {}
 221
 222                 void outerStep() { ++list; }
 223                 void innerStep() { ++elem; }
 224         };
 225
 226         /** Rotated 90deg. cw., transposed: x<-, y-> */
 227         template<class T> struct Rotation_1_T: public RBase_it<T> {
 228                 using RBase_it<T>::list; using RBase_it<T>::elem;
 229
 230                 Rotation_1_T( T **pixels, const Block &block )
 231                 : RBase_it<T>( pixels+block.xend-1, block.y0 ) {}
 232
 233                 void outerStep() { --list; }
 234                 void innerStep() { ++elem; }
 235         };
 236
 237         /** Rotated 180deg. cw.: x<-, y<- */
 238         template<class T> struct Rotation_2: public RBase_it<T> {
 239                 using RBase_it<T>::list; using RBase_it<T>::elem;
 240
 241                 Rotation_2( T **pixels, const Block &block )
 242                 : RBase_it<T>( pixels+block.xend-1, block.yend-1 ) {}
 243
 244                 void outerStep() { --list; }
 245                 void innerStep() { --elem; }
 246         };
 247
 248         /** Rotated 270deg. cw., transposed: x->, y<- */
 249         template<class T> struct Rotation_3_T: public RBase_it<T> {
 250                 using RBase_it<T>::list; using RBase_it<T>::elem;
 251
 252                 Rotation_3_T( T **pixels, const Block &block )
 253                 : RBase_it<T>( pixels+block.x0, block.yend-1 ) {}
 254
 255                 void outerStep() { ++list; }
 256                 void innerStep() { --elem; }
 257         };
 258
 259
 260         namespace NOSPACE {
 261                 /** Base structure for walkers changing 'y' in the outer and 'x' in the inner loop */
 262                 template<class T> class RBase_ix {
 263                 protected:
 264                         T *lineBegin; ///< points to the first element of the current line (first to be walked)
 265                         const size_t elemStep; ///< the number of elements to get the next column (positive)
 266                         T *elem; ///< points to the current element
 267
 268                 public:
 269                         RBase_ix( T **pixels, size_t x_start, size_t y_start )
 270                         : lineBegin(pixels[x_start]+y_start), elemStep(pixels[1]-pixels[0]) {
 271                                 DEBUG_ONLY(elem=0;)
 272                                 assert( size_t(pixels[2]-pixels[1]) == elemStep );
 273                         }
 274
 275                         void innerInit()        { elem= lineBegin; }
 276                         T& get()                        { return *elem; }
 277                 };
 278         }
 279
 280         /** No rotation, transposed: y->, x-> */
 281         template<class T> struct Rotation_0_T: public RBase_ix<T> {
 282                 using RBase_ix<T>::lineBegin; using RBase_ix<T>::elemStep; using RBase_ix<T>::elem;
 283
 284                 Rotation_0_T( T **pixels, const Block &block )
 285                 : RBase_ix<T>( pixels, block.x0, block.y0 ) {}
 286
 287                 void outerStep() { ++lineBegin; }
 288                 void innerStep() { elem+= elemStep; }
 289         };
 290
 291         /** Rotated 90deg. cw.: y->, x<- */
 292         template<class T> struct Rotation_1: public RBase_ix<T> {
 293                 using RBase_ix<T>::lineBegin; using RBase_ix<T>::elemStep; using RBase_ix<T>::elem;
 294
 295                 Rotation_1( T **pixels, const Block &block )
 296                 : RBase_ix<T>( pixels, block.xend-1, block.y0 ) {}
 297
 298                 void outerStep() { ++lineBegin; }
 299                 void innerStep() { elem-= elemStep; }
 300         };
 301
 302         /** Rotated 180deg. cw., transposed: y<-, x<- */
 303         template<class T> struct Rotation_2_T: public RBase_ix<T> {
 304                 using RBase_ix<T>::lineBegin; using RBase_ix<T>::elemStep; using RBase_ix<T>::elem;
 305
 306                 Rotation_2_T(T **pixels,const Block &block)
 307                 : RBase_ix<T>( pixels, block.xend-1, block.yend-1 ) {}
 308
 309                 void outerStep() { --lineBegin; }
 310                 void innerStep() { elem-= elemStep; }
 311         };
 312
 313         /** Rotated 270deg. cw.: y<-, x-> */
 314         template<class T> struct Rotation_3: public RBase_ix<T> {
 315                 using RBase_ix<T>::lineBegin; using RBase_ix<T>::elemStep; using RBase_ix<T>::elem;
 316
 317                 Rotation_3( T **pixels, const Block &block )
 318                 : RBase_ix<T>( pixels, block.x0, block.yend-1 ) {}
 319
 320                 void outerStep() { --lineBegin; }
 321                 void innerStep() { elem+= elemStep; }
 322         };
 323
 324
 325         template < class Check, class Unchecked, class Operator >
 326         Operator walkOperate( Check checked, Unchecked unchecked, Operator oper ) {
 327         //      outer cycle start - to be always run at least once
 328                 assert( checked.outerCond() );
 329                 do {
 330                 //      inner initialization
 331                         checked.innerInit();
 332                         unchecked.innerInit();
 333                 // inner cycle start - to be always run at least once
 334                         assert( checked.innerCond() );
 335                         do {
 336                         //      perform the operation and do the inner step for both iterators
 337                                 oper( checked.get(), unchecked.get() );
 338                                 checked.innerStep();
 339                                 unchecked.innerStep();
 340                         } while ( checked.innerCond() );
 341
 342                 //      signal the end of inner cycle to the operator and do the outer step for both iterators
 343                         oper.innerEnd();
 344                         checked.outerStep();
 345                         unchecked.outerStep();
 346
 347                 } while ( checked.outerCond() );
 348
 349                 return oper;
 350         }
 351
 352         template<class Check,class U,class Operator>// inline
 353         Operator walkOperateCheckRotate( Check checked, Operator oper
 354         , U **pixels2, const Block &block2, char rotation ) {
 355                 switch (rotation) {
 356                         case 0: return walkOperate( checked, Rotation_0  <U>(pixels2,block2) , oper );
 357                         case 1: return walkOperate( checked, Rotation_0_T<U>(pixels2,block2) , oper );
 358                         case 2: return walkOperate( checked, Rotation_1  <U>(pixels2,block2) , oper );
 359                         case 3: return walkOperate( checked, Rotation_1_T<U>(pixels2,block2) , oper );
 360                         case 4: return walkOperate( checked, Rotation_2  <U>(pixels2,block2) , oper );
 361                         case 5: return walkOperate( checked, Rotation_2_T<U>(pixels2,block2) , oper );
 362                         case 6: return walkOperate( checked, Rotation_3  <U>(pixels2,block2) , oper );
 363                         case 7: return walkOperate( checked, Rotation_3_T<U>(pixels2,block2) , oper );
 364                         default: assert(false); return oper;
 365                 }
 366         }
 367
 368         struct OperatorBase {
 369                 typedef MTypes::SReal SReal;
 370
 371                 void innerEnd() {}
 372         };
 373
 374         template<class CReal> struct RDSummer: public OperatorBase {
 375                 CReal totalSum, lineSum;
 376
 377                 RDSummer()
 378                 : totalSum(0), lineSum(0) {}
 379                 void operator()(const SReal &num1,const SReal& num2)
 380                         { lineSum+= CReal(num1) * CReal(num2); }
 381                 void innerEnd()
 382                         { totalSum+= lineSum; lineSum= 0; }
 383                 CReal result()
 384                         { assert(!lineSum); return totalSum; }
 385         };
 386
 387         template<class T> struct AddMulCopy: public OperatorBase {
 388                 const T toAdd, toMul;
 389
 390                 AddMulCopy(T add,T mul)
 391                 : toAdd(add), toMul(mul) {}
 392
 393                 template<class R1,class R2> void operator()(R1 &res,R2 f) const
 394                         { res= (f+toAdd)*toMul; }
 395         };
 396
 397         template<class T> struct MulAddCopyChecked: public OperatorBase {
 398                 const T toMul, toAdd, min, max;
 399
 400                 MulAddCopyChecked(T mul,T add,T minVal,T maxVal)
 401                 : toMul(mul), toAdd(add), min(minVal), max(maxVal) {}
 402                 template<class R1,class R2> void operator()(R1 &res,R2 f) const
 403                         { res= checkBoundsFunc( min, f*toMul+toAdd, max ); }
 404         };
 405 }
 406
 407 #endif // MATRIXUTIL_HEADER_