structured_edge_detection.cpp

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                           License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's 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 Intel Corporation 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 Intel Corporation or 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.
//
//M*/

#include <vector>
#include <algorithm>
#include <iterator>
#include <iostream>
#include <cmath>

#include "precomp.hpp"

#include "advanced_types.hpp"

/********************* Helper functions *********************/

/*!
 * Lightweight wrapper over cv::resize
 *
 * \param src : source image to resize
 * \param dst : destination image size
 * \return resized image
 */
static cv::Mat imresize(const cv::Mat &src, const cv::Size &nSize)
{
    cv::Mat dst;
    if (nSize.width < src.size().width
    &&  nSize.height < src.size().height)
        cv::resize(src, dst, nSize, 0.0, 0.0, cv::INTER_AREA);
    else
        cv::resize(src, dst, nSize, 0.0, 0.0, cv::INTER_LINEAR);

    return dst;
}

/*!
 * The function filters src with triangle filter with radius equal rad
 *
 * \param src : source image to filter
 * \param rad : radius of filtering kernel
 * \return filtering result
 */
static cv::Mat imsmooth(const cv::Mat &src, const int rad)
{
    if (rad == 0)
        return src;
    else
    {
        const float p = 12.0f/rad/(rad + 2) - 2;
        cv::Mat dst;

        if (rad <= 1)
        {
            CV_INIT_VECTOR(kernelXY, float, {1/(p + 2), p/(p + 2), 1/(p + 2)});
            cv::sepFilter2D(src, dst, -1, kernelXY, kernelXY);
        }
        else
        {
            float nrml = CV_SQR(rad + 1.0f);

            std::vector <float> kernelXY(2*rad + 1);
            for (int i = 0; i <= rad; ++i)
            {
                kernelXY[2*rad - i] = (i + 1) / nrml;
                kernelXY[i] = (i + 1) / nrml;
            }
            sepFilter2D(src, dst, -1, kernelXY, kernelXY);
        }

        return dst;
    }
}

/*!
 *  The function implements rgb to luv conversion in a way similar
 *  to UCSD computer vision toolbox
 *
 * \param src : source image (RGB, float, in [0;1]) to convert
 * \return converted image in luv colorspace
 */
static cv::Mat rgb2luv(const cv::Mat &src)
{
    cv::Mat dst(src.size(), src.type());

    const float a  = CV_CUBE(29.0f)/27;
    const float y0 = 8.0f/a;

    const float mX[] = {0.430574f, 0.341550f, 0.178325f};
    const float mY[] = {0.222015f, 0.706655f, 0.071330f};
    const float mZ[] = {0.020183f, 0.129553f, 0.939180f};

    const float maxi= 1.0f/270;
    const float minu=  -88*maxi;
    const float minv= -134*maxi;

    const float un = 0.197833f;
    const float vn = 0.468331f;

    // build (padded) lookup table for y->l conversion assuming y in [0,1]
    std::vector <float> lTable(1024);
    for (int i = 0; i < 1024; ++i)
    {
        float y = i/1024.0f;
        float l = y > y0 ? 116*powf(y, 1.0f/3.0f) - 16 : y*a;

        lTable[i] = l*maxi;
    }
    for (int i = 0; i < 40; ++i)
        lTable.push_back(*--lTable.end());

    const int nchannels = 3;

    for (int i = 0; i < src.rows; ++i)
    {
        const float *pSrc = src.ptr<float>(i);
        float *pDst = dst.ptr<float>(i);

        for (int j = 0; j < src.cols*nchannels; j += nchannels)
        {
            const float rgb[] = {pSrc[j + 0], pSrc[j + 1], pSrc[j + 2]};

            const float xyz[] = {mX[0]*rgb[0] + mX[1]*rgb[1] + mX[2]*rgb[2],
                                 mY[0]*rgb[0] + mY[1]*rgb[1] + mY[2]*rgb[2],
                                 mZ[0]*rgb[0] + mZ[1]*rgb[1] + mZ[2]*rgb[2]};
            const float nz = 1.0f / float(xyz[0] + 15*xyz[1] + 3*xyz[2] + 1e-35);

            const float l = pDst[j] = lTable[cvFloor(1024*xyz[1])];

            pDst[j + 1] = l * (13*4*xyz[0]*nz - 13*un) - minu;;
            pDst[j + 2] = l * (13*9*xyz[1]*nz - 13*vn) - minv;
        }
    }

    return dst;
}

/*!
 * The function computes gradient magnitude and weighted (with magnitude)
 * orientation histogram. Magnitude is additionally normalized
 * by dividing on imsmooth(M, gnrmRad) + 0.01;
 *
 * \param src : source image
 * \param magnitude : gradient magnitude
 * \param histogram : gradient orientation nBins-channels histogram
 * \param nBins : number of gradient orientations
 * \param pSize : factor to downscale histogram
 * \param gnrmRad : radius for magnitude normalization
 */
static void gradientHist(const cv::Mat &src, cv::Mat &magnitude, cv::Mat &histogram,
                         const int nBins, const int pSize, const int gnrmRad)
{
    cv::Mat phase, Dx, Dy;

    magnitude.create( src.size(), cv::DataType<float>::type );
    phase.create( src.size(), cv::DataType<float>::type );
    histogram.create( cv::Size( cvRound(src.size().width/float(pSize)),
                                cvRound(src.size().height/float(pSize)) ),
        CV_MAKETYPE(cv::DataType<float>::type, nBins) );

    histogram.setTo(0);

    cv::Sobel( src, Dx, cv::DataType<float>::type,
        1, 0, 1, 1.0, 0.0, cv::BORDER_REFLECT );
    cv::Sobel( src, Dy, cv::DataType<float>::type,
        0, 1, 1, 1.0, 0.0, cv::BORDER_REFLECT );

    int nchannels = src.channels();

    for (int i = 0; i < src.rows; ++i)
    {
        const float *pDx = Dx.ptr<float>(i);
        const float *pDy = Dy.ptr<float>(i);

        float *pMagnitude = magnitude.ptr<float>(i);
        float *pPhase = phase.ptr<float>(i);

        for (int j = 0; j < src.cols*nchannels; j += nchannels)
        {
            float fMagn = float(-1e-5), fdx = 0, fdy = 0;
            for (int k = 0; k < nchannels; ++k)
            {
                float cMagn = CV_SQR( pDx[j + k] ) + CV_SQR( pDy[j + k] );
                if (cMagn > fMagn)
                {
                    fMagn = cMagn;
                    fdx = pDx[j + k];
                    fdy = pDy[j + k];
                }
            }

            pMagnitude[j/nchannels] = sqrtf(fMagn);

            float angle = cv::fastAtan2(fdy, fdx) / 180.0f - 1.0f * (fdy < 0);
            if (std::fabs(fdx) + std::fabs(fdy) < 1e-5)
                angle = 0.5f;
            pPhase[j/nchannels] = angle;
        }
    }

    magnitude /= imsmooth( magnitude, gnrmRad )
        + 0.01*cv::Mat::ones( magnitude.size(), magnitude.type() );

    for (int i = 0; i < phase.rows; ++i)
    {
        const float *pPhase = phase.ptr<float>(i);
        const float *pMagn  = magnitude.ptr<float>(i);

        float *pHist = histogram.ptr<float>(i/pSize);

        for (int j = 0; j < phase.cols; ++j)
            pHist[cvRound((j/pSize + pPhase[j])*nBins)] += pMagn[j] / CV_SQR(pSize);
    }
}

/********************* RFFeatureGetter class *********************/

namespace cv
{
namespace ximgproc
{

class RFFeatureGetterImpl : public RFFeatureGetter
{
public:
    /*!
     * Default constructor
     */
    RFFeatureGetterImpl() : name("RFFeatureGetter"){}

    /*!
     * The method extracts features from img and store them to features.
     * Extracted features are appropriate for StructuredEdgeDetection::predictEdges.
     *
     * \param src : source image (RGB, float, in [0;1]) to extract features
     * \param features : destination feature image
     *
     * \param gnrmRad : __rf.options.gradientNormalizationRadius
     * \param gsmthRad : __rf.options.gradientSmoothingRadius
     * \param shrink : __rf.options.shrinkNumber
     * \param outNum : __rf.options.numberOfOutputChannels
     * \param gradNum : __rf.options.numberOfGradientOrientations
     */
    virtual void getFeatures(const Mat &src, Mat &features, const int gnrmRad, const int gsmthRad,
                             const int shrink, const int outNum, const int gradNum) const
    {
        cv::Mat luvImg = rgb2luv(src);

        std::vector <cv::Mat> featureArray;

        cv::Size nSize = src.size() / float(shrink);
        split( imresize(luvImg, nSize), featureArray );

        CV_INIT_VECTOR(scales, float, {1.0f, 0.5f});

        for (size_t i = 0; i < scales.size(); ++i)
        {
            int pSize = std::max( 1, int(shrink*scales[i]) );

            cv::Mat magnitude, histogram;
            gradientHist(/**/ imsmooth(imresize(luvImg, scales[i]*src.size()), gsmthRad),
                magnitude, histogram, gradNum, pSize, gnrmRad /**/);

            featureArray.push_back(/**/ imresize( magnitude, nSize ).clone() /**/);
            featureArray.push_back(/**/ imresize( histogram, nSize ).clone() /**/);
        }

        // Mixing
        int resType = CV_MAKETYPE(cv::DataType<float>::type, outNum);
        features.create(nSize, resType);

        std::vector <int> fromTo;
        for (int i = 0; i < 2*outNum; ++i)
            fromTo.push_back(i/2);

        mixChannels(featureArray, features, fromTo);
    }

protected:
    /*! algorithm name */
    String name;
};

Ptr<RFFeatureGetter> createRFFeatureGetter()
{
        return makePtr<RFFeatureGetterImpl>();
}

}
}

/********************* StructuredEdgeDetection class *********************/

namespace cv
{
namespace ximgproc
{

class StructuredEdgeDetectionImpl : public StructuredEdgeDetection
{
public:
    /*!
     * This constructor loads __rf model from filename
     *
     * \param filename : name of the file where the model is stored
     */
    StructuredEdgeDetectionImpl(const cv::String &filename,
        Ptr<const RFFeatureGetter> _howToGetFeatures)
        : name("StructuredEdgeDetection"),
          howToGetFeatures( (!_howToGetFeatures.empty())
                          ? _howToGetFeatures
                          : createRFFeatureGetter().staticCast<const RFFeatureGetter>() )
    {
        cv::FileStorage modelFile(filename, FileStorage::READ);
        CV_Assert( modelFile.isOpened() );

        __rf.options.stride
            = modelFile["options"]["stride"];
        __rf.options.shrinkNumber
            = modelFile["options"]["shrinkNumber"];
        __rf.options.patchSize
            = modelFile["options"]["patchSize"];
        __rf.options.patchInnerSize
            = modelFile["options"]["patchInnerSize"];

        __rf.options.numberOfGradientOrientations
            = modelFile["options"]["numberOfGradientOrientations"];
        __rf.options.gradientSmoothingRadius
            = modelFile["options"]["gradientSmoothingRadius"];
        __rf.options.regFeatureSmoothingRadius
            = modelFile["options"]["regFeatureSmoothingRadius"];
        __rf.options.ssFeatureSmoothingRadius
            = modelFile["options"]["ssFeatureSmoothingRadius"];
        __rf.options.gradientNormalizationRadius
            = modelFile["options"]["gradientNormalizationRadius"];

        __rf.options.selfsimilarityGridSize
            = modelFile["options"]["selfsimilarityGridSize"];

        __rf.options.numberOfTrees
            = modelFile["options"]["numberOfTrees"];
        __rf.options.numberOfTreesToEvaluate
            = modelFile["options"]["numberOfTreesToEvaluate"];

        __rf.options.numberOfOutputChannels =
            2*(__rf.options.numberOfGradientOrientations + 1) + 3;
        //--------------------------------------------

        cv::FileNode childs = modelFile["childs"];
        cv::FileNode featureIds = modelFile["featureIds"];

        std::vector <int> currentTree;

        for(cv::FileNodeIterator it = childs.begin();
            it != childs.end(); ++it)
        {
            (*it) >> currentTree;
            std::copy(currentTree.begin(), currentTree.end(),
                std::back_inserter(__rf.childs));
        }

        for(cv::FileNodeIterator it = featureIds.begin();
            it != featureIds.end(); ++it)
        {
            (*it) >> currentTree;
            std::copy(currentTree.begin(), currentTree.end(),
                std::back_inserter(__rf.featureIds));
        }

        cv::FileNode thresholds = modelFile["thresholds"];
        std::vector <float> fcurrentTree;

        for(cv::FileNodeIterator it = thresholds.begin();
            it != thresholds.end(); ++it)
        {
            (*it) >> fcurrentTree;
            std::copy(fcurrentTree.begin(), fcurrentTree.end(),
                std::back_inserter(__rf.thresholds));
        }

        cv::FileNode edgeBoundaries = modelFile["edgeBoundaries"];
        cv::FileNode edgeBins = modelFile["edgeBins"];

        for(cv::FileNodeIterator it = edgeBoundaries.begin();
            it != edgeBoundaries.end(); ++it)
        {
            (*it) >> currentTree;
            std::copy(currentTree.begin(), currentTree.end(),
                std::back_inserter(__rf.edgeBoundaries));
        }

        for(cv::FileNodeIterator it = edgeBins.begin();
            it != edgeBins.end(); ++it)
        {
            (*it) >> currentTree;
            std::copy(currentTree.begin(), currentTree.end(),
                std::back_inserter(__rf.edgeBins));
        }

        __rf.numberOfTreeNodes = int( __rf.childs.size() ) / __rf.options.numberOfTrees;
    }

    /*!
     * The function detects edges in src and draw them to dst
     *
     * \param src : source image (RGB, float, in [0;1]) to detect edges
     * \param dst : destination image (grayscale, float, in [0;1])
     *              where edges are drawn
     */
    void detectEdges(const cv::Mat &src, cv::Mat &dst) const
    {
        CV_Assert( src.type() == CV_32FC3 );

        dst.create( src.size(), cv::DataType<float>::type );

        int padding = ( __rf.options.patchSize
            - __rf.options.patchInnerSize )/2;

        cv::Mat nSrc;
        copyMakeBorder( src, nSrc, padding, padding,
            padding, padding, BORDER_REFLECT );

        NChannelsMat features;
        createRFFeatureGetter()->getFeatures( nSrc, features,
            __rf.options.gradientNormalizationRadius,
            __rf.options.gradientSmoothingRadius,
            __rf.options.shrinkNumber,
            __rf.options.numberOfOutputChannels,
            __rf.options.numberOfGradientOrientations );
        predictEdges( features, dst );
    }

protected:
    /*!
     * Private method used by process method. The function
     * predict edges in n-channel feature image and store them to dst.
     *
     * \param features : source image (n-channels, float) to detect edges
     * \param dst : destination image (grayscale, float, in [0;1]) where edges are drawn
     */
    void predictEdges(const NChannelsMat &features, cv::Mat &dst) const
    {
        int shrink = __rf.options.shrinkNumber;
        int rfs = __rf.options.regFeatureSmoothingRadius;
        int sfs = __rf.options.ssFeatureSmoothingRadius;

        int nTreesEval = __rf.options.numberOfTreesToEvaluate;
        int nTrees = __rf.options.numberOfTrees;
        int nTreesNodes = __rf.numberOfTreeNodes;

        const int nchannels = features.channels();
        int pSize  = __rf.options.patchSize;

        int nFeatures = CV_SQR(pSize/shrink)*nchannels;
        int outNum = __rf.options.numberOfOutputChannels;

        int stride = __rf.options.stride;
        int ipSize = __rf.options.patchInnerSize;
        int gridSize = __rf.options.selfsimilarityGridSize;

        const int height = cvCeil( double(features.rows*shrink - pSize) / stride );
        const int width  = cvCeil( double(features.cols*shrink - pSize) / stride );
        // image size in patches with overlapping

        //-------------------------------------------------------------------------

        NChannelsMat regFeatures = imsmooth(features, cvRound(rfs / float(shrink)));
        NChannelsMat  ssFeatures = imsmooth(features, cvRound(sfs / float(shrink)));

        NChannelsMat indexes(height, width, CV_MAKETYPE(DataType<int>::type, nTreesEval));

        std::vector <int> offsetI(/**/ CV_SQR(pSize/shrink)*nchannels, 0);
        for (int i = 0; i < CV_SQR(pSize/shrink)*nchannels; ++i)
        {
            int z = i / CV_SQR(pSize/shrink);
            int y = ( i % CV_SQR(pSize/shrink) )/(pSize/shrink);
            int x = ( i % CV_SQR(pSize/shrink) )%(pSize/shrink);

            offsetI[i] = x*features.cols*nchannels + y*nchannels + z;
        }
        // lookup table for mapping linear index to offsets

        std::vector <int> offsetE(/**/ CV_SQR(ipSize)*outNum, 0);
        for (int i = 0; i < CV_SQR(ipSize)*outNum; ++i)
        {
            int z = i / CV_SQR(ipSize);
            int y = ( i % CV_SQR(ipSize) )/ipSize;
            int x = ( i % CV_SQR(ipSize) )%ipSize;

            offsetE[i] = x*dst.cols*outNum + y*outNum + z;
        }
        // lookup table for mapping linear index to offsets

        std::vector <int> offsetX( CV_SQR(gridSize)*(CV_SQR(gridSize) - 1)/2 * nchannels, 0);
        std::vector <int> offsetY( CV_SQR(gridSize)*(CV_SQR(gridSize) - 1)/2 * nchannels, 0);

        int hc = cvRound( (pSize/shrink) / (2.0*gridSize) );
        // half of cell
        std::vector <int> gridPositions;
        for(int i = 0; i < gridSize; i++)
            gridPositions.push_back( int( (i+1)*(pSize/shrink + 2*hc - 1)/(gridSize + 1.0) - hc + 0.5f ) );

        for (int i = 0, n = 0; i < CV_SQR(gridSize)*nchannels; ++i)
            for (int j = (i%CV_SQR(gridSize)) + 1; j < CV_SQR(gridSize); ++j, ++n)
            {
                int z = i / CV_SQR(gridSize);

                int x1 = gridPositions[i%CV_SQR(gridSize)%gridSize];
                int y1 = gridPositions[i%CV_SQR(gridSize)/gridSize];

                int x2 = gridPositions[j%gridSize];
                int y2 = gridPositions[j/gridSize];

                offsetX[n] = x1*features.cols*nchannels + y1*nchannels + z;
                offsetY[n] = x2*features.cols*nchannels + y2*nchannels + z;
            }
            // lookup tables for mapping linear index to offset pairs

        for (int i = 0; i < height; ++i)
        {
            float *regFeaturesPtr = regFeatures.ptr<float>(i*stride/shrink);
            float  *ssFeaturesPtr = ssFeatures.ptr<float>(i*stride/shrink);

            int *indexPtr = indexes.ptr<int>(i);

            for (int j = 0, k = 0; j < width; ++k, j += !(k %= nTreesEval))
                // for j,k in [0;width)x[0;nTreesEval)
            {
                int baseNode = ( ((i + j)%(2*nTreesEval) + k)%nTrees )*nTreesNodes;
                int currentNode = baseNode;
                // select root node of the tree to evaluate

                int offset = (j*stride/shrink)*nchannels;
                while ( __rf.childs[currentNode] != 0 )
                {
                    int currentId = __rf.featureIds[currentNode];
                    float currentFeature;

                    if (currentId >= nFeatures)
                    {
                        int xIndex = offsetX[currentId - nFeatures];
                        float A = ssFeaturesPtr[offset + xIndex];

                        int yIndex = offsetY[currentId - nFeatures];
                        float B = ssFeaturesPtr[offset + yIndex];

                        currentFeature = A - B;
                    }
                    else
                        currentFeature = regFeaturesPtr[offset + offsetI[currentId]];

                    // compare feature to threshold and move left or right accordingly
                    if (currentFeature < __rf.thresholds[currentNode])
                        currentNode = baseNode + __rf.childs[currentNode] - 1;
                    else
                        currentNode = baseNode + __rf.childs[currentNode];
                }

                indexPtr[j*nTreesEval + k] = currentNode;
            }
        }

        NChannelsMat dstM(dst.size(),
            CV_MAKETYPE(DataType<float>::type, outNum));
        dstM.setTo(0);

        float step = 2.0f * CV_SQR(stride) / CV_SQR(ipSize) / nTreesEval;
        for (int i = 0; i < height; ++i)
        {
            int *pIndex = indexes.ptr<int>(i);
            float *pDst = dstM.ptr<float>(i*stride);

                for (int j = 0, k = 0; j < width; ++k, j += !(k %= nTreesEval))
                {// for j,k in [0;width)x[0;nTreesEval)

                    int currentNode = pIndex[j*nTreesEval + k];

                    int start  = __rf.edgeBoundaries[currentNode];
                    int finish = __rf.edgeBoundaries[currentNode + 1];

                    if (start == finish)
                        continue;

                    int offset = j*stride*outNum;
                    for (int p = start; p < finish; ++p)
                        pDst[offset + offsetE[__rf.edgeBins[p]]] += step;
                }
        }

        cv::reduce( dstM.reshape(1, int( dstM.total() ) ), dstM, 2, CV_REDUCE_SUM);
        imsmooth( dstM.reshape(1, dst.rows), 1 ).copyTo(dst);
    }

/********************* Members *********************/
protected:
    /*! algorithm name */
    String name;

    /*! optional feature getter (getFeatures method) */
    Ptr<const RFFeatureGetter> howToGetFeatures;

    /*! random forest used to detect edges */
    struct RandomForest
    {
        /*! random forest options, e.g. number of trees */
        struct RandomForestOptions
        {
            // model params

            int numberOfOutputChannels; /*!< number of edge orientation bins for output */

            int patchSize;              /*!< width of image patches */
            int patchInnerSize;         /*!< width of predicted part inside patch*/

            // feature params

            int regFeatureSmoothingRadius;    /*!< radius for smoothing of regular features
                                               *   (using convolution with triangle filter) */

            int ssFeatureSmoothingRadius;     /*!< radius for smoothing of additional features
                                               *   (using convolution with triangle filter) */

            int shrinkNumber;                 /*!< amount to shrink channels */

            int numberOfGradientOrientations; /*!< number of orientations per gradient scale */

            int gradientSmoothingRadius;      /*!< radius for smoothing of gradients
                                               *   (using convolution with triangle filter) */

            int gradientNormalizationRadius;  /*!< gradient normalization radius */
            int selfsimilarityGridSize;       /*!< number of self similarity cells */

            // detection params
            int numberOfTrees;            /*!< number of trees in forest to train */
            int numberOfTreesToEvaluate;  /*!< number of trees to evaluate per location */

            int stride;                   /*!< stride at which to compute edges */

        } options;

        int numberOfTreeNodes;

        std::vector <int> featureIds;     /*!< feature coordinate thresholded at k-th node */
        std::vector <float> thresholds;   /*!< threshold applied to featureIds[k] at k-th node */
        std::vector <int> childs;         /*!< k --> child[k] - 1, child[k] */

        std::vector <int> edgeBoundaries; /*!< ... */
        std::vector <int> edgeBins;       /*!< ... */
    } __rf;
};

Ptr<StructuredEdgeDetection> createStructuredEdgeDetection(const String &model,
    Ptr<const RFFeatureGetter> howToGetFeatures)
{
        return makePtr<StructuredEdgeDetectionImpl>(model, howToGetFeatures);
}

}
}