/* Redistribution and use in source and binary forms, with or
 * without modification, are permitted provided that the following
 * conditions are met:
 * 	Redistributions of source code must retain the above
 * 	copyright notice, this list of conditions and the following
 * 	disclaimer.
 * 	Redistributions in binary form must reproduce the above
 * 	copyright notice, this list of conditions and the following
 * 	disclaimer in the documentation and/or other materials
 * 	provided with the distribution.
 * 	The name of Contributor may not be used to endorse or
 * 	promote products derived from this software without
 * 	specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE CONTRIBUTORS BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
 * OF SUCH DAMAGE.
 * Copyright© 2009, Liu Liu All rights reserved.
 *
 * OpenCV functions for MSER extraction
 *
 * 1. there are two different implementation of MSER, one for grey image, one for color image
 * 2. the grey image algorithm is taken from: Linear Time Maximally Stable Extremal Regions;
 *    the paper claims to be faster than union-find method;
 *    it actually get 1.5~2m/s on my centrino L7200 1.2GHz laptop.
 * 3. the color image algorithm is taken from: Maximally Stable Colour Regions for Recognition and Match;
 *    it should be much slower than grey image method ( 3~4 times );
 *    the chi_table.h file is taken directly from paper's source code which is distributed under GPL.
 * 4. though the name is *contours*, the result actually is a list of point set.
 */

#include "precomp.hpp"
#include "opencv2/imgproc/imgproc_c.h"
#include <limits>

namespace cv
{

using std::vector;

class MSER_Impl : public MSER
{
public:
    struct Params
    {
        explicit Params( int _delta=5, double _maxVariation=0.25,
                        int _minArea=60, int _maxArea=14400,
                        double _minDiversity=.2, int _maxEvolution=200,
                        double _areaThreshold=1.01,
                        double _minMargin=0.003, int _edgeBlurSize=5 )
        {
            delta = _delta;
            minArea = _minArea;
            maxArea = _maxArea;
            maxVariation = _maxVariation;
            minDiversity = _minDiversity;
            pass2Only = false;
            maxEvolution = _maxEvolution;
            areaThreshold = _areaThreshold;
            minMargin = _minMargin;
            edgeBlurSize = _edgeBlurSize;
        }
        int delta;
        int minArea;
        int maxArea;
        double maxVariation;
        double minDiversity;
        bool pass2Only;

        int maxEvolution;
        double areaThreshold;
        double minMargin;
        int edgeBlurSize;
    };

    explicit MSER_Impl(const Params& _params) : params(_params) {}

    virtual ~MSER_Impl() {}

    void set(int propId, double value)
    {
        if( propId == DELTA )
            params.delta = cvRound(value);
        else if( propId == MIN_AREA )
            params.minArea = cvRound(value);
        else if( propId == MAX_AREA )
            params.maxArea = cvRound(value);
        else if( propId == PASS2_ONLY )
            params.pass2Only = value != 0;
        else
            CV_Error(CV_StsBadArg, "Unknown parameter id");
    }

    double get(int propId) const
    {
        double value = 0;
        if( propId == DELTA )
            value = params.delta;
        else if( propId == MIN_AREA )
            value = params.minArea;
        else if( propId == MAX_AREA )
            value = params.maxArea;
        else if( propId == PASS2_ONLY )
            value = params.pass2Only;
        else
            CV_Error(CV_StsBadArg, "Unknown parameter id");
        return value;
    }

    enum { DIR_SHIFT = 29, NEXT_MASK = ((1<<DIR_SHIFT)-1)  };

    struct Pixel
    {
        Pixel() : val(0) {}
        Pixel(int _val) : val(_val) {}

        int getGray(const Pixel* ptr0, const uchar* imgptr0, int mask) const
        {
            return imgptr0[this - ptr0] ^ mask;
        }
        int getNext() const { return (val & NEXT_MASK); }
        void setNext(int next) { val = (val & ~NEXT_MASK) | next; }

        int getDir() const { return (int)((unsigned)val >> DIR_SHIFT); }
        void setDir(int dir) { val = (val & NEXT_MASK) | (dir << DIR_SHIFT); }
        bool isVisited() const { return (val & ~NEXT_MASK) != 0; }

        int val;
    };
    typedef int PPixel;

    // the history of region grown
    struct CompHistory
    {
        CompHistory() { shortcut = child = 0; stable = val = size = 0; }
        CompHistory* shortcut;
        CompHistory* child;
        int stable; // when it ever stabled before, record the size
        int val;
        int size;
    };

    struct ConnectedComp
    {
        ConnectedComp()
        {
            init(0);
        }

        void init(int gray)
        {
            head = tail = 0;
            history = 0;
            size = 0;
            grey_level = gray;
            dvar = false;
            var = 0;
        }

        // add history chunk to a connected component
        void growHistory( CompHistory* h )
        {
            h->child = h;
            if( !history )
            {
                h->shortcut = h;
                h->stable = 0;
            }
            else
            {
                history->child = h;
                h->shortcut = history->shortcut;
                h->stable = history->stable;
            }
            h->val = grey_level;
            h->size = size;
            history = h;
        }

        // merging two connected components
        static void
        merge( const ConnectedComp* comp1,
              const ConnectedComp* comp2,
              ConnectedComp* comp,
              CompHistory* h,
              Pixel* pix0 )
        {
            comp->grey_level = comp2->grey_level;
            h->child = h;
            // select the winner by size
            if ( comp1->size < comp2->size )
                std::swap(comp1, comp2);

            if( !comp1->history )
            {
                h->shortcut = h;
                h->stable = 0;
            }
            else
            {
                comp1->history->child = h;
                h->shortcut = comp1->history->shortcut;
                h->stable = comp1->history->stable;
            }
            if( comp2->history && comp2->history->stable > h->stable )
                h->stable = comp2->history->stable;
            h->val = comp1->grey_level;
            h->size = comp1->size;
            // put comp1 to history
            comp->var = comp1->var;
            comp->dvar = comp1->dvar;
            if( comp1->size > 0 && comp2->size > 0 )
                pix0[comp1->tail].setNext(comp2->head);
            PPixel head = comp1->size > 0 ? comp1->head : comp2->head;
            PPixel tail = comp2->size > 0 ? comp2->tail : comp1->tail;
            // always made the newly added in the last of the pixel list (comp1 ... comp2)
            comp->head = head;
            comp->tail = tail;
            comp->history = h;
            comp->size = comp1->size + comp2->size;
        }

        float calcVariation( int delta ) const
        {
            if( !history )
                return 1.f;
            int val = grey_level;
            CompHistory* shortcut = history->shortcut;
            while( shortcut != shortcut->shortcut && shortcut->val + delta > val )
                shortcut = shortcut->shortcut;
            CompHistory* child = shortcut->child;
            while( child != child->child && child->val + delta <= val )
            {
                shortcut = child;
                child = child->child;
            }
            // get the position of history where the shortcut->val <= delta+val and shortcut->child->val >= delta+val
            history->shortcut = shortcut;
            return (float)(size - shortcut->size)/(float)shortcut->size;
            // here is a small modification of MSER where cal ||R_{i}-R_{i-delta}||/||R_{i-delta}||
            // in standard MSER, cal ||R_{i+delta}-R_{i-delta}||/||R_{i}||
            // my calculation is simpler and much easier to implement
        }

        bool isStable(const Params& p)
        {
            // tricky part: it actually check the stablity of one-step back
            if( !history || history->size <= p.minArea || history->size >= p.maxArea )
                return false;
            float div = (float)(history->size - history->stable)/(float)history->size;
            float _var = calcVariation( p.delta );
            bool _dvar = (var < _var) || (history->val + 1 < grey_level);
            bool stable = _dvar && !dvar && _var < p.maxVariation && div > p.minDiversity;
            var = _var;
            dvar = _dvar;
            if( stable )
                history->stable = history->size;
            return stable;
        }

        // convert the point set to CvSeq
        Rect label( Mat& labels, int lval, const Pixel* pix0, int step ) const
        {
            int* lptr = labels.ptr<int>();
            int lstep = labels.step/sizeof(lptr[0]);
            int xmin = INT_MAX, ymin = INT_MAX, xmax = INT_MIN, ymax = INT_MIN;

            for( PPixel pix = head; pix != 0; pix = pix0[pix].getNext() )
            {
                int y = pix/step;
                int x = pix - y*step;

                xmin = std::min(xmin, x);
                xmax = std::max(xmax, x);
                ymin = std::min(ymin, y);
                ymax = std::max(ymax, y);

                lptr[lstep*y + x] = lval;
            }

            return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1);
        }

        PPixel head;
        PPixel tail;
        CompHistory* history;
        int grey_level;
        int size;
        float var; // the current variation (most time is the variation of one-step back)
        bool dvar; // the derivative of last var
    };

    int detectAndLabel( InputArray _src, OutputArray _labels, OutputArray _bboxes );
    void detect( InputArray _src, vector<KeyPoint>& keypoints, InputArray _mask );

    void preprocess1( const Mat& img, int* level_size )
    {
        memset(level_size, 0, 256*sizeof(level_size[0]));

        int i, j, cols = img.cols, rows = img.rows;
        int step = cols;
        pixbuf.resize(step*rows);
        heapbuf.resize(cols*rows + 256);
        histbuf.resize(cols*rows);
        Pixel borderpix;
        borderpix.setDir(4);

        for( j = 0; j < step; j++ )
        {
            pixbuf[j] = pixbuf[j + (rows-1)*step] = borderpix;
        }

        for( i = 1; i < rows-1; i++ )
        {
            const uchar* imgptr = img.ptr(i);
            Pixel* pptr = &pixbuf[i*step];
            pptr[0] = pptr[cols-1] = borderpix;
            for( j = 1; j < cols-1; j++ )
            {
                int val = imgptr[j];
                level_size[val]++;
                pptr[j].val = 0;
            }
        }
    }

    void preprocess2( const Mat& img, int* level_size )
    {
        int i;

        for( i = 0; i < 128; i++ )
            std::swap(level_size[i], level_size[255-i]);

        if( !params.pass2Only )
        {
            int j, cols = img.cols, rows = img.rows;
            int step = cols;
            for( i = 1; i < rows-1; i++ )
            {
                Pixel* pptr = &pixbuf[i*step + 1];
                for( j = 1; j < cols-1; j++ )
                {
                    pptr[j].val = 0;
                }
            }
        }
    }

    void pass( const Mat& img, Mat& labels, int& lval, vector<Rect>& bboxvec,
              Size size, const int* level_size, int mask )
    {
        CompHistory* histptr = &histbuf[0];
        int step = size.width;
        Pixel *ptr0 = &pixbuf[0], *ptr = &ptr0[step+1];
        const uchar* imgptr0 = img.ptr();
        Pixel** heap[256];
        ConnectedComp comp[257];
        ConnectedComp* comptr = &comp[0];

        heap[0] = &heapbuf[0];
        heap[0][0] = 0;

        for( int i = 1; i < 256; i++ )
        {
            heap[i] = heap[i-1] + level_size[i-1] + 1;
            heap[i][0] = 0;
        }

        comptr->grey_level = 256;
        comptr++;
        comptr->grey_level = ptr->getGray(ptr0, imgptr0, mask);
        ptr->setDir(1);
        int dir[] = { 0, 1, step, -1, -step };
        for( ;; )
        {
            int curr_gray = ptr->getGray(ptr0, imgptr0, mask);
            int nbr_idx = ptr->getDir();
            // take tour of all the 4 directions
            for( ; nbr_idx <= 4; nbr_idx++ )
            {
                // get the neighbor
                Pixel* ptr_nbr = ptr + dir[nbr_idx];
                if( !ptr_nbr->isVisited() )
                {
                    // set dir=1, next=0
                    ptr_nbr->val = 1 << DIR_SHIFT;
                    int nbr_gray = ptr_nbr->getGray(ptr0, imgptr0, mask);
                    if( nbr_gray < curr_gray )
                    {
                        // when the value of neighbor smaller than current
                        // push current to boundary heap and make the neighbor to be the current one
                        // create an empty comp
                        *(++heap[curr_gray]) = ptr;
                        ptr->val = (nbr_idx+1) << DIR_SHIFT;
                        ptr = ptr_nbr;
                        comptr++;
                        comptr->init(nbr_gray);
                        curr_gray = nbr_gray;
                        nbr_idx = 0;
                        continue;
                    }
                    // otherwise, push the neighbor to boundary heap
                    *(++heap[nbr_gray]) = ptr_nbr;
                }
            }

            // set dir = nbr_idx, next = 0
            ptr->val = nbr_idx << DIR_SHIFT;
            int ptrofs = (int)(ptr - ptr0);
            CV_Assert(ptrofs != 0);

            // add a pixel to the pixel list
            if( comptr->tail )
                ptr0[comptr->tail].setNext(ptrofs);
            else
                comptr->head = ptrofs;
            comptr->tail = ptrofs;
            comptr->size++;
            // get the next pixel from boundary heap
            if( *heap[curr_gray] )
            {
                ptr = *heap[curr_gray];
                heap[curr_gray]--;
            }
            else
            {
                for( curr_gray++; curr_gray < 256; curr_gray++ )
                {
                    if( *heap[curr_gray] )
                        break;
                }
                if( curr_gray >= 256 )
                    break;

                ptr = *heap[curr_gray];
                heap[curr_gray]--;
                if( curr_gray < comptr[-1].grey_level )
                {
                    // check the stablity and push a new history, increase the grey level
                    if( comptr->isStable(params) )
                    {
                        Rect box = comptr->label( labels, lval++, ptr0, step );
                        bboxvec.push_back(box);
                    }
                    comptr->growHistory( histptr++ );
                    comptr[0].grey_level = curr_gray;
                }
                else
                {
                    // keep merging top two comp in stack until the grey level >= pixel_val
                    for(;;)
                    {
                        comptr--;
                        ConnectedComp::merge(comptr+1, comptr, comptr, histptr++, ptr0);
                        if( curr_gray <= comptr[0].grey_level )
                            break;
                        if( curr_gray < comptr[-1].grey_level )
                        {
                            // check the stablity here otherwise it wouldn't be an ER
                            if( comptr->isStable(params) )
                            {
                                Rect box = comptr->label( labels, lval++, ptr0, step );
                                bboxvec.push_back(box);
                            }
                            comptr->growHistory( histptr++ );
                            comptr[0].grey_level = curr_gray;
                            break;
                        }
                    }
                }
            }
        }
    }

    Mat tempsrc;
    vector<Pixel> pixbuf;
    vector<Pixel*> heapbuf;
    vector<CompHistory> histbuf;

    Params params;
};

/*

TODO:
the color MSER has not been completely refactored yet. We leave it mostly as-is,
with just enough changes to convert C structures to C++ ones and
add support for color images into MSER_Impl::detectAndLabel.
*/

const int TABLE_SIZE = 400;

static const float chitab3[]=
{
    0.f,  0.0150057f,  0.0239478f,  0.0315227f,
    0.0383427f,  0.0446605f,  0.0506115f,  0.0562786f,
    0.0617174f,  0.0669672f,  0.0720573f,  0.0770099f,
    0.081843f,  0.0865705f,  0.0912043f,  0.0957541f,
    0.100228f,  0.104633f,  0.108976f,  0.113261f,
    0.117493f,  0.121676f,  0.125814f,  0.12991f,
    0.133967f,  0.137987f,  0.141974f,  0.145929f,
    0.149853f,  0.15375f,  0.15762f,  0.161466f,
    0.165287f,  0.169087f,  0.172866f,  0.176625f,
    0.180365f,  0.184088f,  0.187794f,  0.191483f,
    0.195158f,  0.198819f,  0.202466f,  0.2061f,
    0.209722f,  0.213332f,  0.216932f,  0.220521f,
    0.2241f,  0.22767f,  0.231231f,  0.234783f,
    0.238328f,  0.241865f,  0.245395f,  0.248918f,
    0.252435f,  0.255947f,  0.259452f,  0.262952f,
    0.266448f,  0.269939f,  0.273425f,  0.276908f,
    0.280386f,  0.283862f,  0.287334f,  0.290803f,
    0.29427f,  0.297734f,  0.301197f,  0.304657f,
    0.308115f,  0.311573f,  0.315028f,  0.318483f,
    0.321937f,  0.32539f,  0.328843f,  0.332296f,
    0.335749f,  0.339201f,  0.342654f,  0.346108f,
    0.349562f,  0.353017f,  0.356473f,  0.35993f,
    0.363389f,  0.366849f,  0.37031f,  0.373774f,
    0.377239f,  0.380706f,  0.384176f,  0.387648f,
    0.391123f,  0.3946f,  0.39808f,  0.401563f,
    0.405049f,  0.408539f,  0.412032f,  0.415528f,
    0.419028f,  0.422531f,  0.426039f,  0.429551f,
    0.433066f,  0.436586f,  0.440111f,  0.44364f,
    0.447173f,  0.450712f,  0.454255f,  0.457803f,
    0.461356f,  0.464915f,  0.468479f,  0.472049f,
    0.475624f,  0.479205f,  0.482792f,  0.486384f,
    0.489983f,  0.493588f,  0.4972f,  0.500818f,
    0.504442f,  0.508073f,  0.511711f,  0.515356f,
    0.519008f,  0.522667f,  0.526334f,  0.530008f,
    0.533689f,  0.537378f,  0.541075f,  0.54478f,
    0.548492f,  0.552213f,  0.555942f,  0.55968f,
    0.563425f,  0.56718f,  0.570943f,  0.574715f,
    0.578497f,  0.582287f,  0.586086f,  0.589895f,
    0.593713f,  0.597541f,  0.601379f,  0.605227f,
    0.609084f,  0.612952f,  0.61683f,  0.620718f,
    0.624617f,  0.628526f,  0.632447f,  0.636378f,
    0.64032f,  0.644274f,  0.648239f,  0.652215f,
    0.656203f,  0.660203f,  0.664215f,  0.668238f,
    0.672274f,  0.676323f,  0.680384f,  0.684457f,
    0.688543f,  0.692643f,  0.696755f,  0.700881f,
    0.70502f,  0.709172f,  0.713339f,  0.717519f,
    0.721714f,  0.725922f,  0.730145f,  0.734383f,
    0.738636f,  0.742903f,  0.747185f,  0.751483f,
    0.755796f,  0.760125f,  0.76447f,  0.768831f,
    0.773208f,  0.777601f,  0.782011f,  0.786438f,
    0.790882f,  0.795343f,  0.799821f,  0.804318f,
    0.808831f,  0.813363f,  0.817913f,  0.822482f,
    0.827069f,  0.831676f,  0.836301f,  0.840946f,
    0.84561f,  0.850295f,  0.854999f,  0.859724f,
    0.864469f,  0.869235f,  0.874022f,  0.878831f,
    0.883661f,  0.888513f,  0.893387f,  0.898284f,
    0.903204f,  0.908146f,  0.913112f,  0.918101f,
    0.923114f,  0.928152f,  0.933214f,  0.938301f,
    0.943413f,  0.94855f,  0.953713f,  0.958903f,
    0.964119f,  0.969361f,  0.974631f,  0.979929f,
    0.985254f,  0.990608f,  0.99599f,  1.0014f,
    1.00684f,  1.01231f,  1.01781f,  1.02335f,
    1.02891f,  1.0345f,  1.04013f,  1.04579f,
    1.05148f,  1.05721f,  1.06296f,  1.06876f,
    1.07459f,  1.08045f,  1.08635f,  1.09228f,
    1.09826f,  1.10427f,  1.11032f,  1.1164f,
    1.12253f,  1.1287f,  1.1349f,  1.14115f,
    1.14744f,  1.15377f,  1.16015f,  1.16656f,
    1.17303f,  1.17954f,  1.18609f,  1.19269f,
    1.19934f,  1.20603f,  1.21278f,  1.21958f,
    1.22642f,  1.23332f,  1.24027f,  1.24727f,
    1.25433f,  1.26144f,  1.26861f,  1.27584f,
    1.28312f,  1.29047f,  1.29787f,  1.30534f,
    1.31287f,  1.32046f,  1.32812f,  1.33585f,
    1.34364f,  1.3515f,  1.35943f,  1.36744f,
    1.37551f,  1.38367f,  1.39189f,  1.4002f,
    1.40859f,  1.41705f,  1.42561f,  1.43424f,
    1.44296f,  1.45177f,  1.46068f,  1.46967f,
    1.47876f,  1.48795f,  1.49723f,  1.50662f,
    1.51611f,  1.52571f,  1.53541f,  1.54523f,
    1.55517f,  1.56522f,  1.57539f,  1.58568f,
    1.59611f,  1.60666f,  1.61735f,  1.62817f,
    1.63914f,  1.65025f,  1.66152f,  1.67293f,
    1.68451f,  1.69625f,  1.70815f,  1.72023f,
    1.73249f,  1.74494f,  1.75757f,  1.77041f,
    1.78344f,  1.79669f,  1.81016f,  1.82385f,
    1.83777f,  1.85194f,  1.86635f,  1.88103f,
    1.89598f,  1.91121f,  1.92674f,  1.94257f,
    1.95871f,  1.97519f,  1.99201f,  2.0092f,
    2.02676f,  2.04471f,  2.06309f,  2.08189f,
    2.10115f,  2.12089f,  2.14114f,  2.16192f,
    2.18326f,  2.2052f,  2.22777f,  2.25101f,
    2.27496f,  2.29966f,  2.32518f,  2.35156f,
    2.37886f,  2.40717f,  2.43655f,  2.46709f,
    2.49889f,  2.53206f,  2.56673f,  2.60305f,
    2.64117f,  2.6813f,  2.72367f,  2.76854f,
    2.81623f,  2.86714f,  2.92173f,  2.98059f,
    3.04446f,  3.1143f,  3.19135f,  3.27731f,
    3.37455f,  3.48653f,  3.61862f,  3.77982f,
    3.98692f,  4.2776f,  4.77167f,  133.333f
};


struct MSCRNode;

struct TempMSCR
{
    MSCRNode* head;
    MSCRNode* tail;
    double m; // the margin used to prune area later
    int size;
};

struct MSCRNode
{
    // the stable mscr should be:
    // bigger than minArea and smaller than maxArea
    // differ from its ancestor more than minDiversity
    bool isStable( const MSER_Impl::Params& params ) const
    {
        if( size <= params.minArea || size >= params.maxArea )
            return 0;
        if( gmsr == NULL )
            return 1;
        double div = (double)(size - gmsr->size)/(double)size;
        return div > params.minDiversity;
    }

    void init( int _index )
    {
        gmsr = tmsr = NULL;
        reinit = 0xffff;
        rank = 0;
        sizei = size = 1;
        prev = next = shortcut = this;
        index = _index;
    }

    // to find the root of one region
    MSCRNode* findRoot()
    {
        MSCRNode* x = this;
        MSCRNode* _prev = x;
        MSCRNode* _next;
        for(;;)
        {
            _next = x->shortcut;
            x->shortcut = _prev;
            if( _next == x )
                break;
            _prev = x;
            x = _next;
        }
        MSCRNode* root = x;
        for(;;)
        {
            _prev = x->shortcut;
            x->shortcut = root;
            if( _prev == x )
                break;
            x = _prev;
        }
        return root;
    }

    MSCRNode* shortcut;
    // to make the finding of root less painful
    MSCRNode* prev;
    MSCRNode* next;
    // a point double-linked list
    TempMSCR* tmsr;
    // the temporary msr (set to NULL at every re-initialise)
    TempMSCR* gmsr;
    // the global msr (once set, never to NULL)
    int index;
    // the index of the node, at this point, it should be x at the first 16-bits, and y at the last 16-bits.
    int rank;
    int reinit;
    int size, sizei;
    double dt, di;
    double s;
};

struct MSCREdge
{
    double init(float _chi, MSCRNode* _left, MSCRNode* _right)
    {
        chi = _chi;
        left = _left;
        right = _right;
        return chi;
    }
    float chi;
    MSCRNode* left;
    MSCRNode* right;
};

static float ChiSquaredDistance( const uchar* x, const uchar* y )
{
    return (float)((x[0]-y[0])*(x[0]-y[0]))/(float)(x[0]+y[0]+FLT_EPSILON)+
           (float)((x[1]-y[1])*(x[1]-y[1]))/(float)(x[1]+y[1]+FLT_EPSILON)+
           (float)((x[2]-y[2])*(x[2]-y[2]))/(float)(x[2]+y[2]+FLT_EPSILON);
}


// the preprocess to get the edge list with proper gaussian blur
static int preprocessMSER_8UC3( MSCRNode* node, MSCREdge* edge,
                                double& total, const Mat& src,
                                Mat& dx, Mat& dy, int Ne, int edgeBlurSize )
{
    int nch = src.channels();
    int i, j, nrows = src.rows, ncols = src.cols;
    float* dxptr = 0;
    float* dyptr = 0;
    for( i = 0; i < nrows; i++ )
    {
        const uchar* srcptr = src.ptr(i);
        const uchar* nextsrc = src.ptr(std::min(i+1, nrows-1));
        dxptr = dx.ptr<float>(i);
        dyptr = dy.ptr<float>(i);

        for( j = 0; j < ncols-1; j++ )
        {
            dxptr[j] = ChiSquaredDistance( srcptr + j*nch, srcptr + (j+1)*nch );
            dyptr[j] = ChiSquaredDistance( srcptr + j*nch, nextsrc + j*nch );
        }
        dyptr[ncols-1] = ChiSquaredDistance( srcptr + (ncols-1)*nch, nextsrc + (ncols-1)*nch );
    }

    // get dx and dy and blur it
    if( edgeBlurSize >= 1 )
    {
        GaussianBlur(dx, dx, Size(edgeBlurSize, edgeBlurSize), 0);
        GaussianBlur(dy, dy, Size(edgeBlurSize, edgeBlurSize), 0);
    }
    dxptr = dx.ptr<float>();
    dyptr = dy.ptr<float>();
    // assian dx, dy to proper edge list and initialize mscr node
    // the nasty code here intended to avoid extra loops
    MSCRNode* nodeptr = node;
    for( j = 0; j < ncols-1; j++ )
    {
        nodeptr[j].init(j);
        total += edge[j].init(dxptr[j], nodeptr+j, nodeptr+j+1);
    }
    dxptr += ncols - 1;
    edge += ncols - 1;
    nodeptr[ncols-1].init(ncols - 1);
    nodeptr += ncols;
    for( i = 1; i < nrows; i++ )
    {
        for( j = 0; j < ncols-1; j++ )
        {
            nodeptr[j].init( (i<<16)|j );
            total += edge[j*2].init(dyptr[j], nodeptr + j - ncols, nodeptr + j);
            total += edge[j*2+1].init(dxptr[j], nodeptr + j, nodeptr + j + 1);
        }
        nodeptr[ncols-1].init((i<<16)|(ncols - 1));
        total += edge[(ncols-1)*2].init(dyptr[ncols-1], nodeptr - 1, nodeptr + ncols-1);
        dxptr += ncols-1;
        dyptr += ncols;
        edge += 2*ncols - 1;
        nodeptr += ncols;
    }

    return Ne;
}

class LessThanEdge
{
public:
    bool operator()(const MSCREdge& a, const MSCREdge& b) const { return a.chi < b.chi; }
};


static void
extractMSER_8uC3( const Mat& src, Mat& labels,
                  vector<Rect>& bboxvec,
                  const MSER_Impl::Params& params )
{
    int npixels = src.cols*src.rows;
    int currlabel = 0;
    int* lptr = labels.ptr<int>();
    int lstep = (int)(labels.step/sizeof(int));

    vector<MSCRNode> mapvec(npixels);
    MSCRNode* map = &mapvec[0];
    int Ne = npixels*2 - src.cols - src.rows;
    vector<MSCREdge> edgevec(Ne+1);
    MSCREdge* edge = &edgevec[0];
    vector<TempMSCR> mscrvec(npixels);
    TempMSCR* mscr = &mscrvec[0];
    double emean = 0;
    Mat dx( src.rows, src.cols-1, CV_32FC1 );
    Mat dy( src.rows, src.cols, CV_32FC1 );
    Ne = preprocessMSER_8UC3( map, edge, emean, src, dx, dy, Ne, params.edgeBlurSize );
    emean = emean / (double)Ne;
    std::sort(edge, edge + Ne, LessThanEdge());
    MSCREdge* edge_ub = edge+Ne;
    MSCREdge* edgeptr = edge;
    TempMSCR* mscrptr = mscr;

    // the evolution process
    for ( int i = 0; i < params.maxEvolution; i++ )
    {
        double k = (double)i/(double)params.maxEvolution*(TABLE_SIZE-1);
        int ti = cvFloor(k);
        double reminder = k-ti;
        double thres = emean*(chitab3[ti]*(1-reminder)+chitab3[ti+1]*reminder);
        // to process all the edges in the list that chi < thres
        while( edgeptr < edge_ub && edgeptr->chi < thres )
        {
            MSCRNode* lr = edgeptr->left->findRoot();
            MSCRNode* rr = edgeptr->right->findRoot();
            // get the region root (who is responsible)
            if ( lr != rr )
            {
                MSCRNode* tmp;
                // rank idea take from: N-tree Disjoint-Set Forests for Maximally Stable Extremal Regions
                if ( rr->rank > lr->rank )
                {
                    CV_SWAP( lr, rr, tmp );
                }
                else if ( lr->rank == rr->rank )
                {
                    // at the same rank, we will compare the size
                    if( lr->size > rr->size )
                    {
                        CV_SWAP( lr, rr, tmp );
                    }
                    lr->rank++;
                }
                rr->shortcut = lr;
                lr->size += rr->size;
                // join rr to the end of list lr (lr is a endless double-linked list)
                lr->prev->next = rr;
                lr->prev = rr->prev;
                rr->prev->next = lr;
                rr->prev = lr;
                // area threshold force to reinitialize
                if ( lr->size > (lr->size-rr->size)*params.areaThreshold )
                {
                    lr->sizei = lr->size;
                    lr->reinit = i;
                    if ( lr->tmsr != NULL )
                    {
                        lr->tmsr->m = lr->dt-lr->di;
                        lr->tmsr = NULL;
                    }
                    lr->di = edgeptr->chi;
                    lr->s = 1e10;
                }
                lr->dt = edgeptr->chi;
                if ( i > lr->reinit )
                {
                    double s = (double)(lr->size-lr->sizei)/(lr->dt-lr->di);
                    if ( s < lr->s )
                    {
                        // skip the first one and check stablity
                        if ( i > lr->reinit+1 && lr->isStable( params ) )
                        {
                            if ( lr->tmsr == NULL )
                            {
                                lr->gmsr = lr->tmsr = mscrptr;
                                mscrptr++;
                            }
                            lr->tmsr->size = lr->size;
                            lr->tmsr->head = lr;
                            lr->tmsr->tail = lr->prev;
                            lr->tmsr->m = 0;
                        }
                        lr->s = s;
                    }
                }
            }
            edgeptr++;
        }
        if ( edgeptr >= edge_ub )
            break;
    }

    for( TempMSCR* ptr = mscr; ptr < mscrptr; ptr++ )
    {
        // to prune area with margin less than minMargin
        if( ptr->m > params.minMargin )
        {
            int xmin = INT_MAX, ymin = INT_MAX, xmax = INT_MIN, ymax = INT_MIN;
            currlabel++;
            MSCRNode* lpt = ptr->head;
            for( int i = 0; i < ptr->size; i++ )
            {
                int x = (lpt->index)&0xffff;
                int y = (lpt->index)>>16;
                lpt = lpt->next;

                xmin = std::min(xmin, x);
                xmax = std::max(xmax, x);
                ymin = std::min(ymin, y);
                ymax = std::max(ymax, y);

                lptr[lstep*y + x] = currlabel;
            }
            bboxvec.push_back(Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1));
        }
    }
}

int MSER_Impl::detectAndLabel( InputArray _src, OutputArray _labels, OutputArray _bboxes )
{
    Mat src = _src.getMat();
    size_t npix = src.total();
    vector<Rect> bboxvec;

    if( npix == 0 )
    {
        _labels.release();
        return 0;
    }

    Size size = src.size();
    _labels.create( size, CV_32S );
    Mat labels = _labels.getMat();
    labels.setTo(Scalar::all(0));

    if( src.type() == CV_8U )
    {
        int level_size[256];
        int lval = 1;

        if( !src.isContinuous() )
        {
            src.copyTo(tempsrc);
            src = tempsrc;
        }

        // darker to brighter (MSER+)
        preprocess1( src, level_size );
        if( !params.pass2Only )
            pass( src, labels, lval, bboxvec, size, level_size, 0 );
        // brighter to darker (MSER-)
        preprocess2( src, level_size );
        pass( src, labels, lval, bboxvec, size, level_size, 255 );
    }
    else
    {
        CV_Assert( src.type() == CV_8UC3 || src.type() == CV_8UC4 );
        extractMSER_8uC3( src, labels, bboxvec, params );
    }
    if( _bboxes.needed() )
        Mat(bboxvec).copyTo(_bboxes);
    return (int)bboxvec.size();
}

void MSER_Impl::detect( InputArray _image, vector<KeyPoint>& keypoints, InputArray _mask )
{
    vector<Rect> bboxes;
    vector<Point> reg;
    Mat labels, mask = _mask.getMat();

    int i, x, y, ncomps = detectAndLabel(_image, labels, bboxes);
    CV_Assert( ncomps == (int)bboxes.size() );

    keypoints.clear();
    for( i = 0; i < ncomps; i++ )
    {
        Rect r = bboxes[i];
        reg.reserve(r.area());
        reg.clear();

        for( y = r.y; y < r.y + r.height; y++ )
        {
            const int* lptr = labels.ptr<int>(y);
            for( x = r.x; x < r.x + r.width; x++ )
            {
                if( lptr[x] == i+1 )
                    reg.push_back(Point(x, y));
            }
        }
        // TODO check transformation from MSER region to KeyPoint
        RotatedRect rect = fitEllipse(Mat(reg));
        float diam = std::sqrt(rect.size.height*rect.size.width);

        if( diam > std::numeric_limits<float>::epsilon() && r.contains(rect.center) &&
            (mask.empty() || mask.at<uchar>(cvRound(rect.center.y), cvRound(rect.center.x)) != 0) )
            keypoints.push_back( KeyPoint(rect.center, diam) );
    }
}

Ptr<MSER> MSER::create( int _delta, int _min_area, int _max_area,
      double _max_variation, double _min_diversity,
      int _max_evolution, double _area_threshold,
      double _min_margin, int _edge_blur_size )
{
    return makePtr<MSER_Impl>(
        MSER_Impl::Params(_delta, _min_area, _max_area,
                          _max_variation, _min_diversity,
                          _max_evolution, _area_threshold,
                          _min_margin, _edge_blur_size));
}

}