00001 #ifndef AUGMENTED_AREA_SUM_TABLE_QUADTREE_H
00002 #define AUGMENTED_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
00032 template <typename T>
00033 class AugmentedAreaSumTableQuadtree
00034 {
00035
00036 BOOST_CONCEPT_ASSERT((Order2StatisticalConcept<T>));
00037
00038
00039
00040 #ifdef CODE_SPOT_CO_ZA_DEBUG
00041 BOOST_CONCEPT_ASSERT((QuadtreeConcept<NaiveQuadtree<T>, T>));
00042 #endif
00043
00044 public:
00045 AugmentedAreaSumTableQuadtree(const T ar[], unsigned int width, unsigned int height, const T & threshold);
00046 ~AugmentedAreaSumTableQuadtree();
00047
00048 inline unsigned int getWidth() const;
00049 inline unsigned int getHeight() const;
00050 const T& operator()(unsigned int x, unsigned int y) const;
00051 unsigned int getLevel(unsigned int x, unsigned int y) const;
00052 unsigned int getNodeCount() const;
00053
00054 private:
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 AugmentedAreaSumTableQuadtree<T>::Node::Node()
00135 {
00136
00137 }
00138
00139 template <typename T>
00140 AugmentedAreaSumTableQuadtree<T>::Node::~Node()
00141 {
00142
00143 }
00144
00145 template <typename T>
00146 template<unsigned int childCount>
00147 AugmentedAreaSumTableQuadtree<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 AugmentedAreaSumTableQuadtree<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 AugmentedAreaSumTableQuadtree<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 AugmentedAreaSumTableQuadtree<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] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, y0, mMidX, mMidY, threshold_sqr);
00191 mChildren[1] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, mMidX, y0, x1, mMidY, threshold_sqr);
00192 mChildren[2] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, mMidY, mMidX, y1, threshold_sqr);
00193 mChildren[3] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, mMidX, mMidY, x1, y1, threshold_sqr);
00194 }
00195
00196 template <typename T>
00197 const T& AugmentedAreaSumTableQuadtree<T>::RectangularBranchNode::getValue(unsigned int x, unsigned int y) const
00198 {
00199 if (x < mMidX)
00200 {
00201 if (y < mMidY)
00202 {
00203 return mChildren[0]->getValue(x,y);
00204 }
00205 else
00206 {
00207 return mChildren[2]->getValue(x,y);
00208 }
00209 }
00210 else
00211 {
00212 if (y < mMidY)
00213 {
00214 return mChildren[1]->getValue(x,y);
00215 }
00216 else
00217 {
00218 return mChildren[3]->getValue(x,y);
00219 }
00220 }
00221 }
00222
00223 template <typename T>
00224 unsigned int AugmentedAreaSumTableQuadtree<T>::RectangularBranchNode::getLevel(unsigned int x, unsigned int y) const
00225 {
00226 if (x < mMidX)
00227 {
00228 if (y < mMidY)
00229 {
00230 return mChildren[0]->getLevel(x,y) + 1;
00231 }
00232 else
00233 {
00234 return mChildren[2]->getLevel(x,y) + 1;
00235 }
00236 }
00237 else
00238 {
00239 if (y < mMidY)
00240 {
00241 return mChildren[1]->getLevel(x,y) + 1;
00242 }
00243 else
00244 {
00245 return mChildren[3]->getLevel(x,y) + 1;
00246 }
00247 }
00248 }
00249
00250
00251 template <typename T>
00252 AugmentedAreaSumTableQuadtree<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):
00253 BranchNode<2>(),
00254 mMidX((x0 + x1) / 2)
00255 {
00256 assert(x0 != x1 + 1);
00257 assert(y0 != y1 + 1);
00258
00259 mChildren[0] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, y0, mMidX, y1, threshold_sqr);
00260 mChildren[1] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, mMidX, y0, x1, y1, threshold_sqr);
00261 }
00262
00263 template <typename T>
00264 const T& AugmentedAreaSumTableQuadtree<T>::HorizontalBranchNode::getValue(unsigned int x, unsigned int y) const
00265 {
00266 if (x < mMidX)
00267 {
00268 return mChildren[0]->getValue(x,y);
00269 }
00270 else
00271 {
00272 return mChildren[1]->getValue(x,y);
00273 }
00274 }
00275
00276 template <typename T>
00277 unsigned int AugmentedAreaSumTableQuadtree<T>::HorizontalBranchNode::getLevel(unsigned int x, unsigned int y) const
00278 {
00279 if (x < mMidX)
00280 {
00281 return mChildren[0]->getLevel(x,y) + 1;
00282 }
00283 else
00284 {
00285 return mChildren[1]->getLevel(x,y) + 1;
00286 }
00287 }
00288
00289
00290 template <typename T>
00291 AugmentedAreaSumTableQuadtree<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):
00292 BranchNode<2>(),
00293 mMidY((y0 + y1) / 2)
00294 {
00295 assert(y0 != y1 + 1);
00296
00297 mChildren[0] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, y0, x1, mMidY, threshold_sqr);
00298 mChildren[1] = AugmentedAreaSumTableQuadtree<T>::initializeChild(ar, sum_table, sqr_sum_table, width, x0, mMidY, x1, y1, threshold_sqr);
00299 }
00300
00301 template <typename T>
00302 const T& AugmentedAreaSumTableQuadtree<T>::VerticalBranchNode::getValue(unsigned int x, unsigned int y) const
00303 {
00304 if (y < mMidY)
00305 {
00306 return mChildren[0]->getValue(x,y);
00307 }
00308 else
00309 {
00310 return mChildren[1]->getValue(x,y);
00311 }
00312 }
00313
00314 template <typename T>
00315 unsigned int AugmentedAreaSumTableQuadtree<T>::VerticalBranchNode::getLevel(unsigned int x, unsigned int y) const
00316 {
00317 if (y < mMidY)
00318 {
00319 return mChildren[0]->getLevel(x,y) + 1;
00320 }
00321 else
00322 {
00323 return mChildren[1]->getLevel(x,y) + 1;
00324 }
00325 }
00326
00327
00328 template <typename T>
00329 AugmentedAreaSumTableQuadtree<T>::LeafNode::LeafNode(const T & value):
00330 Node(),
00331 mValue(value)
00332 {}
00333
00334 template <typename T>
00335 const T& AugmentedAreaSumTableQuadtree<T>::LeafNode::getValue(unsigned int x, unsigned int y) const
00336 {
00337 return mValue;
00338 }
00339
00340 template <typename T>
00341 unsigned int AugmentedAreaSumTableQuadtree<T>::LeafNode::getLevel(unsigned int x, unsigned int y) const
00342 {
00343 return 0;
00344 }
00345
00346 template <typename T>
00347 unsigned int AugmentedAreaSumTableQuadtree<T>::LeafNode::getNodeCount() const
00348 {
00349 return 1;
00350 }
00351
00352
00353 template <typename T>
00354 AugmentedAreaSumTableQuadtree<T>::AugmentedAreaSumTableQuadtree(const T ar[], unsigned int width, unsigned int height, const T & threshold):
00355 mWidth(width),
00356 mHeight(height)
00357 {
00358 T * sum_table = new T[(width + 1) * (height + 1)];
00359 T * sqr_sum_table = new T[(width + 1) * (height + 1)];
00360 make_sum_tables(ar, sum_table, sqr_sum_table, width, height);
00361
00362 mRoot = initializeChild(ar, sum_table, sqr_sum_table, width, 0, 0, width, height, sqr(threshold));
00363
00364 delete[] sum_table;
00365 delete[] sqr_sum_table;
00366 }
00367
00368 template <typename T>
00369 AugmentedAreaSumTableQuadtree<T>::~AugmentedAreaSumTableQuadtree()
00370 {
00371 delete mRoot;
00372 mRoot = 0;
00373 }
00374
00375 template <typename T>
00376 unsigned int AugmentedAreaSumTableQuadtree<T>::getWidth() const
00377 {
00378 return mWidth;
00379 }
00380
00381 template <typename T>
00382 unsigned int AugmentedAreaSumTableQuadtree<T>::getHeight() const
00383 {
00384 return mHeight;
00385 }
00386
00387 #ifndef access
00388 template <typename T>
00389 const T& AugmentedAreaSumTableQuadtree<T>::access(const T ar[], unsigned int x, unsigned int y, unsigned int width)
00390 {
00391 return ar[y*width + x];
00392 }
00393
00394 template <typename T>
00395 T& AugmentedAreaSumTableQuadtree<T>::access(T ar[], unsigned int x, unsigned int y, unsigned int width)
00396 {
00397 return ar[y*width + x];
00398 }
00399 #endif
00400
00401 template <typename T>
00402 void AugmentedAreaSumTableQuadtree<T>::make_sum_tables(const T ar[], T sum_table[], T sqr_sum_table[], unsigned int width, unsigned int height)
00403 {
00404 unsigned int new_table_width = width + 1;
00405
00406 for(unsigned int i = 0; i < width + 1; i++)
00407 {
00408 access(sum_table, i, 0, new_table_width) = 0;
00409 access(sqr_sum_table, i, 0, new_table_width) = 0;
00410 }
00411
00412 for(unsigned int j = 1; j < height + 1; j++)
00413 {
00414 access(sum_table, 0, j, new_table_width) = 0;
00415 access(sqr_sum_table, 0, j, new_table_width) = 0;
00416 }
00417
00418 for(unsigned int i = 1; i < width + 1; i++)
00419 {
00420 for(unsigned int j = 1; j < height + 1; j++)
00421 {
00422 access(sum_table, i, j, new_table_width) =
00423 access(sum_table, i, j - 1, new_table_width)
00424 + access(sum_table, i - 1, j, new_table_width)
00425 - access(sum_table, i - 1, j - 1, new_table_width)
00426 + access(ar, i - 1, j - 1, width);
00427
00428 access(sqr_sum_table, i, j, new_table_width) =
00429 access(sqr_sum_table, i, j - 1, new_table_width)
00430 + access(sqr_sum_table, i - 1, j, new_table_width)
00431 - access(sqr_sum_table, i - 1, j - 1, new_table_width)
00432 + sqr(access(ar, i - 1, j - 1, width));
00433 }
00434 }
00435 }
00436
00437 template <typename T>
00438 T AugmentedAreaSumTableQuadtree<T>::getMean(const T sum_table[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1)
00439 {
00440 assert(x0 <= x1);
00441 assert(y0 <= y1);
00442
00443 T sum;
00444 unsigned int new_table_width = width + 1;
00445
00446 sum = access(sum_table, x1+1, y1+1, new_table_width) - access(sum_table, x1+1, y0, new_table_width)
00447 - access(sum_table, x0, y1+1, new_table_width) + access(sum_table, x0, y0, new_table_width);
00448
00449 return sum / ((x1 - x0 + 1) * (y1 - y0 + 1));
00450 }
00451
00452 template <typename T>
00453 T AugmentedAreaSumTableQuadtree<T>::getVariance(const T sqr_sum_table[], unsigned int width, T mean, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1)
00454 {
00455 assert(x0 <= x1);
00456 assert(y0 <= y1);
00457
00458 T sum;
00459 unsigned int new_table_width = width + 1;
00460
00461 sum = access(sqr_sum_table, x1+1, y1+1, new_table_width) - access(sqr_sum_table, x1+1, y0, new_table_width)
00462 - access(sqr_sum_table, x0, y1+1, new_table_width) + access(sqr_sum_table, x0, y0, new_table_width);
00463
00464 return sum / ((x1 - x0 + 1) * (y1 - y0 + 1)) - sqr(mean);
00465 }
00466
00467 template <typename T>
00468 const T& AugmentedAreaSumTableQuadtree<T>::operator()(unsigned int x, unsigned int y) const
00469 {
00470 return mRoot->getValue(x,y);
00471 }
00472
00473 template <typename T>
00474 unsigned int AugmentedAreaSumTableQuadtree<T>::getLevel(unsigned int x, unsigned int y) const
00475 {
00476 return mRoot->getLevel(x, y);
00477 }
00478
00479 template<typename T>
00480 typename AugmentedAreaSumTableQuadtree<T>::Node * AugmentedAreaSumTableQuadtree<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)
00481 {
00482 if (x0 + 1 == x1)
00483 {
00484 if (y0 + 1 == y1)
00485 {
00486 return new LeafNode(access(ar, x0, y0, width));
00487 }
00488 else
00489 {
00490 return new VerticalBranchNode(ar, sum_table, sqr_sum_table, width, x0, y0, x1, y1, threshold_sqr);
00491 }
00492 }
00493 else
00494 {
00495 if (y0 + 1 == y1)
00496 {
00497 return new HorizontalBranchNode(ar, sum_table, sqr_sum_table, width, x0, y0, x1, y1, threshold_sqr);
00498 }
00499 }
00500
00501 T mean(getMean(sum_table, width, x0, y0, x1 - 1, y1 - 1));
00502 T variance(getVariance(sqr_sum_table, width, mean, x0, y0, x1 - 1, y1 - 1));
00503
00504 if(variance <= threshold_sqr)
00505 {
00506 return new LeafNode(mean);
00507 }
00508 else
00509 {
00510 return new RectangularBranchNode(ar, sum_table, sqr_sum_table, width, x0, y0, x1, y1, threshold_sqr);
00511 }
00512
00513 assert(0);
00514 }
00515
00516 template<typename T>
00517 unsigned int AugmentedAreaSumTableQuadtree<T>::getNodeCount() const
00518 {
00519 return mRoot->getNodeCount();
00520 }
00521
00522 }
00523 }
00524
00525 #endif //AUGMENTED_AREA_SUM_TABLE_QUADTREE_H