1 #ifndef MATRIXUTIL_HEADER_
2 #define MATRIXUTIL_HEADER_
4 /** A simple structure representing a rectangle */
6 short x0
, y0
, xend
, yend
;
8 int width() const { return xend
-x0
; }
9 int height() const { return yend
-y0
; }
10 int size() const { return width()*height(); }
12 bool contains(short x
,short y
) const
13 { return x0
<=x
&& x
<xend
&& y0
<=y
&& y
<yend
; }
16 Block(short x0_
,short y0_
,short xend_
,short yend_
)
17 : x0(x0_
), y0(y0_
), xend(xend_
), yend(yend_
) {}
22 /** A simple generic template for matrices of fixed size, uses shallow copying
23 * and manual memory management */
24 template<class T
,class I
=PtrInt
> struct MatrixSlice
{
25 typedef MatrixSlice
<const T
,I
> Const
; ///< the class is convertible to Const type
27 T
*start
; ///< pointer to the top-left pixel
28 I colSkip
; ///< how many pixels to skip to get in the next column
30 /** Initializes an empty slice */
31 MatrixSlice(): start(0) {}
33 /** Creates a shallow copy (deleting one destroys all copies) */
34 MatrixSlice(const MatrixSlice
&m
): start(m
.start
), colSkip(m
.colSkip
) {}
36 /** Converts to a matrix of constant objects (shallow copy) */
37 operator Const() const {
40 result
.colSkip
= colSkip
;
44 /** Indexing operator - returns pointer to a column */
45 T
* operator[](I column
) {
47 return start
+column
*colSkip
;
49 /** Const version of indexing operator - doesn't allow to change the elements */
50 const T
* operator[](I column
) const {
51 return constCast(*this)[column
];
54 /** Reallocates the matrix for a new size. If \p memory parameter is given,
55 * it is used for storage (the user is responsible that the matrix fits in it, etc.) */
56 void allocate( I width
, I height
, T
*memory
=0 ) {
57 ASSERT( width
>0 && height
>0 );
59 start
= memory
? memory
: new T
[width
*height
];
62 /** Releases the memory */
67 /** Returns whether the matrix is allocated (and thus usable for indexing) */
68 bool isValid() const {
72 /** Fills a submatrix of a valid matrix with a value */
73 void fillSubMatrix(const Block
&block
,T value
) {
75 // compute begin and end column starts
76 T
*begin
= start
+block
.y0
+colSkip
*block
.x0
77 , *end
= start
+block
.y0
+colSkip
*block
.xend
;
79 for (T
*it
= begin
; it
!=end
; it
+= colSkip
)
80 std::fill( it
, it
+block
.height(), value
);
82 /** Shifts the indexing of this matrix - dangerous.
83 * After calling this, only addressing or more shifts can be done (not checked).
84 * Also shifts out of the allocated matrix aren't detected */
85 MatrixSlice
& shiftMatrix(I x0
,I y0
) {
87 start
+= x0
*colSkip
+y0
;
90 /** Computes relative position of a pointer in the matrix (always 0 <= \p y < #colSkip) */
91 void getPosition(const T
*elem
,int &x
,int &y
) const {
93 PtrInt diff
= elem
-start
;
98 x
= -((-diff
-1)/colSkip
) -1;
99 y
= colSkip
- (-diff
-1)%colSkip
-1;
102 }; // MatrixSlice class template
105 /** MatrixSummer objects store partial sums of a matrix, allowing quick computation
106 * of sum of any rectangle in the matrix. It's parametrized by "type of the result",
107 * "type of the input" and "indexing type" (defaults to Int) */
108 template<class T
,class I
=PtrInt
> struct MatrixSummer
{
111 MatrixSlice
<T
,I
> sums
; ///< Internal matrix containing precomputed partial sums
114 /** Creates an empty summer */
116 /** Only empty objects are allowed to be copied (assertion) */
117 MatrixSummer( const MatrixSummer
&other
)
118 { ASSERT( !other
.isValid() ); }
119 /** Only empty objects are allowed to be assigned (assertion) */
120 MatrixSummer
& operator=( const MatrixSummer
&other
)
121 { ASSERT( !other
.isValid() && !isValid() ); return *this; }
124 /** Returns whether the object is filled with data */
125 bool isValid() const { return sums
.isValid(); }
126 /** Clears the object */
127 void free() { sums
.free(); };
129 /** Computes the sum of a rectangle (in constant time) */
130 Result
getSum(I x0
,I y0
,I xend
,I yend
) const {
131 ASSERT( sums
.isValid() );
132 return sums
[xend
][yend
] -sums
[x0
][yend
] -sums
[xend
][y0
] +sums
[x0
][y0
];
134 /** A shortcut to get the sum of a block */
135 Result
getSum(const Block
&b
) const
136 { return getSum( b
.x0
, b
.y0
, b
.xend
, b
.yend
); }
138 /** Prepares object to make sums for a matrix. If the summer has already been
139 * used before, the method assumes it was for a matrix of the same size */
140 template<class Input
> void fill(Input inp
,I width
,I height
) {
141 if ( !sums
.isValid() )
142 sums
.allocate(width
+1,height
+1);
144 // fill the edges with zeroes
145 for (I i
=0; i
<=width
; ++i
)
147 for (I j
=1; j
<=height
; ++j
)
149 // acummulate in the y-growing direction
150 for (I i
=1; i
<=width
; ++i
)
151 for (I j
=1; j
<=height
; ++j
)
152 sums
[i
][j
]= sums
[i
][j
-1] + Result(inp
[i
-1][j
-1]);
153 // acummulate in the x-growing direction
154 for (I i
=2; i
<=width
; ++i
)
155 for (I j
=1; j
<=height
; ++j
)
156 sums
[i
][j
]+= sums
[i
-1][j
];
158 }; // MatrixSummer class template
160 /** Helper structure for computing with value and squared sums at once */
161 template<class Num
> struct DoubleNum
{
165 { DEBUG_ONLY( value
= square
= std::numeric_limits
<Num
>::quiet_NaN(); ) }
168 : value(val
), square(sqr(val
)) {}
170 DoubleNum(const DoubleNum
&other
)
171 : value(other
.value
), square(other
.square
) {}
173 void unpack(Num
&val
,Num
&sq
) const { val
= value
; sq
= square
; }
175 DoubleNum
& operator+=(const DoubleNum
&other
) {
177 square
+= other
.square
;
180 DoubleNum
& operator-=(const DoubleNum
&other
) {
182 square
-= other
.square
;
185 friend DoubleNum
operator+(const DoubleNum
&a
,const DoubleNum
&b
)
186 { return DoubleNum(a
)+= b
; }
187 friend DoubleNum
operator-(const DoubleNum
&a
,const DoubleNum
&b
)
188 { return DoubleNum(a
)-= b
; }
189 }; // DoubleNum template struct
192 /** Structure for a block of pixels - also contains summers and dimensions */
193 template< class SumT
, class PixT
, class I
=PtrInt
>
194 struct SummedMatrix
{
195 typedef DoubleNum
<SumT
> BSumRes
;
196 typedef MatrixSummer
<BSumRes
> BSummer
;
198 I width
/// The width of #pixels
199 , height
; ///< The height of #pixels
200 MatrixSlice
<PixT
> pixels
; ///< The matrix of pixels
201 BSummer summer
; ///< Summer for values and squares of #pixels
202 bool sumsValid
; ///< Indicates whether the summer values are valid
204 /** Sets the size of #pixels, optionally allocates memory */
205 void setSize( I width_
, I height_
) {
209 pixels
.allocate(width
,height
);
211 /** Frees the memory */
212 void free(bool freePixels
=true) {
221 /** Just validates both summers (if needed) */
222 void summers_makeValid() const {
223 ASSERT(pixels
.isValid());
225 constCast(summer
).fill(pixels
,width
,height
);
226 constCast(sumsValid
)= true;
229 /** Justs invalidates both summers (to be called after changes in the pixel-matrix) */
230 void summers_invalidate()
231 { sumsValid
= false; }
232 /** A shortcut for getting sums of a block */
233 BSumRes
getSums(const Block
&block
) const
234 { return getSums( block
.x0
, block
.y0
, block
.xend
, block
.yend
); }
235 /** Gets both sums of a nonempty rectangle in #pixels, the summer isn't validated */
236 BSumRes
getSums( I x0
, I y0
, I xend
, I yend
) const {
237 ASSERT( sumsValid
&& x0
>=0 && y0
>=0 && xend
>x0
&& yend
>y0
238 && xend
<=width
&& yend
<=height
);
239 return summer
.getSum(x0
,y0
,xend
,yend
);
241 }; // SummedPixels template struct
245 /** Contains various iterators for matrices (see Matrix)
246 * to be used in walkOperate() and walkOperateCheckRotate() */
247 namespace MatrixWalkers
{
249 /** Iterates two matrix iterators and performs an action.
250 * The loop is controled by the first iterator (\p checked)
251 * and on every corresponding pair (a,b) \p oper(a,b) is invoked. Returns \p oper. */
252 template < class Check
, class Unchecked
, class Operator
>
253 Operator
walkOperate( Check checked
, Unchecked unchecked
, Operator oper
) {
254 // outer cycle start - to be always run at least once
255 ASSERT( checked
.outerCond() );
257 // inner initialization
259 unchecked
.innerInit();
260 // inner cycle start - to be always run at least once
261 ASSERT( checked
.innerCond() );
263 // perform the operation and do the inner step for both iterators
264 oper( checked
.get(), unchecked
.get() );
266 unchecked
.innerStep();
267 } while ( checked
.innerCond() );
269 // signal the end of inner cycle to the operator and do the outer step for both iterators
272 unchecked
.outerStep();
274 } while ( checked
.outerCond() );
280 /** Base structure for walkers */
281 template<class T
,class I
> struct RotBase
{
283 typedef MatrixSlice
<T
,I
> TMatrix
;
285 TMatrix current
; ///< matrix starting on the current element
286 T
*lastStart
; ///< the place of the last enter of the inner loop
289 RotBase( TMatrix matrix
, int x0
, int y0
)
290 : current( matrix
.shiftMatrix(x0
,y0
) ), lastStart(current
.start
) {
291 DEBUG_ONLY( current
.start
= 0; )
292 ASSERT( matrix
.isValid() );
295 void innerInit() { current
.start
= lastStart
; }
296 T
& get() { return *current
.start
; }
297 }; // RotBase class template
299 #define ROTBASE_INHERIT \
300 typedef typename RotBase<T,I>::TMatrix TMatrix; \
301 using RotBase<T,I>::current; \
302 using RotBase<T,I>::lastStart;
304 /** No rotation: x->, y-> */
305 template<class T
,class I
> struct Rotation_0
: public RotBase
<T
,I
> { ROTBASE_INHERIT
306 Rotation_0( TMatrix matrix
, const Block
&block
)
307 : RotBase
<T
,I
>( matrix
, block
.x0
, block
.y0
) {}
309 void outerStep() { lastStart
+= current
.colSkip
; }
310 void innerStep() { ++current
.start
; }
313 /** Rotated 90deg\. cw\., transposed: x<-, y-> */
314 template<class T
,class I
> struct Rotation_1_T
: public RotBase
<T
,I
> { ROTBASE_INHERIT
315 Rotation_1_T( TMatrix matrix
, const Block
&block
)
316 : RotBase
<T
,I
>( matrix
, block
.xend
-1, block
.y0
) {}
318 void outerStep() { lastStart
-= current
.colSkip
; }
319 void innerStep() { ++current
.start
; }
322 /** Rotated 180deg\. cw\.: x<-, y<- */
323 template<class T
,class I
> struct Rotation_2
: public RotBase
<T
,I
> { ROTBASE_INHERIT
324 Rotation_2( TMatrix matrix
, const Block
&block
)
325 : RotBase
<T
,I
>( matrix
, block
.xend
-1, block
.yend
-1 ) {}
327 void outerStep() { lastStart
-= current
.colSkip
; }
328 void innerStep() { --current
.start
; }
331 /** Rotated 270deg\. cw\., transposed: x->, y<- */
332 template<class T
,class I
> struct Rotation_3_T
: public RotBase
<T
,I
> { ROTBASE_INHERIT
333 Rotation_3_T( TMatrix matrix
, const Block
&block
)
334 : RotBase
<T
,I
>( matrix
, block
.x0
, block
.yend
-1 ) {}
336 void outerStep() { lastStart
+= current
.colSkip
; }
337 void innerStep() { --current
.start
; }
340 /** No rotation, transposed: y->, x-> */
341 template<class T
,class I
> struct Rotation_0_T
: public RotBase
<T
,I
> { ROTBASE_INHERIT
342 Rotation_0_T( TMatrix matrix
, const Block
&block
)
343 : RotBase
<T
,I
>( matrix
, block
.x0
, block
.y0
) {}
345 void outerStep() { ++lastStart
; }
346 void innerStep() { current
.start
+= current
.colSkip
; }
349 /** Rotated 90deg\. cw\.: y->, x<- */
350 template<class T
,class I
> struct Rotation_1
: public RotBase
<T
,I
> { ROTBASE_INHERIT
351 Rotation_1( TMatrix matrix
, const Block
&block
)
352 : RotBase
<T
,I
>( matrix
, block
.xend
-1, block
.y0
) {}
354 void outerStep() { ++lastStart
; }
355 void innerStep() { current
.start
-= current
.colSkip
; }
358 /** Rotated 180deg\. cw\., transposed: y<-, x<- */
359 template<class T
,class I
> struct Rotation_2_T
: public RotBase
<T
,I
> { ROTBASE_INHERIT
360 Rotation_2_T( TMatrix matrix
, const Block
&block
)
361 : RotBase
<T
,I
>( matrix
, block
.xend
-1, block
.yend
-1 ) {}
363 void outerStep() { --lastStart
; }
364 void innerStep() { current
.start
-= current
.colSkip
; }
367 /** Rotated 270deg\. cw\.: y<-, x-> */
368 template<class T
,class I
> struct Rotation_3
: public RotBase
<T
,I
> { ROTBASE_INHERIT
369 Rotation_3( TMatrix matrix
, const Block
&block
)
370 : RotBase
<T
,I
>( matrix
, block
.x0
, block
.yend
-1 ) {}
372 void outerStep() { --lastStart
; }
373 void innerStep() { current
.start
+= current
.colSkip
; }
377 /** A flavour of walkOperate() choosing the right Rotation_* iterator based on \p rotation */
378 template<class Check
,class U
,class Operator
>
379 inline Operator
walkOperateCheckRotate( Check checked
, Operator oper
380 , MatrixSlice
<U
> pixels2
, const Block
&block2
, char rotation
) {
383 case 0: return walkOperate( checked
, Rotation_0
<U
,I
>(pixels2
,block2
) , oper
);
384 case 1: return walkOperate( checked
, Rotation_0_T
<U
,I
>(pixels2
,block2
) , oper
);
385 case 2: return walkOperate( checked
, Rotation_1
<U
,I
>(pixels2
,block2
) , oper
);
386 case 3: return walkOperate( checked
, Rotation_1_T
<U
,I
>(pixels2
,block2
) , oper
);
387 case 4: return walkOperate( checked
, Rotation_2
<U
,I
>(pixels2
,block2
) , oper
);
388 case 5: return walkOperate( checked
, Rotation_2_T
<U
,I
>(pixels2
,block2
) , oper
);
389 case 6: return walkOperate( checked
, Rotation_3
<U
,I
>(pixels2
,block2
) , oper
);
390 case 7: return walkOperate( checked
, Rotation_3_T
<U
,I
>(pixels2
,block2
) , oper
);
391 default: ASSERT(false); return oper
;
396 /** Performs manipulations with rotation codes 0-7 (dihedral group of order eight) */
398 /** Asserts the parameter is within 0-7 */
399 static void check(int DEBUG_ONLY(r
)) {
400 ASSERT( 0<=r
&& r
<8 );
402 /** Returns inverted rotation (the one that takes this one back to identity) */
403 static int invert(int r
) {
405 return (4-r
/2)%4 *2 +r
%2;
407 /** Computes rotation equal to projecting at first with \p r1 and the result with \p r2 */
408 static int compose(int r1
,int r2
) {
409 check(r1
); check(r2
);
412 return (r1
/2 + r2
/2) %4 *2+ ( (r1
%2)^(r2
%2) );
414 }; // Rotation struct
416 /** Checked_ iterator for a rectangle in a matrix, no rotation */
417 template<class T
,class I
=PtrInt
> struct Checked
: public Rotation_0
<T
,I
> {
418 typedef MatrixSlice
<T
,I
> TMatrix
;
419 typedef Rotation_0
<T
,I
> Base
;
421 using Base::lastStart
;
423 T
*colEnd
/// the end of the current column
424 , *colsEndStart
; ///< the start of the end column
426 /** Initializes a new iterator for a \p block of \p pixels */
427 Checked( TMatrix pixels
, const Block
&block
)
428 : Base( pixels
, block
), colEnd( pixels
.start
+pixels
.colSkip
*block
.x0
+block
.yend
)
429 , colsEndStart( pixels
.start
+pixels
.colSkip
*block
.xend
+block
.y0
) {
430 ASSERT( block
.xend
>block
.x0
&& block
.yend
>block
.y0
);
434 ASSERT(lastStart
<=colsEndStart
);
435 return lastStart
!=colsEndStart
;
438 colEnd
+= current
.colSkip
;
442 ASSERT(current
.start
<=colEnd
);
443 return current
.start
!=colEnd
;
445 }; // Checked class template
447 /** Checked_ iterator for a whole QImage; no rotation, but always transposed */
448 template<class T
,class U
> struct CheckedImage
{
452 QImage
&img
; ///< reference to the image
453 int lineIndex
/// index of the current line
454 , width
/// the width of the image
455 , height
; ///< the height of the image
456 QRgb
*line
/// pointer to the current pixel
457 , *lineEnd
; ///< pointer to the end of the line
460 /** Initializes a new iterator for an instance of QImage (Qt class) */
461 CheckedImage(QImage
&image
)
462 : img(image
), lineIndex(0), width(image
.width()), height(image
.height())
463 { DEBUG_ONLY(line
=lineEnd
=0;) }
465 bool outerCond() { return lineIndex
<height
; }
466 void outerStep() { ++lineIndex
; }
468 void innerInit() { line
= (QRgb
*)img
.scanLine(lineIndex
); lineEnd
= line
+width
; }
469 bool innerCond() { return line
!=lineEnd
; }
470 void innerStep() { ++line
; }
472 QRgb
& get() { return *line
; }
473 }; // CheckedImage class template
476 /** A convenience base type for operators to use with walkOperate() */
477 struct OperatorBase
{
481 /** An operator computing the sum of products */
482 template<class TOut
,class TIn
> struct RDSummer
: public OperatorBase
{
483 TOut totalSum
, lineSum
;
486 : totalSum(0), lineSum(0) {}
487 void operator()(const TIn
&num1
,const TIn
& num2
)
488 { lineSum
+= TOut(num1
) * TOut(num2
); }
490 { totalSum
+= lineSum
; lineSum
= 0; }
491 TOut
result() ///< returns the result
492 { ASSERT(!lineSum
); return totalSum
; }
495 /** An operator performing a= (b+::toAdd)*::toMul */
496 template<class T
> struct AddMulCopy
: public OperatorBase
{
497 const T toAdd
, toMul
;
499 AddMulCopy(T add
,T mul
)
500 : toAdd(add
), toMul(mul
) {}
502 template<class R1
,class R2
> void operator()(R1
&res
,R2 f
) const
503 { res
= (f
+toAdd
)*toMul
; }
506 /** An operator performing a= b*::toMul+::toAdd
507 * and moving the result into [::min,::max] bounds */
508 template<class T
> struct MulAddCopyChecked
: public OperatorBase
{
509 const T toMul
, toAdd
, min
, max
;
511 MulAddCopyChecked(T mul
,T add
,T minVal
,T maxVal
)
512 : toMul(mul
), toAdd(add
), min(minVal
), max(maxVal
) {}
514 template<class R1
,class R2
> void operator()(R1
&res
,R2 f
) const
515 { res
= checkBoundsFunc( min
, f
*toMul
+toAdd
, max
); }
517 } // MatrixWalkers namespace
519 #endif // MATRIXUTIL_HEADER_