/************************************************************************
            Blob.cpp

FUNCTIONALITY: Implementation of the Blob class and some helper classes to perform
               some calculations on it
AUTHOR: Inspecta S.L.
MODIFICATIONS (Modification, Author, Date):

**************************************************************************/


#include "blob.h"
#include <list>
#include "opencv2/core/types_c.h"
#include "opencv2/imgproc/imgproc_c.h"

Blob::Blob()
{
    m_area = m_perimeter = -1;
    m_externalPerimeter = m_meanGray = m_stdDevGray = -1;
    m_boundingBox.width = -1;
    m_ellipse.size.width = -1;
    m_id = -1;
    m_toBeDeleted = 0;
    m_deleteRequestOwnerBlob = NULL;
    m_isJoined = false;
    m_startPassed = false;
}

Blob::Blob(LabelID id, const cv::Point& startPoint, const cv::Size& originalImageSize) : m_externalContour(startPoint, originalImageSize)
{
    m_externalContour.m_parent = this;
    m_id = id;
    m_area = m_perimeter = -1;
    m_externalPerimeter = m_meanGray = m_stdDevGray = -1;
    m_boundingBox.width = -1;
    m_ellipse.size.width = -1;
    m_originalImageSize = originalImageSize;
    m_toBeDeleted = 0;
    m_deleteRequestOwnerBlob = NULL;
    m_isJoined = false;
    m_startPassed = false;
}

//! Copy constructor
Blob::Blob(const Blob& src)
{
    *this = src;
}

Blob::Blob(const Blob *src)
{
    if (src != NULL) {
        *this = *src;
    }
}

Blob& Blob::operator=(const Blob& src)
{
    if (this != &src) {
        m_id = src.m_id;
        m_area = src.m_area;
        m_perimeter = src.m_perimeter;
        m_externalPerimeter = src.m_externalPerimeter;
        m_meanGray = src.m_meanGray;
        m_stdDevGray = src.m_stdDevGray;
        m_boundingBox = src.m_boundingBox;
        m_ellipse = src.m_ellipse;
        m_originalImageSize = src.m_originalImageSize;
        m_toBeDeleted = src.m_toBeDeleted;
        m_deleteRequestOwnerBlob = src.m_deleteRequestOwnerBlob;
        m_startPassed = false;
        //! clear all current blob contours
        clearContours();
        m_externalContour = src.m_externalContour;

        //! copy all internal contours
        if (src.m_internalContours.size() != 0) {
            //m_internalContours = t_contourList(src.m_internalContours.size());
            BlobContourList::const_iterator itSrc,enSrc;
            BlobContourList::iterator it;

            itSrc = src.m_internalContours.begin();
            enSrc = src.m_internalContours.end();

            while (itSrc != enSrc) {
                m_internalContours.push_back(new BlobContour(*itSrc));
                itSrc++;
            }
        }
    }

    m_isJoined = src.m_isJoined;
    m_joinedBlobs.clear();
    if (m_isJoined) {
        std::list<Blob*>::const_iterator it,en = src.m_joinedBlobs.end();
        for (it = src.m_joinedBlobs.begin(); it != en; ++it) {
            m_joinedBlobs.push_back(new Blob(*it));
        }
    }

    return *this;
}

Blob::~Blob()
{
    if (m_isJoined) {
        std::list<Blob*>::iterator it,en = m_joinedBlobs.end();
        for (it = m_joinedBlobs.begin(); it != en; ++it) {
            delete (*it);
        }
    }

    clearContours();
}

void Blob::clearContours()
{	
    BlobContourList::iterator it = m_internalContours.begin(), en = m_internalContours.end();
    for (; it != en; ++it) {
        delete (*it);
    }
    m_internalContours.clear();
    m_externalContour.reset();
}

void Blob::addInternalContour(const BlobContour &newContour)
{
    m_internalContours.push_back(new BlobContour(newContour));
}

//! Shows if the blob has associated information
bool Blob::isEmpty()
{
    return getExternalContour()->m_contour.size() == 0;
}

/**
- FUNCTIONS: Area
- FUNCTIONALITY: Get blob area, ie. external contour area minus internal contours area
- PARAMETERS: -
- RESULT: -
- RESTRICTIONS: -
- AUTHOR: trees
- DATE OF CREATION: 30/04/2008
- MODIFICATION: Date. Author. Description.
*/
double Blob::area(AreaMode areaCompMode)
{
    double area=0;
    switch (areaCompMode)
    {
    case GREEN:
    {
        BlobContourList::iterator itContour;
        if (m_isJoined) {
            std::list<Blob *>::iterator it, en = m_joinedBlobs.end();
            for (it = m_joinedBlobs.begin(); it != en; ++it) {
                area += (*it)->area();
            }
        }
        else {
            area = m_externalContour.getArea();
            itContour = m_internalContours.begin();
            while (itContour != m_internalContours.end()) {
                if (*itContour) {
                    area -= (*itContour)->getArea();
                }
                itContour++;
            }
        }
        break;
    }
    case PIXELWISE:
    {
        cv::Rect bbox = getBoundingBox();
        cv::Mat image = cv::Mat::zeros(bbox.height, bbox.width, CV_8UC1);
        fillBlob(image, cv::Scalar(255), -bbox.x, -bbox.y, true);
        area = countNonZero(image);
        break;
    }
    }

    return area;
}

/**
- FUNCTIONS: Perimeter
- FUNCTIONALITY: Get perimeter blob, ie. sum of the length of all the contours
- PARAMETERS: -
- RESULT: -
- RESTRICTIONS: -
- AUTHOR: trees
- DATE OF CREATION: 30/04/2008
- MODIFICATION: Date. Author. Description.
*/
double Blob::perimeter()
{
    double perimeter = 0;
    BlobContourList::iterator itContour;

    if (m_isJoined) {
        std::list<Blob *>::iterator it, en = m_joinedBlobs.end();
        for (it = m_joinedBlobs.begin(); it != en; ++it) {
            perimeter += (*it)->perimeter();
        }
    }
    else {
        perimeter = m_externalContour.getPerimeter();
        itContour = m_internalContours.begin();
        while (itContour != m_internalContours.end()) {
            if (*itContour) {
                perimeter += (*itContour)->getPerimeter();
            }
            itContour++;
        }
    }

    return perimeter;
}

/**
- FUNCTIONS: exterior
- FUNCIONALITAT: Return true for external blobs
- PARAMETERS:
    - xBorder: true to consider blobs touching horizontal borders as external
    - yBorder: true to consider blobs touching vertical borders as external
- RESULTAT:
    -
- RESTRICCIONS:
    -
- AUTOR: rborras
- DATE OF CREATION: 2008/05/06
- MODIFICATION: Data. Autor. Description.
*/
int	Blob::exterior(IplImage *mask, bool xBorder /* = true */, bool yBorder /* = true */)
{
    int result = 0;
    if (externalPerimeter(mask, xBorder, yBorder) > 0)
    {
        result = 1;
    }

    return result;
}

int	Blob::exterior(const cv::Mat& mask, bool xBorder /* = true */, bool yBorder /* = true */)
{
    IplImage temp = cvIplImage(mask);
    return exterior(&temp, xBorder, yBorder);
}

/**
- FUNCTIONS: externalPerimeter
- FUNCIONALITAT: Get external perimeter (perimeter touching image borders)
- PARAMETERS:
    - maskImage: if != NULL, counts maskImage black pixels as external pixels and contour points touching them are counted as external contour points.
    - xBorder: true to consider blobs touching horizontal borders as external
    - yBorder: true to consider blobs touching vertical borders as external
- RESULTAT:
    -
- RESTRICCIONS:
    -
- AUTOR: rborras
- DATE OF CREATION: 2008/05/05
- MODIFICATION: Data. Autor. Description.
- NOTA: If BlobContour::getContourPoints aproximates contours with a method different that NONE, this function will not give correct results
*/
double Blob::externalPerimeter(IplImage *maskImage, bool xBorder /* = true */, bool yBorder /* = true */)
{
    PointList externalContour, externalPoints;
    cv::Point actualPoint, previousPoint;
    bool find = false;
    //int i,j;
    int delta = 0;

    //! it is calculated?
    /*if (m_externalPerimeter != -1)
    {
        return m_externalPerimeter;
    }*/


    if (m_isJoined) {
        //! it an en are always different at first assignment (if m_isJoined is true I have at least one joined blob).
        std::list<Blob *>::iterator it = m_joinedBlobs.begin(),en = m_joinedBlobs.end();
        m_externalPerimeter = (*it)->externalPerimeter(maskImage, xBorder, yBorder);
        it++;
        for (; it != en; ++it) {
            m_externalPerimeter += (*it)->externalPerimeter(maskImage, xBorder, yBorder);
        }
        return m_externalPerimeter;
    }


    //! get contour pixels
    externalContour = m_externalContour.getContourPoints();
    m_externalPerimeter = 0;

    //! there are contour pixels?
    if (externalContour.size() == 0) {
        return m_externalPerimeter;
    }

    PointList::iterator it = externalContour.begin(), en = externalContour.end();

    previousPoint.x = -1;

    //! which contour pixels touch border?
    for (; it != en; ++it) {
        actualPoint = *(it);
        find = false;

        //! pixel is touching border?
        if (xBorder & ((actualPoint.x == 0) || (actualPoint.x == m_originalImageSize.width - 1)) ||
                yBorder & ((actualPoint.y == 0) || (actualPoint.y == m_originalImageSize.height - 1))) {
            find = true;
        }
        else {
            if (maskImage != NULL) {
                //! verify if some of 8-connected neighbours is black in mask
                char *pMask;

                pMask = (maskImage->imageData + actualPoint.x - 1 + (actualPoint.y - 1) * maskImage->widthStep);

                for (int i = 0; i < 3; i++, pMask++) {
                    if (*pMask == 0 && !find) {
                        find = true;
                        break;
                    }
                }

                if (!find) {
                    pMask = (maskImage->imageData + actualPoint.x - 1 + (actualPoint.y) * maskImage->widthStep);

                    for (int i = 0; i < 3; i++, pMask++) {
                        if (*pMask == 0 && !find) {
                            find = true;
                            break;
                        }
                    }
                }

                if (!find) {
                    pMask = (maskImage->imageData + actualPoint.x - 1 + (actualPoint.y + 1) * maskImage->widthStep);

                    for (int i = 0; i < 3; i++, pMask++) {
                        if (*pMask == 0 && !find) {
                            find = true;
                            break;
                        }
                    }
                }
            }
        }

        if (find) {
            if (previousPoint.x > 0) {
                delta = abs(previousPoint.x - actualPoint.x) + abs(previousPoint.y - actualPoint.y);
            }

            //! calculate separately each external contour segment
            if (delta > 2) {
                m_externalPerimeter += arcLength(externalPoints, false);

                externalPoints.clear();
                delta = 0;
                previousPoint.x = -1;
            }
            externalPoints.push_back(actualPoint);
            previousPoint = actualPoint;
        }

    }

    if (externalPoints.size() != 0) {
        m_externalPerimeter += arcLength(externalPoints, false);
    }
    //! divide by two because external points have one side inside the blob and the other outside
    //! Perimeter of external points counts both sides, so it must be divided
    m_externalPerimeter /= 2.0;
    return m_externalPerimeter;
}

double Blob::externalPerimeter(const cv::Mat& maskImage, bool xBorder /* = true */, bool yBorder /* = true */) {
    if (!maskImage.data) {
        return externalPerimeter(NULL, xBorder /* = true */, yBorder /* = true */);
    }
    else {
        IplImage temp = cvIplImage(maskImage);
        return externalPerimeter(&temp, xBorder /* = true */, yBorder /* = true */);
    }
}

//! Compute blob's moment (p,q up to MAX_CALCULATED_MOMENTS)
double Blob::moment(int p, int q, bool intContours /*= true*/)
{
    double moment=0;
    BlobContourList::iterator itContour;

    if (m_isJoined) {
        std::list<Blob *>::iterator it,en = m_joinedBlobs.end();
        for (it = m_joinedBlobs.begin(); it != en; ++it) {
            moment += (*it)->moment(p, q);
        }
    }
    else {
        moment = m_externalContour.getMoment(p, q);
        if (intContours) {
            itContour = m_internalContours.begin();
            while (itContour != m_internalContours.end()) {
                moment -= (*itContour)->getMoment(p, q);
                itContour++;
            }
        }
    }
    return moment;
}

/**
- FUNCTIONS: mean
- FUNCIONALITAT: Get blob mean color in input image
- PARAMETERS:
    - image: image from gray color are extracted
- RESULTAT:
    -
- RESTRICCIONS:
    -
- AUTOR: rborras
- DATE OF CREATION: 2008/05/06
- MODIFICATION: Data. Autor. Description.
*/
double Blob::mean(IplImage *image)
{
    //! Create a mask with same size as blob bounding box
    IplImage *mask;
    CvScalar mean, std;
    cv::Point offset;

    getBoundingBox();

    if (m_boundingBox.height == 0 ||m_boundingBox.width == 0 || !CV_IS_IMAGE(image)) {
        m_meanGray = 0;
        return m_meanGray;
    }

    //! apply ROI and mask to input image to compute mean gray and standard deviation
    mask = cvCreateImage(cvSize(m_boundingBox.width, m_boundingBox.height), IPL_DEPTH_8U, 1);
    cvSetZero(mask);

    offset.x = -m_boundingBox.x;
    offset.y = -m_boundingBox.y;

    cv::Mat mask_mat = cv::cvarrToMat(mask);

    //! If joined
    if (m_isJoined) {
        std::list<Blob *>::iterator it,en = m_joinedBlobs.end();
        for (it = m_joinedBlobs.begin(); it != en; ++it) {
            //cvDrawContours(mask, (*it)->m_externalContour.getContourPoints(), CV_RGB(255, 255, 255), CV_RGB(255, 255, 255),0, CV_FILLED, 8, offset);
            //vector<PointList> conts;
            //conts.push_back((*it)->m_externalContour.getContourPoints());
            drawContours(mask_mat,(*it)->m_externalContour.getContours(), -1, CV_RGB(255, 255, 255), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
            BlobContourList::iterator itint = (*it)->m_internalContours.begin();
            while(itint != (*it)->m_internalContours.end()) {
                //cvDrawContours(mask, (*itint).getContourPoints(), CV_RGB(0, 0, 0), CV_RGB(0, 0, 0),0, CV_FILLED, 8, offset);
                drawContours(mask_mat, (*itint)->getContours(), -1, CV_RGB(0, 0, 0), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
                itint++;
            }
        }
    }
    else {
        //! draw contours on mask
        //cvDrawContours(mask, m_externalContour.getContourPoints(), CV_RGB(255, 255, 255), CV_RGB(255, 255, 255), 0, CV_FILLED, 8, offset);
        //vector<PointList> conts;
        //conts.push_back(m_externalContour.getContourPoints());
        drawContours(mask_mat, m_externalContour.getContours(), -1, CV_RGB(255, 255, 255), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
        //! draw internal contours
        BlobContourList::iterator it = m_internalContours.begin();
        while(it != m_internalContours.end())
        {
            //cvDrawContours(mask, (*it).getContourPoints(), CV_RGB(0, 0, 0), CV_RGB(0, 0, 0),0, CV_FILLED, 8, offset);
            drawContours(mask_mat, (*it)->getContours(), -1, CV_RGB(0, 0, 0), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
            it++;
        }
    }

    cvSetImageROI(image, cvRect(m_boundingBox));
    cvAvgSdv(image, &mean, &std, mask);

    m_meanGray = mean.val[0];
    m_stdDevGray = std.val[0];

    cvReleaseImage(&mask);
    cvResetImageROI(image);

    return m_meanGray;
}

double Blob::mean(const cv::Mat& image) {
    IplImage temp = cvIplImage(image);
    return mean(&temp);
}

double Blob::stdDev(IplImage *image) {
    //! call mean calculation (where also standard deviation is calculated)
    mean(image);

    return m_stdDevGray;
}

double Blob::stdDev(const cv::Mat& image) {
    IplImage temp = cvIplImage(image);
    return stdDev(&temp);
}

//void Blob::meanStdDev(Mat image, double *mean, double *stddev)
//{
//	IplImage temp = cvIplImage(image);
//	mean(&temp);
//	*mean = m_meanGray;
//	*stddev = m_stdDevGray;
//	return;
//}

void Blob::meanStdDev(const cv::Mat& image, cv::Scalar& mean, cv::Scalar& stddev)
{
    getBoundingBox();
    //! Create a mask with same size as blob bounding box and set it to 0
    cv::Mat mask_mat = cv::Mat(m_boundingBox.height, m_boundingBox.width, CV_8UC1, cv::Scalar(0));
    cv::Point offset(-m_boundingBox.x, -m_boundingBox.y);

    //! If joined
    if (m_isJoined) {
        std::list<Blob *>::iterator it,en = m_joinedBlobs.end();
        for (it = m_joinedBlobs.begin();it!=en;it++) {
            drawContours(mask_mat, (*it)->m_externalContour.getContours(), -1, CV_RGB(255, 255, 255), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
            BlobContourList::iterator itint = (*it)->m_internalContours.begin();
            while (itint != (*it)->m_internalContours.end()) {
                drawContours(mask_mat, (*itint)->getContours(), -1, CV_RGB(0, 0, 0), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
                itint++;
            }
        }
    }
    else {
        //! draw contours on mask
        drawContours(mask_mat, m_externalContour.getContours(), -1, CV_RGB(255, 255, 255), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
        //! draw internal contours
        BlobContourList::iterator it = m_internalContours.begin();
        while (it != m_internalContours.end()) {
            drawContours(mask_mat, (*it)->getContours(), -1, CV_RGB(0, 0, 0), CV_FILLED, 8, cv::noArray(), 2147483647, offset);
            it++;
        }
    }

    cv::meanStdDev(image(m_boundingBox), mean, stddev, mask_mat);
}

/**
- FUNCTIONS: getBoundingBox
- FUNCIONALITAT: Get bounding box (without rotation) of a blob
- PARAMETERS:
    -
- RESULTAT:
    -
- RESTRICCIONS:
    -
- AUTOR: rborras
- DATE OF CREATION: 2008/05/06
- MODIFICATION: Data. Autor. Description.
*/
cv::Rect Blob::getBoundingBox()
{
    //! it is calculated?
    if (m_boundingBox.width != -1) {
        return m_boundingBox;
    }

    if (m_isJoined) {
        cv::Rect bigRect;
        bigRect.x = 1000000;
        bigRect.y = 1000000;
        bigRect.height = 0;
        bigRect.width = 0;
        std::list<Blob *>::iterator it, en = m_joinedBlobs.end();
        int maxX = 0, maxY = 0;
        for (it = m_joinedBlobs.begin(); it != en; ++it) {
            cv::Rect temp = (*it)->getBoundingBox();
            if (bigRect.x > temp.x) {
                bigRect.x = temp.x;
            }
            if (bigRect.y > temp.y) {
                bigRect.y = temp.y;
            }
            //! -1 in order to obtain the correct measure (I'm looking for the points)
            if (maxX < temp.x + temp.width) {
                maxX = temp.x + temp.width - 1;
            }
            if (maxY < temp.y + temp.height) {
                maxY = temp.y + temp.height - 1;
            }
        }
        //! +1 in order to obtain the correct measure
        bigRect.width = maxX - bigRect.x + 1;
        bigRect.height = maxY - bigRect.y + 1;
        m_boundingBox = bigRect;
        return m_boundingBox;
    }

    PointList externalContour;

    //! get contour pixels
    externalContour = m_externalContour.getContourPoints();

    //! it is an empty blob?
    m_boundingBox.x = 1000000;
    m_boundingBox.y = 1000000;
    m_boundingBox.width = 0;
    m_boundingBox.height = 0;


    PointList::iterator it = externalContour.begin(), en = externalContour.end();
    for (; it != en; ++it) {
        cv::Point& actualPoint = *it;

        m_boundingBox.x = MIN(actualPoint.x, m_boundingBox.x);
        m_boundingBox.y = MIN(actualPoint.y, m_boundingBox.y);

        m_boundingBox.width = MAX(actualPoint.x, m_boundingBox.width);
        m_boundingBox.height = MAX(actualPoint.y, m_boundingBox.height);
    }

    //! +1 in order to take into account single pixel blobs.
    //! In this case, a single pixel has a bounding box like Rect(x,y,1,1), which is ok with opencv functions
    m_boundingBox.width = m_boundingBox.width - m_boundingBox.x + 1;
    m_boundingBox.height = m_boundingBox.height - m_boundingBox.y + 1;

    return m_boundingBox;
}

/**
- FUNCTIONS: getEllipse
- FUNCIONALITAT: Calculates bounding ellipse of external contour points
- PARAMETERS:
    -
- RESULTAT:
    -
- RESTRICCIONS:
    -
- AUTOR: rborras
- DATE OF CREATION: 2008/05/06
- MODIFICATION: Data. Autor. Description.
- NOTA: Calculation is made using second order moment aproximation
*/
cv::RotatedRect Blob::getEllipse()
{
    //! it is calculated?
    if (m_ellipse.size.width != -1) {
        return m_ellipse;
    }

    double u00, u11, u01, u10, u20, u02, delta, num, den, temp;

    //! central moments calculation
    u00 = moment(0,0);

    //! empty blob?
    if (u00 <= 0) {
        m_ellipse.size.width = 0;
        m_ellipse.size.height = 0;
        m_ellipse.center.x = 0;
        m_ellipse.center.y = 0;
        m_ellipse.angle = 0;
        return m_ellipse;
    }
    u10 = moment(1, 0) / u00;
    u01 = moment(0, 1) / u00;

    u11 = -(moment(1, 1) - moment(1, 0) * moment(0, 1) / u00) / u00;
    u20 = (moment(2, 0) - moment(1, 0) * moment(1, 0) / u00) / u00;
    u02 = (moment(0, 2) - moment(0, 1) * moment(0, 1) / u00) / u00;


    //! elipse calculation
    delta = sqrt(4*u11*u11 + (u20-u02)*(u20-u02));
    m_ellipse.center.x = (float)u10;
    m_ellipse.center.y = (float)u01;

    temp = u20 + u02 + delta;
    if (temp > 0) {
        m_ellipse.size.width = (float)sqrt(2 * (u20 + u02 + delta));
    }
    else {
        m_ellipse.size.width = 0;
        return m_ellipse;
    }

    temp = u20 + u02 - delta;
    if (temp > 0) {
        m_ellipse.size.height = (float)sqrt(2 * (u20 + u02 - delta));
    }
    else {
        m_ellipse.size.height = 0;
        return m_ellipse;
    }

    //! elipse orientation
    if (u20 > u02) {
        num = u02 - u20 + sqrt((u02 - u20) * (u02 - u20) + 4 * u11 * u11);
        den = 2*u11;
    }
    else {
        num = 2*u11;
        den = u20 - u02 + sqrt((u20 - u02) * (u20 - u02) + 4 * u11 * u11);
    }

    if (num != 0 && den != 00) {
        m_ellipse.angle = (float)(180.0 + (180.0 / CV_PI) * atan(num / den));
    }
    else {
        m_ellipse.angle = 0;
    }

    return m_ellipse;
}

/**
- FUNCTION: fillBlob
- FUNCTIONALITY: 
    - Fills the blob with a specified colour
- PARAMETERS:
    - image: where to paint
    - color: colour to paint the blob
    - offset: point offset for drawing
    - intContours: do not paint the internal holes (leave them transm_parent)
    - srcImage: image from where to copy the internal holes contents
- RESULT: - modified input image
- RESTRICTIONS:
- AUTHOR: Ricard Borràs
- CREATION DATE: 25-05-2005.
- MODIFICATION: - sep/2013. Luca Nardelli. Added functionality to consider internal contours when filling the blob.
*/
void Blob::fillBlob(IplImage* image, cv::Scalar color, int offsetX , int offsetY, bool intContours, const IplImage* srcImage)
{
    if (srcImage == NULL) {
        fillBlob(cv::cvarrToMat(image), color, offsetX, offsetY, intContours, cv::Mat());
    }
    else {
        fillBlob(cv::cvarrToMat(image), color, offsetX, offsetY, intContours, cv::cvarrToMat(srcImage));
    }
}

void Blob::fillBlob(const cv::Mat& image, cv::Scalar color, int offsetX, int offsetY, bool intContours, const cv::Mat& srcImage) {
    CV_FUNCNAME("Blob::fillBlob");
    __CV_BEGIN__;
    if (srcImage.data && intContours) CV_ASSERT(image.size() == srcImage.size() && image.type() == srcImage.type());
    {
        cv::Rect bbox = getBoundingBox();
        cv::Point drawOffset(offsetX, offsetY);
        cv::Size imSz = image.size();
        if (bbox.x + offsetX + bbox.width >= imSz.width) {
            bbox.width = imSz.width - bbox.x - offsetX;
        }
        else if (bbox.x + offsetX < 0) {
            bbox.x = -offsetX;
            bbox.width = bbox.width + offsetX;
        }
        if (bbox.y+offsetY+bbox.height >= imSz.height) {
            bbox.height = imSz.height - bbox.y - offsetY;
        }
        else if (bbox.y + offsetY < 0) {
            bbox.y = -offsetY;
            bbox.height = bbox.height + offsetY;
        }
        if (bbox.width <0 || bbox.height < 0) {
            return;
        }
        if (bbox.width == 0) {
            bbox.width++;
        }
        if (bbox.height == 0) {
            bbox.height++;
        }
        if (m_isJoined) {
            std::list<Blob *>::iterator itBlob = m_joinedBlobs.begin(), enBlob = m_joinedBlobs.end();
            for (itBlob = m_joinedBlobs.begin(); itBlob != enBlob; ++itBlob) {
                (*itBlob)->fillBlob(image, color, offsetX, offsetY, intContours, srcImage);
            }
            //if (intContours) {
            //	Point offset(-bbox.x, -bbox.y);
            //	Size sz(bbox.width, bbox.height);
            //	Mat temp(sz,image.type());
            //	Mat mask(sz, CV_8UC1);
            //	mask.setTo(0);
            //	for (itBlob = m_joinedBlobs.begin(); itBlob != enBlob; ++itBlob) {
            //		Blob *curBlob = *itBlob;
            //		BlobContourList::iterator it = curBlob->m_internalContours.begin(),en = curBlob->m_internalContours.end();
            //		for (it; it != en; ++it) {
            //			drawContours(mask,(*it)->getContours(),-1,255,CV_FILLED,8,noArray(),2147483647,offset);
            //		}
            //	}
            //	srcImage(bbox).copyTo(temp,mask);
            //	for (itBlob = m_joinedBlobs.begin(); itBlob != enBlob; ++itBlob) {
            //		drawContours(image,(*itBlob)->m_externalContour.getContours(),-1,color,CV_FILLED,8,noArray(),2147483647,drawOffset);
            //	}
            //	temp.copyTo(image(bbox+drawOffset),mask);
            //	for (itBlob = m_joinedBlobs.begin(); itBlob != enBlob; ++itBlob) {
            //		Blob *curBlob = *itBlob;
            //		BlobContourList::const_iterator it = curBlob->m_internalContours.begin(),en = curBlob->m_internalContours.end();
            //		for (it;it!=en;it++) {
            //			drawContours(image,(*it)->getContours(),-1,color,1,8,noArray(),2147483647,drawOffset);
            //		}
            //	}
            //}
            //else {
            //	for (itBlob = m_joinedBlobs.begin();itBlob!=enBlob;itBlob++) {
            //		drawContours(image,(*itBlob)->m_externalContour.getContours(),-1,color,CV_FILLED,8,noArray(),2147483647,drawOffset);
            //	}
            //}
        }
        else {
            if (intContours) {
                cv::Point offset(-bbox.x, -bbox.y);
                BlobContourList::iterator it = m_internalContours.begin(), en = m_internalContours.end();
                cv::Mat temp(bbox.height, bbox.width, image.type());
                drawContours(image, m_externalContour.getContours(), -1, color, CV_FILLED, 8, cv::noArray(), 2147483647, drawOffset);
                if (srcImage.data) {
                    cv::Mat mask(bbox.height, bbox.width, CV_8UC1);
                    mask.setTo(0);
                    for (; it != en; ++it) {
                        drawContours(mask, (*it)->getContours(), -1, 255, CV_FILLED, 8, cv::noArray(), 2147483647, offset);
                    }
                    srcImage(bbox).copyTo(temp,mask);
                    cv::Mat image_roi = image(bbox + drawOffset);
                    temp.copyTo(image_roi, mask);
                }
                else {
                    for (; it != en; ++it) {
                        drawContours(image, (*it)->getContours(), -1,CV_RGB(0, 0, 0), CV_FILLED, 8, cv::noArray(), 2147483647, drawOffset);
                    }
                }

                for (it = m_internalContours.begin(); it != en; ++it) {
                    drawContours(image, (*it)->getContours(), -1, color, 1, 8, cv::noArray(), 2147483647, drawOffset);
                }
            }
            else {
                drawContours(image,m_externalContour.getContours(), -1, color, CV_FILLED, 8, cv::noArray(), 2147483647, drawOffset);
            }
        }
    }
    __CV_END__;
}

/**
- FUNCTION: getConvexHull
- FUNCTIONALITY: Calculates the convex hull polygon of the blob
- PARAMETERS:
    - dst: where to store the result
- RESULT:
    - true if no error ocurred
- RESTRICTIONS:
- AUTHOR: Ricard Borràs
- CREATION DATE: 25-05-2005.
- MODIFICATION: Date. Author. Description.
*/
void Blob::getConvexHull(Contours& hull)
{
    hull.clear();
    PointList extCont;
    hull.push_back(PointList());
    if (m_isJoined) {
        int numPts = 0;
        BlobList::iterator it, en = m_joinedBlobs.end();
        for (it = m_joinedBlobs.begin(); it != en; ++it) {
            numPts += (*it)->getExternalContour()->getContourPoints().size();
        }
        hull[0].reserve(numPts);
        extCont.reserve(numPts);
        for (it = m_joinedBlobs.begin(); it != en; ++it) {
            const PointList& pts = (*it)->getExternalContour()->getContourPoints();
            extCont.insert(extCont.end(), pts.begin(), pts.end());
        }
    }
    else {
        extCont = m_externalContour.getContourPoints();
    }
    cv::convexHull(extCont, hull[0], true, true);
}

/**
- FUNCTION: joinBlob
- FUNCTIONALITY: Joins the 2 blobs, creating another blob which contains the 2 joined ones
- PARAMETERS:
    - blob: blob to join with the calling one
- RESULT:
    - Joined blob
- RESTRICTIONS: Only external contours are added
- AUTHOR: Ricard Borràs
- CREATION DATE: 25-05-2005.
- MODIFICATION: 08-2013, Luca Nardelli & Saverio Murgia, Created a working version of the join blob function
*/
void Blob::joinBlob(Blob* blob)
{
    //! Check on m_storage in order to not add empty blobs.
    if (!m_isJoined && !isEmpty()) {
        this->m_joinedBlobs.push_back(new Blob(this));
    }
    if (blob->m_isJoined) {
        std::list<Blob *>::iterator it, en = blob->m_joinedBlobs.end();
        for (it = blob->m_joinedBlobs.begin(); it != en; ++it) {
            this->m_joinedBlobs.push_back(new Blob(*it));
        }
    }
    else {
        this->m_joinedBlobs.push_back(new Blob(blob));
    }

    this->m_isJoined = true;
    this->m_boundingBox.width = -1;
    this->m_externalPerimeter = -1;
    this->m_meanGray = -1;
}

void Blob::requestDeletion(Blob* blob)
{
    //! If the blob has already been reported for deletion, then I also report the blobs that have requested it to be deleted
    while (blob->m_deleteRequestOwnerBlob != NULL) {
        Blob *temp = blob;
        blob->m_deleteRequestOwnerBlob->m_toBeDeleted = 1;
        blob = blob->m_deleteRequestOwnerBlob;
        temp->m_deleteRequestOwnerBlob = this;
    }
    blob->m_toBeDeleted = 1;
    blob->m_deleteRequestOwnerBlob = this;
}

int Blob::getNumJoinedBlobs()
{
    return m_joinedBlobs.size();
}

void Blob::shiftBlob(int x, int y)
{
    m_externalContour.shiftBlobContour(x,y);
    BlobContourList::iterator it, en = m_internalContours.end();
    for (it = m_internalContours.begin(); it != en; ++it) {
        (*it)->shiftBlobContour(x,y);
    }
    m_boundingBox.x += x;
    m_boundingBox.y += y;
}

cv::Point Blob::getCenter()
{
    return cv::Point((int)(getBoundingBox().x + getBoundingBox().width * 0.5), (int)(getBoundingBox().y + getBoundingBox().height * 0.5));
}

int Blob::overlappingPixels(Blob* blob)
{
    cv::Rect r1 = getBoundingBox(),r2 = blob->getBoundingBox();
    cv::Rect interRect = r1 & r2;
    if (interRect.width == 0 || interRect.height == 0) {
        return 0;
    }
    //! Minimum containing rectangle
    cv::Rect minContainingRect = r1 | r2;
    cv::Mat m1 = cv::Mat::zeros(minContainingRect.height, minContainingRect.width, CV_8UC1);
    cv::Mat m2 = cv::Mat::zeros(minContainingRect.height, minContainingRect.width, CV_8UC1);
    fillBlob(m1, cv::Scalar(255), -minContainingRect.x, -minContainingRect.y, true);
    blob->fillBlob(m2, cv::Scalar(255), -minContainingRect.x, -minContainingRect.y, true);
    return countNonZero(m1 & m2);
}

double Blob::density(AreaMode areaCalculationMode)
{
    double density = 0;
    BlobContourList::iterator itContour;
    switch (areaCalculationMode) {
    case GREEN:
    {
        double blobArea = area();
        Contours cHull, cHullPts;
        getConvexHull(cHull);
        //approxPolyDP(cHull, cHullPts, 0.001, true);
        double cHullArea = fabs(contourArea(cHull[0], false));
        density = blobArea/cHullArea;
        break;
    }
    case PIXELWISE:
    {
        Contours cHull;
        getConvexHull(cHull);
        cv::Rect bbox = getBoundingBox();
        cv::Mat blMat = cv::Mat::zeros(bbox.height, bbox.width, CV_8UC1);
        cv::Mat cHullMat = cv::Mat::zeros(bbox.height, bbox.width, CV_8UC1);
        fillBlob(blMat, cv::Scalar(255), -bbox.x, -bbox.y, true);
        drawContours(cHullMat, cHull, -1, cv::Scalar(255), -1, 8, cv::noArray(), 2147483647, cv::Point(-bbox.x, -bbox.y));
        int totArea = countNonZero(cHullMat);
        int actArea = countNonZero(blMat);
        density = (double)actArea/totArea;
        break;
    }
    }

    return density;
}

void Blob::getExtremes(cv::Point& xmax,cv::Point& xmin, cv::Point& ymax, cv::Point& ymin)
{
    const PointList& pts = m_externalContour.getContourPoints();
    /*Contours hull;
    getConvexHull(hull);
    const PointList &pts = hull[0];*/
    int nPts = pts.size();
    xmax = pts[0];
    xmin = pts[0];
    ymax = pts[0];
    ymin = pts[0];
    for (int i = 1; i < nPts; ++i) {
        if (xmax.x < pts[i].x) {
            xmax = pts[i];
        }
        else if (xmin.x > pts[i].x) {
            xmin = pts[i];
        }
        if (ymax.y < pts[i].y) {
            ymax = pts[i];
        }
        else if (ymin.y > pts[i].y) {
            ymin = pts[i];
        }
    }
}

cv::Point2f Blob::getCentroid(bool useInternalContours)
{
    float m00 = moment(0, 0, useInternalContours), m10 = moment(1, 0, useInternalContours);
    float m01 = moment(0, 1, useInternalContours);
    float x = m10 / m00, y = m01 / m00;
    return cv::Point2f(x, y);
}