00001 #ifndef AREA_SUM_TABLE_QUADTREE_H
00002 #define AREA_SUM_TABLE_QUADTREE_H
00003
00004 #include <math.h>
00005 #include "dmath.h"
00006 #include "debug.h"
00007
00008 #include "Order2StatisticalConcept.h"
00009
00010 namespace za_co_codespot
00011 {
00012 namespace datastructures
00013 {
00014 using math::sqr;
00015
00030 template <typename T>
00031 class AreaSumTableQuadtree
00032 {
00033
00034 BOOST_CONCEPT_ASSERT((Order2StatisticalConcept<T>));
00035
00036
00037
00038 #ifdef CODE_SPOT_CO_ZA_DEBUG
00039 BOOST_CONCEPT_ASSERT((QuadtreeConcept<NaiveQuadtree<T>, T>));
00040 #endif
00041
00042 public:
00043 AreaSumTableQuadtree(const T ar[], unsigned int width, unsigned int height, const T & threshold);
00044 ~AreaSumTableQuadtree();
00045
00046 inline unsigned int getWidth() const;
00047 inline unsigned int getHeight() const;
00048
00049 const T& operator()(unsigned int x, unsigned int y) const;
00050 unsigned int getLevel(unsigned int x, unsigned int y) const;
00051 unsigned int getNodeCount() const;
00052
00053 private:
00054
00055 class Node
00056 {
00057 public:
00058 explicit Node();
00059 virtual ~Node() = 0;
00060 virtual inline const T& getValue(unsigned int x, unsigned int y) const = 0;
00061 virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const = 0;
00062 virtual unsigned int getNodeCount() const = 0;
00063 };
00064
00065 template <unsigned int childCount>
00066 class BranchNode : public Node
00067 {
00068 public:
00069 explicit BranchNode();
00070 virtual ~BranchNode();
00071 virtual unsigned int getNodeCount() const;
00072 protected:
00073 Node * mChildren[childCount];
00074 };
00075
00076 class RectangularBranchNode : public BranchNode<4>
00077 {
00078 public:
00079 RectangularBranchNode(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, const T & threshold_sqr);
00080 virtual inline const T& getValue(unsigned int x, unsigned int y) const;
00081 virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;
00082 private:
00083 unsigned int mMidX;
00084 unsigned int mMidY;
00085 };
00086
00087 class HorizontalBranchNode : public BranchNode<2>
00088 {
00089 public:
00090 HorizontalBranchNode(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, int unsigned x1, int unsigned y1, const T & threshold_sqr);
00091 virtual inline const T& getValue(unsigned int x, unsigned int y) const;
00092 virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;
00093 private:
00094 unsigned int mMidX;
00095 };
00096
00097 class VerticalBranchNode : public BranchNode<2>
00098 {
00099 public:
00100 VerticalBranchNode(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, const T & threshold_sqr);
00101 virtual inline const T& getValue(unsigned int x, unsigned int y) const;
00102 virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;
00103 private:
00104 unsigned int mMidY;
00105 };
00106
00107 class LeafNode : public Node
00108 {
00109 public:
00110 LeafNode(const T & value);
00111 virtual inline const T& getValue(unsigned int x, unsigned int y) const;
00112 virtual inline unsigned int getNodeCount() const;
00113 virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;
00114 private:
00115 T mValue;
00116 };
00117
00118 static inline T getVariance(const T sqr_sum_table[], unsigned int width, T mean, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1);
00119 static inline T getMean(const T sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1);
00120 static inline Node * initializeChild(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, const T & threshold_sqr);
00121 static inline const T& access(const T ar[], unsigned int x, unsigned int y, unsigned int width);
00122 static inline T& access(T ar[], unsigned int x, unsigned int y, unsigned int width);
00123 static void make_sum_tables(const T ar[], T sum_table[], T sqr_sum_table[], unsigned int width, unsigned int height);
00124
00125 unsigned int mWidth;
00126 unsigned int mHeight;
00127 Node * mRoot;
00128 };
00129
00130
00131
00132
00133 template <typename T>
00134 AreaSumTableQuadtree<T>::Node::Node()
00135 {
00136
00137 }
00138
00139 template <typename T>
00140 AreaSumTableQuadtree<T>::Node::~Node()
00141 {
00142
00143 }
00144
00145 template <typename T>
00146 template <unsigned int childCount>
00147 AreaSumTableQuadtree<T>::BranchNode<childCount>::BranchNode():
00148 Node()
00149 {
00150 for(unsigned int k = 0; k < childCount; k++)
00151 {
00152 mChildren[k] = 0;
00153 }
00154 }
00155
00156 template <typename T>
00157 template<unsigned int childCount>
00158 AreaSumTableQuadtree<T>::BranchNode<childCount>::~BranchNode()
00159 {
00160 for(unsigned int k = 0; k < childCount; k++)
00161 {
00162 delete mChildren[k];
00163 mChildren[k] = 0;
00164 }
00165 }
00166
00167 template <typename T>
00168 template<unsigned int childCount>
00169 unsigned int AreaSumTableQuadtree<T>::BranchNode<childCount>::getNodeCount() const
00170 {
00171 unsigned int sum = 1;
00172
00173 for(unsigned int k = 0; k < childCount; k++)
00174 {
00175 sum += mChildren[k]->getNodeCount();
00176 }
00177
00178 return sum;
00179 }
00180
00181 template <typename T>
00182 AreaSumTableQuadtree<T>::RectangularBranchNode::RectangularBranchNode(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, const T & threshold_sqr):
00183 BranchNode<4>(),
00184 mMidX((x0 + x1) / 2),
00185 mMidY((y0 + y1) / 2)
00186 {
00187 assert(x0 != x1 + 1);
00188 assert(y0 != y1 + 1);
00189
00190 mChildren[0] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, y0, mMidX, mMidY, threshold_sqr);
00191 mChildren[1] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, mMidX, y0, x1, mMidY, threshold_sqr);
00192 mChildren[2] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, mMidY, mMidX, y1, threshold_sqr);
00193 mChildren[3] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, mMidX, mMidY, x1, y1, threshold_sqr);
00194 }
00195
00196
00197
00198 template <typename T>
00199 const T& AreaSumTableQuadtree<T>::RectangularBranchNode::getValue(unsigned int x, unsigned int y) const
00200 {
00201 if (x < mMidX)
00202 {
00203 if (y < mMidY)
00204 {
00205 return mChildren[0]->getValue(x,y);
00206 }
00207 else
00208 {
00209 return mChildren[2]->getValue(x,y);
00210 }
00211 }
00212 else
00213 {
00214 if (y < mMidY)
00215 {
00216 return mChildren[1]->getValue(x,y);
00217 }
00218 else
00219 {
00220 return mChildren[3]->getValue(x,y);
00221 }
00222 }
00223 }
00224
00225 template <typename T>
00226 unsigned int AreaSumTableQuadtree<T>::RectangularBranchNode::getLevel(unsigned int x, unsigned int y) const
00227 {
00228 if (x < mMidX)
00229 {
00230 if (y < mMidY)
00231 {
00232 return mChildren[0]->getLevel(x,y) + 1;
00233 }
00234 else
00235 {
00236 return mChildren[2]->getLevel(x,y) + 1;
00237 }
00238 }
00239 else
00240 {
00241 if (y < mMidY)
00242 {
00243 return mChildren[1]->getLevel(x,y) + 1;
00244 }
00245 else
00246 {
00247 return mChildren[3]->getLevel(x,y) + 1;
00248 }
00249 }
00250 }
00251
00252
00253 template <typename T>
00254 AreaSumTableQuadtree<T>::HorizontalBranchNode::HorizontalBranchNode(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, const T & threshold_sqr):
00255 BranchNode<2>(),
00256 mMidX((x0 + x1) / 2)
00257 {
00258 assert(x0 != x1 + 1);
00259 assert(y0 != y1 + 1);
00260
00261 mChildren[0] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, y0, mMidX, y1, threshold_sqr);
00262 mChildren[1] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, mMidX, y0, x1, y1, threshold_sqr);
00263 }
00264
00265 template <typename T>
00266 const T& AreaSumTableQuadtree<T>::HorizontalBranchNode::getValue(unsigned int x, unsigned int y) const
00267 {
00268 if (x < mMidX)
00269 {
00270 return mChildren[0]->getValue(x,y);
00271 }
00272 else
00273 {
00274 return mChildren[1]->getValue(x,y);
00275 }
00276 }
00277
00278 template <typename T>
00279 unsigned int AreaSumTableQuadtree<T>::HorizontalBranchNode::getLevel(unsigned int x, unsigned int y) const
00280 {
00281 if (x < mMidX)
00282 {
00283 return mChildren[0]->getLevel(x,y) + 1;
00284 }
00285 else
00286 {
00287 return mChildren[1]->getLevel(x,y) + 1;
00288 }
00289 }
00290
00291
00292 template <typename T>
00293 AreaSumTableQuadtree<T>::VerticalBranchNode::VerticalBranchNode(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, const T & threshold_sqr):
00294 BranchNode<2>(),
00295 mMidY((y0 + y1) / 2)
00296 {
00297 assert(y0 != y1 + 1);
00298
00299 mChildren[0] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, y0, x1, mMidY, threshold_sqr);
00300 mChildren[1] = AreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, mMidY, x1, y1, threshold_sqr);
00301 }
00302
00303 template <typename T>
00304 const T& AreaSumTableQuadtree<T>::VerticalBranchNode::getValue(unsigned int x, unsigned int y) const
00305 {
00306 if (y < mMidY)
00307 {
00308 return mChildren[0]->getValue(x,y);
00309 }
00310 else
00311 {
00312 return mChildren[1]->getValue(x,y);
00313 }
00314 }
00315
00316 template <typename T>
00317 unsigned int AreaSumTableQuadtree<T>::VerticalBranchNode::getLevel(unsigned int x, unsigned int y) const
00318 {
00319 if (y < mMidY)
00320 {
00321 return mChildren[0]->getLevel(x,y) + 1;
00322 }
00323 else
00324 {
00325 return mChildren[1]->getLevel(x,y) + 1;
00326 }
00327 }
00328
00329
00330 template <typename T>
00331 AreaSumTableQuadtree<T>::LeafNode::LeafNode(const T & value):
00332 Node(),
00333 mValue(value)
00334 {}
00335
00336 template <typename T>
00337 const T& AreaSumTableQuadtree<T>::LeafNode::getValue(unsigned int x, unsigned int y) const
00338 {
00339 return mValue;
00340 }
00341
00342 template <typename T>
00343 unsigned int AreaSumTableQuadtree<T>::LeafNode::getLevel(unsigned int x, unsigned int y) const
00344 {
00345 return 0;
00346 }
00347
00348 template <typename T>
00349 unsigned int AreaSumTableQuadtree<T>::LeafNode::getNodeCount() const
00350 {
00351 return 1;
00352 }
00353
00354
00355 template <typename T>
00356 AreaSumTableQuadtree<T>::AreaSumTableQuadtree(const T ar[], unsigned int width, unsigned int height, const T & threshold):
00357 mWidth(width),
00358 mHeight(height)
00359 {
00360 T * sum_table = new T[width*height];
00361 T * sqr_sum_table = new T[width*height];
00362 make_sum_tables(ar, sum_table, sqr_sum_table, width, height);
00363
00364 mRoot = initializeChild(ar, sum_table, sqr_sum_table, width, 0, 0, width, height, sqr(threshold));
00365
00366 delete[] sum_table;
00367 delete[] sqr_sum_table;
00368 }
00369
00370 template <typename T>
00371 AreaSumTableQuadtree<T>::~AreaSumTableQuadtree()
00372 {
00373 delete mRoot;
00374 mRoot = 0;
00375 }
00376
00377 template <typename T>
00378 unsigned int AreaSumTableQuadtree<T>::getWidth() const
00379 {
00380 return mWidth;
00381 }
00382
00383 template <typename T>
00384 unsigned int AreaSumTableQuadtree<T>::getHeight() const
00385 {
00386 return mHeight;
00387 }
00388
00389 #ifndef access
00390 template <typename T>
00391 const T& AreaSumTableQuadtree<T>::access(const T ar[], unsigned int x, unsigned int y, unsigned int width)
00392 {
00393 return ar[y*width + x];
00394 }
00395
00396 template <typename T>
00397 T& AreaSumTableQuadtree<T>::access(T ar[], unsigned int x, unsigned int y, unsigned int width)
00398 {
00399 return ar[y*width + x];
00400 }
00401 #endif
00402
00403 template <typename T>
00404 void AreaSumTableQuadtree<T>::make_sum_tables(const T ar[], T sum_table[], T sqr_sum_table[], unsigned int width, unsigned int height)
00405 {
00406 access(sum_table, 0, 0, width) = access(ar, 0, 0, width);
00407 access(sqr_sum_table, 0, 0, width) = sqr(access(ar, 0, 0, width));
00408
00409 for(unsigned int i = 1; i < width; i++)
00410 {
00411 access(sum_table, i, 0, width) = access(sum_table, i - 1, 0, width) + access(ar, i, 0, width);
00412 access(sqr_sum_table, i, 0, width) = access(sqr_sum_table, i - 1, 0, width) + sqr(access(ar, i, 0, width));
00413 }
00414
00415 for(unsigned int j = 1; j < height; j++)
00416 {
00417 access(sum_table, 0, j, width) = access(sum_table, 0, j - 1, width) + access(ar, 0, j, width);
00418 access(sqr_sum_table, 0, j, width) = access(sqr_sum_table, 0, j - 1, width) + sqr(access(ar, 0, j, width));
00419 }
00420
00421 for(unsigned int i = 1; i < width; i++)
00422 {
00423 for(unsigned int j = 1; j < height; j++)
00424 {
00425 access(sum_table, i, j, width) =
00426 access(sum_table, i, j - 1, width)
00427 + access(sum_table, i - 1, j, width)
00428 - access(sum_table, i - 1, j - 1, width)
00429 + access(ar, i, j, width);
00430
00431 access(sqr_sum_table, i, j, width) =
00432 access(sqr_sum_table, i, j - 1, width)
00433 + access(sqr_sum_table, i - 1, j, width)
00434 - access(sqr_sum_table, i - 1, j - 1, width)
00435 + sqr(access(ar, i, j, width));
00436 }
00437 }
00438 }
00439
00440 template <typename T>
00441 T AreaSumTableQuadtree<T>::getMean(const T sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1)
00442 {
00443 assert(x0 <= x1);
00444 assert(y0 <= y1);
00445
00446 T sum;
00447
00448 if (x0 == 0)
00449 {
00450 if (y0 == 0)
00451 {
00452 sum = access(sum_table, x1, y1, width);
00453 }
00454 else
00455 {
00456 sum = access(sum_table, x1, y1, width) - access(sum_table, x1, y0 - 1, width);;
00457 }
00458 }
00459 else
00460 {
00461 if (y0 == 0)
00462 {
00463 sum = access(sum_table, x1, y1, width) - access(sum_table, x0 - 1, y1, width);
00464 }
00465 else
00466 {
00467 sum = access(sum_table, x1, y1, width) - access(sum_table, x1, y0 - 1, width)
00468 - access(sum_table, x0 - 1, y1, width) + access(sum_table, x0 - 1, y0 - 1, width);
00469 }
00470 }
00471
00472 return sum / ((x1 - x0 + 1) * (y1 - y0 + 1));
00473 }
00474
00475 template <typename T>
00476 T AreaSumTableQuadtree<T>::getVariance(const T sqr_sum_table[], unsigned int width, T mean, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1)
00477 {
00478 assert(x0 <= x1);
00479 assert(y0 <= y1);
00480
00481 T sum;
00482
00483 if (x0 == 0)
00484 {
00485 if (y0 == 0)
00486 {
00487 sum = access(sqr_sum_table, x1, y1, width);
00488 }
00489 else
00490 {
00491 sum = access(sqr_sum_table, x1, y1, width) - access(sqr_sum_table, x1, y0 - 1, width);;
00492 }
00493 }
00494 else
00495 {
00496 if (y0 == 0)
00497 {
00498 sum = access(sqr_sum_table, x1, y1, width) - access(sqr_sum_table, x0 - 1, y1, width);
00499 }
00500 else
00501 {
00502 sum = access(sqr_sum_table, x1, y1, width) - access(sqr_sum_table, x1, y0 - 1, width)
00503 - access(sqr_sum_table, x0 - 1, y1, width) + access(sqr_sum_table, x0 - 1, y0 - 1, width);
00504 }
00505 }
00506
00507 return sum / ((x1 - x0 + 1) * (y1 - y0 + 1)) - sqr(mean);
00508 }
00509
00510 template <typename T>
00511 const T& AreaSumTableQuadtree<T>::operator()(unsigned int x, unsigned int y) const
00512 {
00513 return mRoot->getValue(x,y);
00514 }
00515
00516 template <typename T>
00517 unsigned int AreaSumTableQuadtree<T>::getLevel(unsigned int x, unsigned int y) const
00518 {
00519 return mRoot->getLevel(x, y);
00520 }
00521
00522 template<typename T>
00523 typename AreaSumTableQuadtree<T>::Node * AreaSumTableQuadtree<T>::initializeChild(const T ar[], const T sum_table[], const T sqr_sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, const T & threshold_sqr)
00524 {
00525 if (x0 + 1 == x1)
00526 {
00527 if (y0 + 1 == y1)
00528 {
00529 return new LeafNode(access(ar, x0, y0, width));
00530 }
00531 else
00532 {
00533 return new VerticalBranchNode(ar, sum_table, sqr_sum_table, width, x0, y0, x1, y1, threshold_sqr);
00534 }
00535 }
00536 else
00537 {
00538 if (y0 + 1 == y1)
00539 {
00540 return new HorizontalBranchNode(ar, sum_table, sqr_sum_table, width, x0, y0, x1, y1, threshold_sqr);
00541 }
00542 }
00543
00544 T mean(getMean(sum_table, width, x0, y0, x1 - 1, y1 - 1));
00545 T variance(getVariance(sqr_sum_table, width, mean, x0, y0, x1 - 1, y1 - 1));
00546
00547 if(variance <= threshold_sqr)
00548 {
00549 return new LeafNode(mean);
00550 }
00551 else
00552 {
00553 return new RectangularBranchNode(ar, sum_table, sqr_sum_table, width, x0, y0, x1, y1, threshold_sqr);
00554 }
00555
00556 assert(0);
00557 }
00558
00559 template<typename T>
00560 unsigned int AreaSumTableQuadtree<T>::getNodeCount() const
00561 {
00562 return mRoot->getNodeCount();
00563 }
00564
00565 }
00566 }
00567
00568 #endif //AREA_SUM_TABLE_QUADTREE_H