Quadtree/src/headers/NaiveQuadtree.h

#ifndef NAIVE_QUADTREE_H
#define NAIVE_QUADTREE_H

#include <math.h>

#include "boost/concept/assert.hpp"

#include "Order2StatisticalConcept.h"
#include "QuadtreeConcept.h"

#include "dmath.h"
#include "debug.h"
#include "profile.h"

namespace za_co_codespot
{
        namespace datastructures
        {
                using math::sqr;

template <typename T>
class NaiveQuadtree
{
        //Make sure T is the correct type
        BOOST_CONCEPT_ASSERT((Order2StatisticalConcept<T>));

        //Make sure this implementation conforms to the compile time interface
        //This is for the developer 
#ifdef CODE_SPOT_CO_ZA_DEBUG
        BOOST_CONCEPT_ASSERT((QuadtreeConcept<NaiveQuadtree<T>, T>));
#endif

public:
        NaiveQuadtree(T ar[], unsigned int width, unsigned int height, T threshold);

        //destructor - releases nodes 
        virtual ~NaiveQuadtree();

        inline unsigned int getWidth() const;
        
        inline unsigned int getHeight() const;

        const T& operator()(unsigned int x, unsigned int y) const;

        unsigned int getLevel(unsigned int x, unsigned int y) const;

        unsigned int getNodeCount() const;

private:
        
        class Node
        {
        public:
                explicit Node();
                virtual ~Node() = 0;

                virtual inline const T& getValue(unsigned int x, unsigned int y) const = 0;

                virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const = 0;

                virtual unsigned int getNodeCount() const = 0;
        };

        template <unsigned int childCount>
        class BranchNode : public Node
        {
        public:
                explicit BranchNode();
                virtual ~BranchNode();  
                virtual inline const T& getValue(unsigned int x, unsigned int y) const = 0;
                virtual unsigned int getNodeCount() const;
        protected:
                Node * mChildren[childCount];
        };

        class RectangularBranchNode : public BranchNode<4>
        {
        public:
                RectangularBranchNode(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, T threshold_sqr);
                
                virtual inline const T& getValue(unsigned int x, unsigned int y) const;
                virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;

        private:
                unsigned int mMidX;
                unsigned int mMidY;
        };

        class HorizontalBranchNode : public BranchNode<2>
        {
        public:
                HorizontalBranchNode(T ar[], unsigned int width, unsigned int x0, unsigned int y0, int unsigned x1, int unsigned y1, T threshold_sqr);
                
                virtual inline const T& getValue(unsigned int x, unsigned int y) const;
                virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;
        private:
                
                unsigned int mMidX;
        };

        class VerticalBranchNode : public BranchNode<2>
        {
        public:
                VerticalBranchNode(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, T threshold_sqr);

                virtual inline const T& getValue(unsigned int x, unsigned int y) const;
                virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;
        private:
                unsigned int mMidY;
        };

        class LeafNode : public Node
        {
        public:
                LeafNode(T value);
                virtual inline const T& getValue(unsigned int x, unsigned int y) const;
                virtual unsigned int getNodeCount() const;
                virtual inline unsigned int getLevel(unsigned int x, unsigned int y) const;
        private:
                /* 
                        The data from the array that this leaf represents.
                        Only leaf nodes have data. The leaf need not represent a
                        single pixel. If more than one pixel is represented, this value
                        is the mean of all the pixels from the original data that this 
                        node represents.
                */
                T mValue;
        };


        static T getMean(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1);
        
        static T getVariance(T ar[], unsigned int width, T mean, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1);

        static Node * initializeChild(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, T threshold_sqr);

        Node * mRoot;

        unsigned int mWidth;
        
        unsigned int mHeight;
};

//-------------------------------------------------------------------------------------
// Implementation

template <typename T>
NaiveQuadtree<T>::Node::Node()
{
        //nothing to initialize;
}

template <typename T>
NaiveQuadtree<T>::Node::~Node()
{
        //nothing to release;
}

template <typename T>
template <unsigned int childCount>
NaiveQuadtree<T>::BranchNode<childCount>::~BranchNode()
{
        for(unsigned int k = 0; k < childCount; k++)
        {
                delete mChildren[k];
                mChildren[k] = 0;
        }
}

template <typename T>
template <unsigned int childCount>
NaiveQuadtree<T>::BranchNode<childCount>::BranchNode():
        Node()
{
        for(unsigned int k = 0; k < childCount; k++)
        {
                // Initilise to 0 and use asserts later to make 
                // sure we don't have initilised to proper values
                // or are not accessing the worng part of the array.
                mChildren[k] = 0;
        }
}

template <typename T>
template <unsigned int childCount>
unsigned int NaiveQuadtree<T>::BranchNode<childCount>::getNodeCount() const
{
        unsigned int sum = 1; // 1 for this node

        for(unsigned int k = 0; k < childCount; k++)
        {
                assert(mChildren[k] != 0); //See constructor

                sum += mChildren[k]->getNodeCount();
        }

        return sum;
}

template <typename T>
NaiveQuadtree<T>::RectangularBranchNode::RectangularBranchNode(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, T threshold_sqr):
        BranchNode<4>(),
        mMidX((x0 + x1) / 2),
        mMidY((y0 + y1) / 2)
{
        assert(x0 != x1 + 1);
        assert(y0 != y1 + 1);

        mChildren[0] = NaiveQuadtree<T>::initializeChild(ar, width, x0, y0, mMidX, mMidY, threshold_sqr);
        mChildren[1] = NaiveQuadtree<T>::initializeChild(ar, width, mMidX, y0, x1, mMidY, threshold_sqr);
        mChildren[2] = NaiveQuadtree<T>::initializeChild(ar, width, x0, mMidY, mMidX, y1, threshold_sqr);
        mChildren[3] = NaiveQuadtree<T>::initializeChild(ar, width, mMidX, mMidY, x1, y1, threshold_sqr);
        
}

template <typename T>
const T& NaiveQuadtree<T>::RectangularBranchNode::getValue(unsigned int x, unsigned int y) const
{
        if (x < mMidX)
        {
                if (y < mMidY)
                {
                        return mChildren[0]->getValue(x,y);
                }
                else
                {
                        return mChildren[2]->getValue(x,y);
                }
        }
        else
        {
                if (y < mMidY)
                {
                        return mChildren[1]->getValue(x,y);
                }
                else
                {
                        return mChildren[3]->getValue(x,y);
                }
        }
}

template <typename T>
unsigned int NaiveQuadtree<T>::RectangularBranchNode::getLevel(unsigned int x, unsigned int y) const
{
        if (x < mMidX)
        {
                if (y < mMidY)
                {
                        return mChildren[0]->getLevel(x,y) + 1;
                }
                else
                {
                        return mChildren[2]->getLevel(x,y) + 1;
                }
        }
        else
        {
                if (y < mMidY)
                {
                        return mChildren[1]->getLevel(x,y) + 1;
                }
                else
                {
                        return mChildren[3]->getLevel(x,y) + 1;
                }
        }
}

//-----------------------------
template <typename T>
NaiveQuadtree<T>::HorizontalBranchNode::HorizontalBranchNode(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, T threshold_sqr):
        BranchNode<2>(),
        mMidX((x0 + x1) / 2)
{
        assert(x0 != x1 + 1);
        assert(y0 != y1 + 1);

        mChildren[0] = NaiveQuadtree<T>::initializeChild(ar, width, x0, y0, mMidX, y1, threshold_sqr);
        mChildren[1] = NaiveQuadtree<T>::initializeChild(ar, width, mMidX, y0, x1, y1, threshold_sqr);  
}

template <typename T>
const T& NaiveQuadtree<T>::HorizontalBranchNode::getValue(unsigned int x, unsigned int y) const
{
        if (x < mMidX)
        {
                return mChildren[0]->getValue(x,y);
        }
        else
        {
                return mChildren[1]->getValue(x,y);             
        }
}

template <typename T>
unsigned int NaiveQuadtree<T>::HorizontalBranchNode::getLevel(unsigned int x, unsigned int y) const
{
        if (x < mMidX)
        {
                return mChildren[0]->getLevel(x,y) + 1;
        }
        else
        {
                return mChildren[1]->getLevel(x,y) + 1;         
        }
}
//------------------------
template <typename T>
NaiveQuadtree<T>::VerticalBranchNode::VerticalBranchNode(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, T threshold_sqr):
        BranchNode<2>(),
        mMidY((y0 + y1) / 2)
{
        assert(x0 != x1 + 1);
        assert(y0 != y1 + 1);

        mChildren[0] = NaiveQuadtree<T>::initializeChild(ar, width, x0, y0, x1, mMidY, threshold_sqr);
        mChildren[1] = NaiveQuadtree<T>::initializeChild(ar, width, x0, mMidY, x1, y1, threshold_sqr);  
}

template <typename T>
const T& NaiveQuadtree<T>::VerticalBranchNode::getValue(unsigned int x, unsigned int y) const
{
        if (y < mMidY)
        {
                return mChildren[0]->getValue(x,y);
        }
        else
        {
                return mChildren[1]->getValue(x,y);
        }
}

template <typename T>
unsigned int NaiveQuadtree<T>::VerticalBranchNode::getLevel(unsigned int x, unsigned int y) const
{
        if (y < mMidY)
        {
                return mChildren[0]->getLevel(x,y) + 1;
        }
        else
        {
                return mChildren[1]->getLevel(x,y) + 1;
        }
}
//--------------------

template <typename T>
NaiveQuadtree<T>::LeafNode::LeafNode(T value):
        Node(),
        mValue(value)
{}

template <typename T>
const T& NaiveQuadtree<T>::LeafNode::getValue(unsigned int x, unsigned int y) const
{
        return mValue;
}

template <typename T>
unsigned int NaiveQuadtree<T>::LeafNode::getLevel(unsigned int x, unsigned int y) const
{
        return 0;
}

template <typename T>
unsigned int NaiveQuadtree<T>::LeafNode::getNodeCount() const
{
        return 1;
}

template <typename T>
NaiveQuadtree<T>::NaiveQuadtree(T ar[], unsigned int width, unsigned int height, T threshold):
        mWidth(width),
        mHeight(height),
        mRoot(initializeChild(ar, width, 0, 0, width, height, sqr(threshold)))
{}

template <typename T>
NaiveQuadtree<T>::~NaiveQuadtree()
{
        delete mRoot;
        mRoot = 0;
}

template <typename T>
unsigned int NaiveQuadtree<T>::getWidth() const
{
        return mWidth;
}

template <typename T>
unsigned int NaiveQuadtree<T>::getHeight() const
{
        return mHeight;
}

template <typename T>
T NaiveQuadtree<T>::getMean(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1)
{
        T sum(0);

        assert(x0 != x1);
        assert(y0 != y1);

        for(unsigned int i = x0; i < x1; i++)
        {
                for(unsigned int j = y0; j < y1; j++)
                {
                        sum = sum + ar[j * width + i];
                }
        }

        return sum / ((x1 - x0) * (y1 - y0));
}

template <typename T>
T NaiveQuadtree<T>::getVariance(T ar[], unsigned int width, T mean, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1)
{
        T sum(0);

        assert(x0 != x1);
        assert(y0 != y1);

        for(unsigned int i = x0; i < x1; i++)
        {
                for(unsigned int j = y0; j < y1; j++)
                {
                        sum = sum + sqr(mean - ar[j * width + i]);
                }
        }
        
        return sum / ((x1 - x0) * (y1 - y0));
}

template <typename T>
const T& NaiveQuadtree<T>::operator()(unsigned int x, unsigned int y) const
{
        return mRoot->getValue(x,y);
}

template <typename T>
unsigned int NaiveQuadtree<T>::getLevel(unsigned int x, unsigned int y) const
{
        return mRoot->getLevel(x, y);
}

template<typename T>
typename NaiveQuadtree<T>::Node * NaiveQuadtree<T>::initializeChild(T ar[], unsigned int width, unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, T threshold_sqr)
{
        if (x0 + 1 == x1)//one pixel wide
        {
                if (y0 + 1 == y1)//one pixel high
                { 
                        return new LeafNode(ar[y0*width + x0]);
                }
                else
                {
                        return new VerticalBranchNode(ar, width, x0, y0, x1, y1, threshold_sqr);
                }
        }
        else
        {
                if (y0 + 1 == y1)//one pixel high
                { 
                        return new HorizontalBranchNode(ar, width, x0, y0, x1, y1, threshold_sqr);
                }
        }

        T mean(getMean(ar, width, x0, y0, x1, y1));
        T variance(getVariance(ar, width, mean, x0, y0, x1, y1));

        /*
                We use the variance in calculations, rather than the standard deviation so that 
                we need not call an expensive square root function here.

                Note that we use <= rather than < so that we can always have 
                a consistent behaviour when threshold is 0 for any (numeric) type.
        */
        if(variance <= threshold_sqr)
        {
                return new LeafNode(mean);
        }
        else
        {
                return new RectangularBranchNode(ar, width, x0, y0, x1, y1, threshold_sqr);
        }

        assert(UNREACHABLE);
}

template<typename T>
unsigned int NaiveQuadtree<T>::getNodeCount() const
{
        return mRoot->getNodeCount();
}

} //end namespace datastructures
} //end namespace za_co_codespot

#endif //NAIVE_QUADTREE_H