added Generalized Hough implementation

modules/gpu/include/opencv2/gpu/gpu.hpp

@@ -770,11 +770,11 @@ CV_EXPORTS void blendLinear(const GpuMat& img1, const GpuMat& img2, const GpuMat
                             GpuMat& result, Stream& stream = Stream::Null());
 
 //! Performa bilateral filtering of passsed image
-CV_EXPORTS void bilateralFilter(const GpuMat& src, GpuMat& dst, int kernel_size, float sigma_color, float sigma_spatial, 
+CV_EXPORTS void bilateralFilter(const GpuMat& src, GpuMat& dst, int kernel_size, float sigma_color, float sigma_spatial,
                                 int borderMode = BORDER_DEFAULT, Stream& stream = Stream::Null());
 
 //! Brute force non-local means algorith (slow but universal)
-CV_EXPORTS void nonLocalMeans(const GpuMat& src, GpuMat& dst, float h, 
+CV_EXPORTS void nonLocalMeans(const GpuMat& src, GpuMat& dst, float h,
                               int search_widow_size = 11, int block_size = 7, int borderMode = BORDER_DEFAULT, Stream& s = Stream::Null());
 
 
@@ -854,6 +854,38 @@ CV_EXPORTS void HoughCircles(const GpuMat& src, GpuMat& circles, int method, flo
 CV_EXPORTS void HoughCircles(const GpuMat& src, GpuMat& circles, HoughCirclesBuf& buf, int method, float dp, float minDist, int cannyThreshold, int votesThreshold, int minRadius, int maxRadius, int maxCircles = 4096);
 CV_EXPORTS void HoughCirclesDownload(const GpuMat& d_circles, OutputArray h_circles);
 
+//! finds arbitrary template in the grayscale image using Generalized Hough Transform
+//! Ballard, D.H. (1981). Generalizing the Hough transform to detect arbitrary shapes. Pattern Recognition 13 (2): 111-122.
+//! Guil, N., González-Linares, J.M. and Zapata, E.L. (1999). Bidimensional shape detection using an invariant approach. Pattern Recognition 32 (6): 1025-1038.
+class CV_EXPORTS GeneralizedHough_GPU : public Algorithm
+{
+public:
+    static Ptr<GeneralizedHough_GPU> create(int method);
+
+    virtual ~GeneralizedHough_GPU();
+
+    //! set template to search
+    void setTemplate(const GpuMat& templ, int cannyThreshold = 100, Point templCenter = Point(-1, -1));
+    void setTemplate(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, Point templCenter = Point(-1, -1));
+
+    //! find template on image
+    void detect(const GpuMat& image, GpuMat& positions, int cannyThreshold = 100);
+    void detect(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, GpuMat& positions);
+
+    void download(const GpuMat& d_positions, OutputArray h_positions, OutputArray h_votes = noArray());
+
+    void release();
+
+protected:
+    virtual void setTemplateImpl(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, Point templCenter) = 0;
+    virtual void detectImpl(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, GpuMat& positions) = 0;
+    virtual void releaseImpl() = 0;
+
+private:
+    GpuMat edges_;
+    CannyBuf cannyBuf_;
+};
+
 ////////////////////////////// Matrix reductions //////////////////////////////
 
 //! computes mean value and standard deviation of all or selected array elements

modules/gpu/perf/perf_imgproc.cpp

@@ -1713,4 +1713,98 @@ PERF_TEST_P(Sz_Dp_MinDist, ImgProc_HoughCircles, Combine(GPU_TYPICAL_MAT_SIZES,
     }
 }
 
+//////////////////////////////////////////////////////////////////////
+// GeneralizedHough
+
+CV_FLAGS(GHMethod, cv::GHT_POSITION, cv::GHT_SCALE, cv::GHT_ROTATION);
+
+DEF_PARAM_TEST(Method_Sz, GHMethod, cv::Size);
+
+PERF_TEST_P(Method_Sz, GeneralizedHough, Combine(
+            Values(GHMethod(cv::GHT_POSITION), GHMethod(cv::GHT_POSITION | cv::GHT_SCALE), GHMethod(cv::GHT_POSITION | cv::GHT_ROTATION), GHMethod(cv::GHT_POSITION | cv::GHT_SCALE | cv::GHT_ROTATION)),
+            GPU_TYPICAL_MAT_SIZES))
+{
+    declare.time(10);
+
+    const int method = GET_PARAM(0);
+    const cv::Size imageSize = GET_PARAM(1);
+
+    const cv::Mat templ = readImage("cv/shared/templ.png", cv::IMREAD_GRAYSCALE);
+    ASSERT_FALSE(templ.empty());
+
+    cv::Mat image(imageSize, CV_8UC1, cv::Scalar::all(0));
+
+    cv::RNG rng(123456789);
+    const int objCount = rng.uniform(5, 15);
+    for (int i = 0; i < objCount; ++i)
+    {
+        double scale = rng.uniform(0.7, 1.3);
+        bool rotate = rng.uniform(0, 2);
+
+        cv::Mat obj;
+        cv::resize(templ, obj, cv::Size(), scale, scale);
+        if (rotate)
+            obj = obj.t();
+
+        cv::Point pos;
+
+        pos.x = rng.uniform(0, image.cols - obj.cols);
+        pos.y = rng.uniform(0, image.rows - obj.rows);
+
+        cv::Mat roi = image(cv::Rect(pos, obj.size()));
+        cv::add(roi, obj, roi);
+    }
+
+    cv::Mat edges;
+    cv::Canny(image, edges, 50, 100);
+
+    cv::Mat dx, dy;
+    cv::Sobel(image, dx, CV_32F, 1, 0);
+    cv::Sobel(image, dy, CV_32F, 0, 1);
+
+    if (runOnGpu)
+    {
+        cv::gpu::GpuMat d_edges(edges);
+        cv::gpu::GpuMat d_dx(dx);
+        cv::gpu::GpuMat d_dy(dy);
+        cv::gpu::GpuMat d_position;
+
+        cv::Ptr<cv::gpu::GeneralizedHough_GPU> d_hough = cv::gpu::GeneralizedHough_GPU::create(method);
+        if (method & cv::GHT_ROTATION)
+        {
+            d_hough->set("maxAngle", 90.0);
+            d_hough->set("angleStep", 2.0);
+        }
+
+        d_hough->setTemplate(cv::gpu::GpuMat(templ));
+
+        d_hough->detect(d_edges, d_dx, d_dy, d_position);
+
+        TEST_CYCLE()
+        {
+            d_hough->detect(d_edges, d_dx, d_dy, d_position);
+        }
+    }
+    else
+    {
+        cv::Mat positions;
+
+        cv::Ptr<cv::GeneralizedHough> hough = cv::GeneralizedHough::create(method);
+        if (method & cv::GHT_ROTATION)
+        {
+            hough->set("maxAngle", 90.0);
+            hough->set("angleStep", 2.0);
+        }
+
+        hough->setTemplate(templ);
+
+        hough->detect(edges, dx, dy, positions);
+
+        TEST_CYCLE()
+        {
+            hough->detect(edges, dx, dy, positions);
+        }
+    }
+}
+
 } // namespace

modules/gpu/src/cuda/hough.cu

@@ -43,6 +43,9 @@
 #include <thrust/sort.h>
 #include "opencv2/gpu/device/common.hpp"
 #include "opencv2/gpu/device/emulation.hpp"
+#include "opencv2/gpu/device/vec_math.hpp"
+#include "opencv2/gpu/device/limits.hpp"
+#include "opencv2/gpu/device/dynamic_smem.hpp"
 
 namespace cv { namespace gpu { namespace device
 {
@@ -53,8 +56,7 @@ namespace cv { namespace gpu { namespace device
         ////////////////////////////////////////////////////////////////////////
         // buildPointList
 
-        const int PIXELS_PER_THREAD = 16;
-
+        template <int PIXELS_PER_THREAD>
         __global__ void buildPointList(const PtrStepSzb src, unsigned int* list)
         {
             __shared__ unsigned int s_queues[4][32 * PIXELS_PER_THREAD];
@@ -113,6 +115,8 @@ namespace cv { namespace gpu { namespace device
 
         int buildPointList_gpu(PtrStepSzb src, unsigned int* list)
         {
+            const int PIXELS_PER_THREAD = 16;
+
             void* counterPtr;
             cudaSafeCall( cudaGetSymbolAddress(&counterPtr, g_counter) );
 
@@ -121,9 +125,9 @@ namespace cv { namespace gpu { namespace device
             const dim3 block(32, 4);
             const dim3 grid(divUp(src.cols, block.x * PIXELS_PER_THREAD), divUp(src.rows, block.y));
 
-            cudaSafeCall( cudaFuncSetCacheConfig(buildPointList, cudaFuncCachePreferShared) );
+            cudaSafeCall( cudaFuncSetCacheConfig(buildPointList<PIXELS_PER_THREAD>, cudaFuncCachePreferShared) );
 
-            buildPointList<<<grid, block>>>(src, list);
+            buildPointList<PIXELS_PER_THREAD><<<grid, block>>>(src, list);
             cudaSafeCall( cudaGetLastError() );
 
             cudaSafeCall( cudaDeviceSynchronize() );
@@ -167,7 +171,7 @@ namespace cv { namespace gpu { namespace device
 
         __global__ void linesAccumShared(const unsigned int* list, const int count, PtrStepi accum, const float irho, const float theta, const int numrho)
         {
-            extern __shared__ int smem[];
+            int* smem = DynamicSharedMem<int>();
 
             for (int i = threadIdx.x; i < numrho + 1; i += blockDim.x)
                 smem[i] = 0;
@@ -410,7 +414,7 @@ namespace cv { namespace gpu { namespace device
                                            float3* circles, const int maxCircles, const float dp,
                                            const int minRadius, const int maxRadius, const int histSize, const int threshold)
         {
-            extern __shared__ int smem[];
+            int* smem = DynamicSharedMem<int>();
 
             for (int i = threadIdx.x; i < histSize + 2; i += blockDim.x)
                 smem[i] = 0;
@@ -481,5 +485,1023 @@ namespace cv { namespace gpu { namespace device
 
             return totalCount;
         }
+
+        ////////////////////////////////////////////////////////////////////////
+        // Generalized Hough
+
+        template <typename T, int PIXELS_PER_THREAD>
+        __global__ void buildEdgePointList(const PtrStepSzb edges, const PtrStep<T> dx, const PtrStep<T> dy, unsigned int* coordList, float* thetaList)
+        {
+            __shared__ unsigned int s_coordLists[4][32 * PIXELS_PER_THREAD];
+            __shared__ float s_thetaLists[4][32 * PIXELS_PER_THREAD];
+            __shared__ int s_sizes[4];
+            __shared__ int s_globStart[4];
+
+            const int x = blockIdx.x * blockDim.x * PIXELS_PER_THREAD + threadIdx.x;
+            const int y = blockIdx.y * blockDim.y + threadIdx.y;
+
+            if (threadIdx.x == 0)
+                s_sizes[threadIdx.y] = 0;
+            __syncthreads();
+
+            if (y < edges.rows)
+            {
+                // fill the queue
+                const uchar* edgesRow = edges.ptr(y);
+                const T* dxRow = dx.ptr(y);
+                const T* dyRow = dy.ptr(y);
+
+                for (int i = 0, xx = x; i < PIXELS_PER_THREAD && xx < edges.cols; ++i, xx += blockDim.x)
+                {
+                    const T dxVal = dxRow[xx];
+                    const T dyVal = dyRow[xx];
+
+                    if (edgesRow[xx] && (dxVal != 0 || dyVal != 0))
+                    {
+                        const unsigned int coord = (y << 16) | xx;
+
+                        float theta = ::atan2f(dyVal, dxVal);
+                        if (theta < 0)
+                            theta += 2.0f * CV_PI_F;
+
+                        const int qidx = Emulation::smem::atomicAdd(&s_sizes[threadIdx.y], 1);
+
+                        s_coordLists[threadIdx.y][qidx] = coord;
+                        s_thetaLists[threadIdx.y][qidx] = theta;
+                    }
+                }
+            }
+
+            __syncthreads();
+
+            // let one thread reserve the space required in the global list
+            if (threadIdx.x == 0 && threadIdx.y == 0)
+            {
+                // find how many items are stored in each list
+                int totalSize = 0;
+                for (int i = 0; i < blockDim.y; ++i)
+                {
+                    s_globStart[i] = totalSize;
+                    totalSize += s_sizes[i];
+                }
+
+                // calculate the offset in the global list
+                const int globalOffset = atomicAdd(&g_counter, totalSize);
+                for (int i = 0; i < blockDim.y; ++i)
+                    s_globStart[i] += globalOffset;
+            }
+
+            __syncthreads();
+
+            // copy local queues to global queue
+            const int qsize = s_sizes[threadIdx.y];
+            int gidx = s_globStart[threadIdx.y] + threadIdx.x;
+            for(int i = threadIdx.x; i < qsize; i += blockDim.x, gidx += blockDim.x)
+            {
+                coordList[gidx] = s_coordLists[threadIdx.y][i];
+                thetaList[gidx] = s_thetaLists[threadIdx.y][i];
+            }
+        }
+
+        template <typename T>
+        int buildEdgePointList_gpu(PtrStepSzb edges, PtrStepSzb dx, PtrStepSzb dy, unsigned int* coordList, float* thetaList)
+        {
+            const int PIXELS_PER_THREAD = 8;
+
+            void* counterPtr;
+            cudaSafeCall( cudaGetSymbolAddress(&counterPtr, g_counter) );
+
+            cudaSafeCall( cudaMemset(counterPtr, 0, sizeof(int)) );
+
+            const dim3 block(32, 4);
+            const dim3 grid(divUp(edges.cols, block.x * PIXELS_PER_THREAD), divUp(edges.rows, block.y));
+
+            cudaSafeCall( cudaFuncSetCacheConfig(buildEdgePointList<T, PIXELS_PER_THREAD>, cudaFuncCachePreferShared) );
+
+            buildEdgePointList<T, PIXELS_PER_THREAD><<<grid, block>>>(edges, (PtrStepSz<T>) dx, (PtrStepSz<T>) dy, coordList, thetaList);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+
+            int totalCount;
+            cudaSafeCall( cudaMemcpy(&totalCount, counterPtr, sizeof(int), cudaMemcpyDeviceToHost) );
+
+            return totalCount;
+        }
+
+        template int buildEdgePointList_gpu<short>(PtrStepSzb edges, PtrStepSzb dx, PtrStepSzb dy, unsigned int* coordList, float* thetaList);
+        template int buildEdgePointList_gpu<int>(PtrStepSzb edges, PtrStepSzb dx, PtrStepSzb dy, unsigned int* coordList, float* thetaList);
+        template int buildEdgePointList_gpu<float>(PtrStepSzb edges, PtrStepSzb dx, PtrStepSzb dy, unsigned int* coordList, float* thetaList);
+
+        __global__ void buildRTable(const unsigned int* coordList, const float* thetaList, const int pointsCount,
+                                    PtrStep<short2> r_table, int* r_sizes, int maxSize,
+                                    const short2 templCenter, const float thetaScale)
+        {
+            const int tid = blockIdx.x * blockDim.x + threadIdx.x;
+
+            if (tid >= pointsCount)
+                return;
+
+            const unsigned int coord = coordList[tid];
+            short2 p;
+            p.x = (coord & 0xFFFF);
+            p.y = (coord >> 16) & 0xFFFF;
+
+            const float theta = thetaList[tid];
+            const int n = __float2int_rn(theta * thetaScale);
+
+            const int ind = ::atomicAdd(r_sizes + n, 1);
+            if (ind < maxSize)
+                r_table(n, ind) = p - templCenter;
+        }
+
+        void buildRTable_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                             PtrStepSz<short2> r_table, int* r_sizes,
+                             short2 templCenter, int levels)
+        {
+            const dim3 block(256);
+            const dim3 grid(divUp(pointsCount, block.x));
+
+            const float thetaScale = levels / (2.0f * CV_PI_F);
+
+            buildRTable<<<grid, block>>>(coordList, thetaList, pointsCount, r_table, r_sizes, r_table.cols, templCenter, thetaScale);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        ////////////////////////////////////////////////////////////////////////
+        // GHT_Ballard_Pos
+
+        __global__ void GHT_Ballard_Pos_calcHist(const unsigned int* coordList, const float* thetaList, const int pointsCount,
+                                                 const PtrStep<short2> r_table, const int* r_sizes,
+                                                 PtrStepSzi hist,
+                                                 const float idp, const float thetaScale)
+        {
+            const int tid = blockIdx.x * blockDim.x + threadIdx.x;
+
+            if (tid >= pointsCount)
+                return;
+
+            const unsigned int coord = coordList[tid];
+            short2 p;
+            p.x = (coord & 0xFFFF);
+            p.y = (coord >> 16) & 0xFFFF;
+
+            const float theta = thetaList[tid];
+            const int n = __float2int_rn(theta * thetaScale);
+
+            const short2* r_row = r_table.ptr(n);
+            const int r_row_size = r_sizes[n];
+
+            for (int j = 0; j < r_row_size; ++j)
+            {
+                short2 c = p - r_row[j];
+
+                c.x = __float2int_rn(c.x * idp);
+                c.y = __float2int_rn(c.y * idp);
+
+                if (c.x >= 0 && c.x < hist.cols - 2 && c.y >= 0 && c.y < hist.rows - 2)
+                    ::atomicAdd(hist.ptr(c.y + 1) + c.x + 1, 1);
+            }
+        }
+
+        void GHT_Ballard_Pos_calcHist_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                          PtrStepSz<short2> r_table, const int* r_sizes,
+                                          PtrStepSzi hist,
+                                          float dp, int levels)
+        {
+            const dim3 block(256);
+            const dim3 grid(divUp(pointsCount, block.x));
+
+            const float idp = 1.0f / dp;
+            const float thetaScale = levels / (2.0f * CV_PI_F);
+
+            GHT_Ballard_Pos_calcHist<<<grid, block>>>(coordList, thetaList, pointsCount, r_table, r_sizes, hist, idp, thetaScale);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        __global__ void GHT_Ballard_Pos_findPosInHist(const PtrStepSzi hist, float4* out, int3* votes, const int maxSize, const float dp, const int threshold)
+        {
+            const int x = blockIdx.x * blockDim.x + threadIdx.x;
+            const int y = blockIdx.y * blockDim.y + threadIdx.y;
+
+            if (x >= hist.cols - 2 || y >= hist.rows - 2)
+                return;
+
+            const int curVotes = hist(y + 1, x + 1);
+
+            if (curVotes > threshold &&
+                curVotes >  hist(y + 1, x) &&
+                curVotes >= hist(y + 1, x + 2) &&
+                curVotes >  hist(y, x + 1) &&
+                curVotes >= hist(y + 2, x + 1))
+            {
+                const int ind = ::atomicAdd(&g_counter, 1);
+
+                if (ind < maxSize)
+                {
+                    out[ind] = make_float4(x * dp, y * dp, 1.0f, 0.0f);
+                    votes[ind] = make_int3(curVotes, 0, 0);
+                }
+            }
+        }
+
+        int GHT_Ballard_Pos_findPosInHist_gpu(PtrStepSzi hist, float4* out, int3* votes, int maxSize, float dp, int threshold)
+        {
+            void* counterPtr;
+            cudaSafeCall( cudaGetSymbolAddress(&counterPtr, g_counter) );
+
+            cudaSafeCall( cudaMemset(counterPtr, 0, sizeof(int)) );
+
+            const dim3 block(32, 8);
+            const dim3 grid(divUp(hist.cols - 2, block.x), divUp(hist.rows - 2, block.y));
+
+            cudaSafeCall( cudaFuncSetCacheConfig(GHT_Ballard_Pos_findPosInHist, cudaFuncCachePreferL1) );
+
+            GHT_Ballard_Pos_findPosInHist<<<grid, block>>>(hist, out, votes, maxSize, dp, threshold);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+
+            int totalCount;
+            cudaSafeCall( cudaMemcpy(&totalCount, counterPtr, sizeof(int), cudaMemcpyDeviceToHost) );
+
+            totalCount = ::min(totalCount, maxSize);
+
+            return totalCount;
+        }
+
+        ////////////////////////////////////////////////////////////////////////
+        // GHT_Ballard_PosScale
+
+        __global__ void GHT_Ballard_PosScale_calcHist(const unsigned int* coordList, const float* thetaList,
+                                                      PtrStep<short2> r_table, const int* r_sizes,
+                                                      PtrStepi hist, const int rows, const int cols,
+                                                      const float minScale, const float scaleStep, const int scaleRange,
+                                                      const float idp, const float thetaScale)
+        {
+            const unsigned int coord = coordList[blockIdx.x];
+            float2 p;
+            p.x = (coord & 0xFFFF);
+            p.y = (coord >> 16) & 0xFFFF;
+
+            const float theta = thetaList[blockIdx.x];
+            const int n = __float2int_rn(theta * thetaScale);
+
+            const short2* r_row = r_table.ptr(n);
+            const int r_row_size = r_sizes[n];
+
+            for (int j = 0; j < r_row_size; ++j)
+            {
+                const float2 d = saturate_cast<float2>(r_row[j]);
+
+                for (int s = threadIdx.x; s < scaleRange; s += blockDim.x)
+                {
+                    const float scale = minScale + s * scaleStep;
+
+                    float2 c = p - scale * d;
+
+                    c.x *= idp;
+                    c.y *= idp;
+
+                    if (c.x >= 0 && c.x < cols && c.y >= 0 && c.y < rows)
+                        ::atomicAdd(hist.ptr((s + 1) * (rows + 2) + __float2int_rn(c.y + 1)) + __float2int_rn(c.x + 1), 1);
+                }
+            }
+        }
+
+        void GHT_Ballard_PosScale_calcHist_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                               PtrStepSz<short2> r_table, const int* r_sizes,
+                                               PtrStepi hist, int rows, int cols,
+                                               float minScale, float scaleStep, int scaleRange,
+                                               float dp, int levels)
+        {
+            const dim3 block(256);
+            const dim3 grid(pointsCount);
+
+            const float idp = 1.0f / dp;
+            const float thetaScale = levels / (2.0f * CV_PI_F);
+
+            GHT_Ballard_PosScale_calcHist<<<grid, block>>>(coordList, thetaList,
+                                                           r_table, r_sizes,
+                                                           hist, rows, cols,
+                                                           minScale, scaleStep, scaleRange,
+                                                           idp, thetaScale);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        __global__ void GHT_Ballard_PosScale_findPosInHist(const PtrStepi hist, const int rows, const int cols, const int scaleRange,
+                                                           float4* out, int3* votes, const int maxSize,
+                                                           const float minScale, const float scaleStep, const float dp, const int threshold)
+        {
+            const int x = blockIdx.x * blockDim.x + threadIdx.x;
+            const int y = blockIdx.y * blockDim.y + threadIdx.y;
+
+            if (x >= cols || y >= rows)
+                return;
+
+            for (int s = 0; s < scaleRange; ++s)
+            {
+                const float scale = minScale + s * scaleStep;
+
+                const int prevScaleIdx = (s) * (rows + 2);
+                const int curScaleIdx = (s + 1) * (rows + 2);
+                const int nextScaleIdx = (s + 2) * (rows + 2);
+
+                const int curVotes = hist(curScaleIdx + y + 1, x + 1);
+
+                if (curVotes > threshold &&
+                    curVotes >  hist(curScaleIdx + y + 1, x) &&
+                    curVotes >= hist(curScaleIdx + y + 1, x + 2) &&
+                    curVotes >  hist(curScaleIdx + y, x + 1) &&
+                    curVotes >= hist(curScaleIdx + y + 2, x + 1) &&
+                    curVotes >  hist(prevScaleIdx + y + 1, x + 1) &&
+                    curVotes >= hist(nextScaleIdx + y + 1, x + 1))
+                {
+                    const int ind = ::atomicAdd(&g_counter, 1);
+
+                    if (ind < maxSize)
+                    {
+                        out[ind] = make_float4(x * dp, y * dp, scale, 0.0f);
+                        votes[ind] = make_int3(curVotes, curVotes, 0);
+                    }
+                }
+            }
+        }
+
+        int GHT_Ballard_PosScale_findPosInHist_gpu(PtrStepi hist, int rows, int cols, int scaleRange, float4* out, int3* votes, int maxSize,
+                                                   float minScale, float scaleStep, float dp, int threshold)
+        {
+            void* counterPtr;
+            cudaSafeCall( cudaGetSymbolAddress(&counterPtr, g_counter) );
+
+            cudaSafeCall( cudaMemset(counterPtr, 0, sizeof(int)) );
+
+            const dim3 block(32, 8);
+            const dim3 grid(divUp(cols, block.x), divUp(rows, block.y));
+
+            cudaSafeCall( cudaFuncSetCacheConfig(GHT_Ballard_PosScale_findPosInHist, cudaFuncCachePreferL1) );
+
+            GHT_Ballard_PosScale_findPosInHist<<<grid, block>>>(hist, rows, cols, scaleRange, out, votes, maxSize, minScale, scaleStep, dp, threshold);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+
+            int totalCount;
+            cudaSafeCall( cudaMemcpy(&totalCount, counterPtr, sizeof(int), cudaMemcpyDeviceToHost) );
+
+            totalCount = ::min(totalCount, maxSize);
+
+            return totalCount;
+        }
+
+        ////////////////////////////////////////////////////////////////////////
+        // GHT_Ballard_PosRotation
+
+        __global__ void GHT_Ballard_PosRotation_calcHist(const unsigned int* coordList, const float* thetaList,
+                                                         PtrStep<short2> r_table, const int* r_sizes,
+                                                         PtrStepi hist, const int rows, const int cols,
+                                                         const float minAngle, const float angleStep, const int angleRange,
+                                                         const float idp, const float thetaScale)
+        {
+            const unsigned int coord = coordList[blockIdx.x];
+            float2 p;
+            p.x = (coord & 0xFFFF);
+            p.y = (coord >> 16) & 0xFFFF;
+
+            const float thetaVal = thetaList[blockIdx.x];
+
+            for (int a = threadIdx.x; a < angleRange; a += blockDim.x)
+            {
+                const float angle = (minAngle + a * angleStep) * (CV_PI_F / 180.0f);
+                float sinA, cosA;
+                sincosf(angle, &sinA, &cosA);
+
+                float theta = thetaVal - angle;
+                if (theta < 0)
+                    theta += 2.0f * CV_PI_F;
+
+                const int n = __float2int_rn(theta * thetaScale);
+
+                const short2* r_row = r_table.ptr(n);
+                const int r_row_size = r_sizes[n];
+
+                for (int j = 0; j < r_row_size; ++j)
+                {
+                    const float2 d = saturate_cast<float2>(r_row[j]);
+
+                    const float2 dr = make_float2(d.x * cosA - d.y * sinA, d.x * sinA + d.y * cosA);
+
+                    float2 c = make_float2(p.x - dr.x, p.y - dr.y);
+                    c.x *= idp;
+                    c.y *= idp;
+
+                    if (c.x >= 0 && c.x < cols && c.y >= 0 && c.y < rows)
+                        ::atomicAdd(hist.ptr((a + 1) * (rows + 2) + __float2int_rn(c.y + 1)) + __float2int_rn(c.x + 1), 1);
+                }
+            }
+        }
+
+        void GHT_Ballard_PosRotation_calcHist_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                                  PtrStepSz<short2> r_table, const int* r_sizes,
+                                                  PtrStepi hist, int rows, int cols,
+                                                  float minAngle, float angleStep, int angleRange,
+                                                  float dp, int levels)
+        {
+            const dim3 block(256);
+            const dim3 grid(pointsCount);
+
+            const float idp = 1.0f / dp;
+            const float thetaScale = levels / (2.0f * CV_PI_F);
+
+            GHT_Ballard_PosRotation_calcHist<<<grid, block>>>(coordList, thetaList,
+                                                              r_table, r_sizes,
+                                                              hist, rows, cols,
+                                                              minAngle, angleStep, angleRange,
+                                                              idp, thetaScale);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        __global__ void GHT_Ballard_PosRotation_findPosInHist(const PtrStepi hist, const int rows, const int cols, const int angleRange,
+                                                              float4* out, int3* votes, const int maxSize,
+                                                              const float minAngle, const float angleStep, const float dp, const int threshold)
+        {
+            const int x = blockIdx.x * blockDim.x + threadIdx.x;
+            const int y = blockIdx.y * blockDim.y + threadIdx.y;
+
+            if (x >= cols || y >= rows)
+                return;
+
+            for (int a = 0; a < angleRange; ++a)
+            {
+                const float angle = minAngle + a * angleStep;
+
+                const int prevAngleIdx = (a) * (rows + 2);
+                const int curAngleIdx = (a + 1) * (rows + 2);
+                const int nextAngleIdx = (a + 2) * (rows + 2);
+
+                const int curVotes = hist(curAngleIdx + y + 1, x + 1);
+
+                if (curVotes > threshold &&
+                    curVotes >  hist(curAngleIdx + y + 1, x) &&
+                    curVotes >= hist(curAngleIdx + y + 1, x + 2) &&
+                    curVotes >  hist(curAngleIdx + y, x + 1) &&
+                    curVotes >= hist(curAngleIdx + y + 2, x + 1) &&
+                    curVotes >  hist(prevAngleIdx + y + 1, x + 1) &&
+                    curVotes >= hist(nextAngleIdx + y + 1, x + 1))
+                {
+                    const int ind = ::atomicAdd(&g_counter, 1);
+
+                    if (ind < maxSize)
+                    {
+                        out[ind] = make_float4(x * dp, y * dp, 1.0f, angle);
+                        votes[ind] = make_int3(curVotes, 0, curVotes);
+                    }
+                }
+            }
+        }
+
+        int GHT_Ballard_PosRotation_findPosInHist_gpu(PtrStepi hist, int rows, int cols, int angleRange, float4* out, int3* votes, int maxSize,
+                                                      float minAngle, float angleStep, float dp, int threshold)
+        {
+            void* counterPtr;
+            cudaSafeCall( cudaGetSymbolAddress(&counterPtr, g_counter) );
+
+            cudaSafeCall( cudaMemset(counterPtr, 0, sizeof(int)) );
+
+            const dim3 block(32, 8);
+            const dim3 grid(divUp(cols, block.x), divUp(rows, block.y));
+
+            cudaSafeCall( cudaFuncSetCacheConfig(GHT_Ballard_PosRotation_findPosInHist, cudaFuncCachePreferL1) );
+
+            GHT_Ballard_PosRotation_findPosInHist<<<grid, block>>>(hist, rows, cols, angleRange, out, votes, maxSize, minAngle, angleStep, dp, threshold);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+
+            int totalCount;
+            cudaSafeCall( cudaMemcpy(&totalCount, counterPtr, sizeof(int), cudaMemcpyDeviceToHost) );
+
+            totalCount = ::min(totalCount, maxSize);
+
+            return totalCount;
+        }
+
+        ////////////////////////////////////////////////////////////////////////
+        // GHT_Guil_Full
+
+        struct FeatureTable
+        {
+            uchar* p1_pos_data;
+            size_t p1_pos_step;
+
+            uchar* p1_theta_data;
+            size_t p1_theta_step;
+
+            uchar* p2_pos_data;
+            size_t p2_pos_step;
+
+            uchar* d12_data;
+            size_t d12_step;
+
+            uchar* r1_data;
+            size_t r1_step;
+
+            uchar* r2_data;
+            size_t r2_step;
+        };
+
+        __constant__ FeatureTable c_templFeatures;
+        __constant__ FeatureTable c_imageFeatures;
+
+        void GHT_Guil_Full_setTemplFeatures(PtrStepb p1_pos, PtrStepb p1_theta, PtrStepb p2_pos, PtrStepb d12, PtrStepb r1, PtrStepb r2)
+        {
+            FeatureTable tbl;
+
+            tbl.p1_pos_data = p1_pos.data;
+            tbl.p1_pos_step = p1_pos.step;
+
+            tbl.p1_theta_data = p1_theta.data;
+            tbl.p1_theta_step = p1_theta.step;
+
+            tbl.p2_pos_data = p2_pos.data;
+            tbl.p2_pos_step = p2_pos.step;
+
+            tbl.d12_data = d12.data;
+            tbl.d12_step = d12.step;
+
+            tbl.r1_data = r1.data;
+            tbl.r1_step = r1.step;
+
+            tbl.r2_data = r2.data;
+            tbl.r2_step = r2.step;
+
+            cudaSafeCall( cudaMemcpyToSymbol(c_templFeatures, &tbl, sizeof(FeatureTable)) );
+        }
+        void GHT_Guil_Full_setImageFeatures(PtrStepb p1_pos, PtrStepb p1_theta, PtrStepb p2_pos, PtrStepb d12, PtrStepb r1, PtrStepb r2)
+        {
+            FeatureTable tbl;
+
+            tbl.p1_pos_data = p1_pos.data;
+            tbl.p1_pos_step = p1_pos.step;
+
+            tbl.p1_theta_data = p1_theta.data;
+            tbl.p1_theta_step = p1_theta.step;
+
+            tbl.p2_pos_data = p2_pos.data;
+            tbl.p2_pos_step = p2_pos.step;
+
+            tbl.d12_data = d12.data;
+            tbl.d12_step = d12.step;
+
+            tbl.r1_data = r1.data;
+            tbl.r1_step = r1.step;
+
+            tbl.r2_data = r2.data;
+            tbl.r2_step = r2.step;
+
+            cudaSafeCall( cudaMemcpyToSymbol(c_imageFeatures, &tbl, sizeof(FeatureTable)) );
+        }
+
+        struct TemplFeatureTable
+        {
+            static __device__ float2* p1_pos(int n)
+            {
+                return (float2*)(c_templFeatures.p1_pos_data + n * c_templFeatures.p1_pos_step);
+            }
+            static __device__ float* p1_theta(int n)
+            {
+                return (float*)(c_templFeatures.p1_theta_data + n * c_templFeatures.p1_theta_step);
+            }
+            static __device__ float2* p2_pos(int n)
+            {
+                return (float2*)(c_templFeatures.p2_pos_data + n * c_templFeatures.p2_pos_step);
+            }
+
+            static __device__ float* d12(int n)
+            {
+                return (float*)(c_templFeatures.d12_data + n * c_templFeatures.d12_step);
+            }
+
+            static __device__ float2* r1(int n)
+            {
+                return (float2*)(c_templFeatures.r1_data + n * c_templFeatures.r1_step);
+            }
+            static __device__ float2* r2(int n)
+            {
+                return (float2*)(c_templFeatures.r2_data + n * c_templFeatures.r2_step);
+            }
+        };
+        struct ImageFeatureTable
+        {
+            static __device__ float2* p1_pos(int n)
+            {
+                return (float2*)(c_imageFeatures.p1_pos_data + n * c_imageFeatures.p1_pos_step);
+            }
+            static __device__ float* p1_theta(int n)
+            {
+                return (float*)(c_imageFeatures.p1_theta_data + n * c_imageFeatures.p1_theta_step);
+            }
+            static __device__ float2* p2_pos(int n)
+            {
+                return (float2*)(c_imageFeatures.p2_pos_data + n * c_imageFeatures.p2_pos_step);
+            }
+
+            static __device__ float* d12(int n)
+            {
+                return (float*)(c_imageFeatures.d12_data + n * c_imageFeatures.d12_step);
+            }
+
+            static __device__ float2* r1(int n)
+            {
+                return (float2*)(c_imageFeatures.r1_data + n * c_imageFeatures.r1_step);
+            }
+            static __device__ float2* r2(int n)
+            {
+                return (float2*)(c_imageFeatures.r2_data + n * c_imageFeatures.r2_step);
+            }
+        };
+
+        __device__ float clampAngle(float a)
+        {
+            float res = a;
+
+            while (res > 2.0f * CV_PI_F)
+                res -= 2.0f * CV_PI_F;
+            while (res < 0.0f)
+                res += 2.0f * CV_PI_F;
+
+            return res;
+        }
+
+        __device__ bool angleEq(float a, float b, float eps)
+        {
+            return (::fabs(clampAngle(a - b)) <= eps);
+        }
+
+        template <class FT, bool isTempl>
+        __global__ void GHT_Guil_Full_buildFeatureList(const unsigned int* coordList, const float* thetaList, const int pointsCount,
+                                                       int* sizes, const int maxSize,
+                                                       const float xi, const float angleEpsilon, const float alphaScale,
+                                                       const float2 center, const float maxDist)
+        {
+            const float p1_theta = thetaList[blockIdx.x];
+            const unsigned int coord1 = coordList[blockIdx.x];
+            float2 p1_pos;
+            p1_pos.x = (coord1 & 0xFFFF);
+            p1_pos.y = (coord1 >> 16) & 0xFFFF;
+
+            for (int i = threadIdx.x; i < pointsCount; i += blockDim.x)
+            {
+                const float p2_theta = thetaList[i];
+                const unsigned int coord2 = coordList[i];
+                float2 p2_pos;
+                p2_pos.x = (coord2 & 0xFFFF);
+                p2_pos.y = (coord2 >> 16) & 0xFFFF;
+
+                if (angleEq(p1_theta - p2_theta, xi, angleEpsilon))
+                {
+                    const float2 d = p1_pos - p2_pos;
+
+                    float alpha12 = clampAngle(::atan2(d.y, d.x) - p1_theta);
+                    float d12 = ::sqrtf(d.x * d.x + d.y * d.y);
+
+                    if (d12 > maxDist)
+                        continue;
+
+                    float2 r1 = p1_pos - center;
+                    float2 r2 = p2_pos - center;
+
+                    const int n = __float2int_rn(alpha12 * alphaScale);
+
+                    const int ind = ::atomicAdd(sizes + n, 1);
+
+                    if (ind < maxSize)
+                    {
+                        if (!isTempl)
+                        {
+                            FT::p1_pos(n)[ind] = p1_pos;
+                            FT::p2_pos(n)[ind] = p2_pos;
+                        }
+
+                        FT::p1_theta(n)[ind] = p1_theta;
+
+                        FT::d12(n)[ind] = d12;
+
+                        if (isTempl)
+                        {
+                            FT::r1(n)[ind] = r1;
+                            FT::r2(n)[ind] = r2;
+                        }
+                    }
+                }
+            }
+        }
+
+        template <class FT, bool isTempl>
+        void GHT_Guil_Full_buildFeatureList_caller(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                                   int* sizes, int maxSize,
+                                                   float xi, float angleEpsilon, int levels,
+                                                   float2 center, float maxDist)
+        {
+            const dim3 block(256);
+            const dim3 grid(pointsCount);
+
+            const float alphaScale = levels / (2.0f * CV_PI_F);
+
+            GHT_Guil_Full_buildFeatureList<FT, isTempl><<<grid, block>>>(coordList, thetaList, pointsCount,
+                                                                         sizes, maxSize,
+                                                                         xi * (CV_PI_F / 180.0f), angleEpsilon * (CV_PI_F / 180.0f), alphaScale,
+                                                                         center, maxDist);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+
+            thrust::device_ptr<int> sizesPtr(sizes);
+            thrust::transform(sizesPtr, sizesPtr + levels + 1, sizesPtr, device::bind2nd(device::minimum<int>(), maxSize));
+        }
+
+        void GHT_Guil_Full_buildTemplFeatureList_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                                     int* sizes, int maxSize,
+                                                     float xi, float angleEpsilon, int levels,
+                                                     float2 center, float maxDist)
+        {
+            GHT_Guil_Full_buildFeatureList_caller<TemplFeatureTable, true>(coordList, thetaList, pointsCount,
+                                                                           sizes, maxSize,
+                                                                           xi, angleEpsilon, levels,
+                                                                           center, maxDist);
+        }
+        void GHT_Guil_Full_buildImageFeatureList_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                                     int* sizes, int maxSize,
+                                                     float xi, float angleEpsilon, int levels,
+                                                     float2 center, float maxDist)
+        {
+            GHT_Guil_Full_buildFeatureList_caller<ImageFeatureTable, false>(coordList, thetaList, pointsCount,
+                                                                            sizes, maxSize,
+                                                                            xi, angleEpsilon, levels,
+                                                                            center, maxDist);
+        }
+
+        __global__ void GHT_Guil_Full_calcOHist(const int* templSizes, const int* imageSizes, int* OHist,
+                                                const float minAngle, const float maxAngle, const float iAngleStep, const int angleRange)
+        {
+            extern __shared__ int s_OHist[];
+            for (int i = threadIdx.x; i <= angleRange; i += blockDim.x)
+                s_OHist[i] = 0;
+            __syncthreads();
+
+            const int tIdx = blockIdx.x;
+            const int level = blockIdx.y;
+
+            const int tSize = templSizes[level];
+
+            if (tIdx < tSize)
+            {
+                const int imSize = imageSizes[level];
+
+                const float t_p1_theta = TemplFeatureTable::p1_theta(level)[tIdx];
+
+                for (int i = threadIdx.x; i < imSize; i += blockDim.x)
+                {
+                    const float im_p1_theta = ImageFeatureTable::p1_theta(level)[i];
+
+                    const float angle = clampAngle(im_p1_theta - t_p1_theta);
+
+                    if (angle >= minAngle && angle <= maxAngle)
+                    {
+                        const int n = __float2int_rn((angle - minAngle) * iAngleStep);
+                        Emulation::smem::atomicAdd(&s_OHist[n], 1);
+                    }
+                }
+            }
+            __syncthreads();
+
+            for (int i = threadIdx.x; i <= angleRange; i += blockDim.x)
+                ::atomicAdd(OHist + i, s_OHist[i]);
+        }
+
+        void GHT_Guil_Full_calcOHist_gpu(const int* templSizes, const int* imageSizes, int* OHist,
+                                         float minAngle, float maxAngle, float angleStep, int angleRange,
+                                         int levels, int tMaxSize)
+        {
+            const dim3 block(256);
+            const dim3 grid(tMaxSize, levels + 1);
+
+            minAngle *= (CV_PI_F / 180.0f);
+            maxAngle *= (CV_PI_F / 180.0f);
+            angleStep *= (CV_PI_F / 180.0f);
+
+            const size_t smemSize = (angleRange + 1) * sizeof(float);
+
+            GHT_Guil_Full_calcOHist<<<grid, block, smemSize>>>(templSizes, imageSizes, OHist,
+                                                               minAngle, maxAngle, 1.0f / angleStep, angleRange);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        __global__ void GHT_Guil_Full_calcSHist(const int* templSizes, const int* imageSizes, int* SHist,
+                                                const float angle, const float angleEpsilon,
+                                                const float minScale, const float maxScale, const float iScaleStep, const int scaleRange)
+        {
+            extern __shared__ int s_SHist[];
+            for (int i = threadIdx.x; i <= scaleRange; i += blockDim.x)
+                s_SHist[i] = 0;
+            __syncthreads();
+
+            const int tIdx = blockIdx.x;
+            const int level = blockIdx.y;
+
+            const int tSize = templSizes[level];
+
+            if (tIdx < tSize)
+            {
+                const int imSize = imageSizes[level];
+
+                const float t_p1_theta = TemplFeatureTable::p1_theta(level)[tIdx] + angle;
+                const float t_d12 = TemplFeatureTable::d12(level)[tIdx] + angle;
+
+                for (int i = threadIdx.x; i < imSize; i += blockDim.x)
+                {
+                    const float im_p1_theta = ImageFeatureTable::p1_theta(level)[i];
+                    const float im_d12 = ImageFeatureTable::d12(level)[i];
+
+                    if (angleEq(im_p1_theta, t_p1_theta, angleEpsilon))
+                    {
+                        const float scale = im_d12 / t_d12;
+
+                        if (scale >= minScale && scale <= maxScale)
+                        {
+                            const int s = __float2int_rn((scale - minScale) * iScaleStep);
+                            Emulation::smem::atomicAdd(&s_SHist[s], 1);
+                        }
+                    }
+                }
+            }
+            __syncthreads();
+
+            for (int i = threadIdx.x; i <= scaleRange; i += blockDim.x)
+                ::atomicAdd(SHist + i, s_SHist[i]);
+        }
+
+        void GHT_Guil_Full_calcSHist_gpu(const int* templSizes, const int* imageSizes, int* SHist,
+                                         float angle, float angleEpsilon,
+                                         float minScale, float maxScale, float iScaleStep, int scaleRange,
+                                         int levels, int tMaxSize)
+        {
+            const dim3 block(256);
+            const dim3 grid(tMaxSize, levels + 1);
+
+            angle *= (CV_PI_F / 180.0f);
+            angleEpsilon *= (CV_PI_F / 180.0f);
+
+            const size_t smemSize = (scaleRange + 1) * sizeof(float);
+
+            GHT_Guil_Full_calcSHist<<<grid, block, smemSize>>>(templSizes, imageSizes, SHist,
+                                                               angle, angleEpsilon,
+                                                               minScale, maxScale, iScaleStep, scaleRange);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        __global__ void GHT_Guil_Full_calcPHist(const int* templSizes, const int* imageSizes, PtrStepSzi PHist,
+                                                const float angle, const float sinVal, const float cosVal, const float angleEpsilon, const float scale,
+                                                const float idp)
+        {
+            const int tIdx = blockIdx.x;
+            const int level = blockIdx.y;
+
+            const int tSize = templSizes[level];
+
+            if (tIdx < tSize)
+            {
+                const int imSize = imageSizes[level];
+
+                const float t_p1_theta = TemplFeatureTable::p1_theta(level)[tIdx] + angle;
+
+                float2 r1 = TemplFeatureTable::r1(level)[tIdx];
+                float2 r2 = TemplFeatureTable::r2(level)[tIdx];
+
+                r1 = r1 * scale;
+                r2 = r2 * scale;
+
+                r1 = make_float2(cosVal * r1.x - sinVal * r1.y, sinVal * r1.x + cosVal * r1.y);
+                r2 = make_float2(cosVal * r2.x - sinVal * r2.y, sinVal * r2.x + cosVal * r2.y);
+
+                for (int i = threadIdx.x; i < imSize; i += blockDim.x)
+                {
+                    const float im_p1_theta = ImageFeatureTable::p1_theta(level)[i];
+
+                    const float2 im_p1_pos = ImageFeatureTable::p1_pos(level)[i];
+                    const float2 im_p2_pos = ImageFeatureTable::p2_pos(level)[i];
+
+                    if (angleEq(im_p1_theta, t_p1_theta, angleEpsilon))
+                    {
+                        float2 c1, c2;
+
+                        c1 = im_p1_pos - r1;
+                        c1 = c1 * idp;
+
+                        c2 = im_p2_pos - r2;
+                        c2 = c2 * idp;
+
+                        if (::fabs(c1.x - c2.x) > 1 || ::fabs(c1.y - c2.y) > 1)
+                            continue;
+
+                        if (c1.y >= 0 && c1.y < PHist.rows - 2 && c1.x >= 0 && c1.x < PHist.cols - 2)
+                            ::atomicAdd(PHist.ptr(__float2int_rn(c1.y) + 1) + __float2int_rn(c1.x) + 1, 1);
+                    }
+                }
+            }
+        }
+
+        void GHT_Guil_Full_calcPHist_gpu(const int* templSizes, const int* imageSizes, PtrStepSzi PHist,
+                                         float angle, float angleEpsilon, float scale,
+                                         float dp,
+                                         int levels, int tMaxSize)
+        {
+            const dim3 block(256);
+            const dim3 grid(tMaxSize, levels + 1);
+
+            angle *= (CV_PI_F / 180.0f);
+            angleEpsilon *= (CV_PI_F / 180.0f);
+
+            const float sinVal = ::sinf(angle);
+            const float cosVal = ::cosf(angle);
+
+            cudaSafeCall( cudaFuncSetCacheConfig(GHT_Guil_Full_calcPHist, cudaFuncCachePreferL1) );
+
+            GHT_Guil_Full_calcPHist<<<grid, block>>>(templSizes, imageSizes, PHist,
+                                                     angle, sinVal, cosVal, angleEpsilon, scale,
+                                                     1.0f / dp);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        __global__ void GHT_Guil_Full_findPosInHist(const PtrStepSzi hist, float4* out, int3* votes, const int maxSize,
+                                                    const float angle, const int angleVotes, const float scale, const int scaleVotes,
+                                                    const float dp, const int threshold)
+        {
+            const int x = blockIdx.x * blockDim.x + threadIdx.x;
+            const int y = blockIdx.y * blockDim.y + threadIdx.y;
+
+            if (x >= hist.cols - 2 || y >= hist.rows - 2)
+                return;
+
+            const int curVotes = hist(y + 1, x + 1);
+
+            if (curVotes > threshold &&
+                curVotes >  hist(y + 1, x) &&
+                curVotes >= hist(y + 1, x + 2) &&
+                curVotes >  hist(y, x + 1) &&
+                curVotes >= hist(y + 2, x + 1))
+            {
+                const int ind = ::atomicAdd(&g_counter, 1);
+
+                if (ind < maxSize)
+                {
+                    out[ind] = make_float4(x * dp, y * dp, scale, angle);
+                    votes[ind] = make_int3(curVotes, scaleVotes, angleVotes);
+                }
+            }
+        }
+
+        int GHT_Guil_Full_findPosInHist_gpu(PtrStepSzi hist, float4* out, int3* votes, int curSize, int maxSize,
+                                             float angle, int angleVotes, float scale, int scaleVotes,
+                                             float dp, int threshold)
+        {
+            void* counterPtr;
+            cudaSafeCall( cudaGetSymbolAddress(&counterPtr, g_counter) );
+
+            cudaSafeCall( cudaMemcpy(counterPtr, &curSize, sizeof(int), cudaMemcpyHostToDevice) );
+
+            const dim3 block(32, 8);
+            const dim3 grid(divUp(hist.cols - 2, block.x), divUp(hist.rows - 2, block.y));
+
+            cudaSafeCall( cudaFuncSetCacheConfig(GHT_Guil_Full_findPosInHist, cudaFuncCachePreferL1) );
+
+            GHT_Guil_Full_findPosInHist<<<grid, block>>>(hist, out, votes, maxSize,
+                                                         angle, angleVotes, scale, scaleVotes,
+                                                         dp, threshold);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+
+            int totalCount;
+            cudaSafeCall( cudaMemcpy(&totalCount, counterPtr, sizeof(int), cudaMemcpyDeviceToHost) );
+
+            totalCount = ::min(totalCount, maxSize);
+
+            return totalCount;
+        }
     }
 }}}

modules/gpu/src/hough.cpp

@@ -42,6 +42,10 @@
 
 #include "precomp.hpp"
 
+using namespace std;
+using namespace cv;
+using namespace cv::gpu;
+
 #if !defined (HAVE_CUDA)
 
 void cv::gpu::HoughLines(const GpuMat&, GpuMat&, float, float, int, bool, int) { throw_nogpu(); }
@@ -52,6 +56,15 @@ void cv::gpu::HoughCircles(const GpuMat&, GpuMat&, int, float, float, int, int,
 void cv::gpu::HoughCircles(const GpuMat&, GpuMat&, HoughCirclesBuf&, int, float, float, int, int, int, int, int) { throw_nogpu(); }
 void cv::gpu::HoughCirclesDownload(const GpuMat&, OutputArray) { throw_nogpu(); }
 
+Ptr<GeneralizedHough_GPU> cv::gpu::GeneralizedHough_GPU::create(int) { throw_nogpu(); return Ptr<GeneralizedHough_GPU>(); }
+cv::gpu::GeneralizedHough_GPU::~GeneralizedHough_GPU() {}
+void cv::gpu::GeneralizedHough_GPU::setTemplate(const GpuMat&, int, Point) { throw_nogpu(); }
+void cv::gpu::GeneralizedHough_GPU::setTemplate(const GpuMat&, const GpuMat&, const GpuMat&, Point) { throw_nogpu(); }
+void cv::gpu::GeneralizedHough_GPU::detect(const GpuMat&, GpuMat&, int) { throw_nogpu(); }
+void cv::gpu::GeneralizedHough_GPU::detect(const GpuMat&, const GpuMat&, const GpuMat&, GpuMat&) { throw_nogpu(); }
+void cv::gpu::GeneralizedHough_GPU::download(const GpuMat&, OutputArray, OutputArray) { throw_nogpu(); }
+void cv::gpu::GeneralizedHough_GPU::release() {}
+
 #else /* !defined (HAVE_CUDA) */
 
 namespace cv { namespace gpu { namespace device
@@ -59,20 +72,21 @@ namespace cv { namespace gpu { namespace device
     namespace hough
     {
         int buildPointList_gpu(PtrStepSzb src, unsigned int* list);
-
-        void linesAccum_gpu(const unsigned int* list, int count, PtrStepSzi accum, float rho, float theta, size_t sharedMemPerBlock, bool has20);
-        int linesGetResult_gpu(PtrStepSzi accum, float2* out, int* votes, int maxSize, float rho, float theta, int threshold, bool doSort);
-
-        void circlesAccumCenters_gpu(const unsigned int* list, int count, PtrStepi dx, PtrStepi dy, PtrStepSzi accum, int minRadius, int maxRadius, float idp);
-        int buildCentersList_gpu(PtrStepSzi accum, unsigned int* centers, int threshold);
-        int circlesAccumRadius_gpu(const unsigned int* centers, int centersCount, const unsigned int* list, int count,
-                                   float3* circles, int maxCircles, float dp, int minRadius, int maxRadius, int threshold, bool has20);
     }
 }}}
 
 //////////////////////////////////////////////////////////
 // HoughLines
 
+namespace cv { namespace gpu { namespace device
+{
+    namespace hough
+    {
+        void linesAccum_gpu(const unsigned int* list, int count, PtrStepSzi accum, float rho, float theta, size_t sharedMemPerBlock, bool has20);
+        int linesGetResult_gpu(PtrStepSzi accum, float2* out, int* votes, int maxSize, float rho, float theta, int threshold, bool doSort);
+    }
+}}}
+
 void cv::gpu::HoughLines(const GpuMat& src, GpuMat& lines, float rho, float theta, int threshold, bool doSort, int maxLines)
 {
     HoughLinesBuf buf;
@@ -144,6 +158,17 @@ void cv::gpu::HoughLinesDownload(const GpuMat& d_lines, OutputArray h_lines_, Ou
 //////////////////////////////////////////////////////////
 // HoughCircles
 
+namespace cv { namespace gpu { namespace device
+{
+    namespace hough
+    {
+        void circlesAccumCenters_gpu(const unsigned int* list, int count, PtrStepi dx, PtrStepi dy, PtrStepSzi accum, int minRadius, int maxRadius, float idp);
+        int buildCentersList_gpu(PtrStepSzi accum, unsigned int* centers, int threshold);
+        int circlesAccumRadius_gpu(const unsigned int* centers, int centersCount, const unsigned int* list, int count,
+                                   float3* circles, int maxCircles, float dp, int minRadius, int maxRadius, int threshold, bool has20);
+    }
+}}}
+
 void cv::gpu::HoughCircles(const GpuMat& src, GpuMat& circles, int method, float dp, float minDist, int cannyThreshold, int votesThreshold, int minRadius, int maxRadius, int maxCircles)
 {
     HoughCirclesBuf buf;
@@ -209,7 +234,7 @@ void cv::gpu::HoughCircles(const GpuMat& src, GpuMat& circles, HoughCirclesBuf&
 
         std::vector< std::vector<ushort2> > grid(gridWidth * gridHeight);
 
-        minDist *= minDist;
+        const float minDist2 = minDist * minDist;
 
         for (int i = 0; i < centersCount; ++i)
         {
@@ -242,7 +267,7 @@ void cv::gpu::HoughCircles(const GpuMat& src, GpuMat& circles, HoughCirclesBuf&
                         float dx = (float)(p.x - m[j].x);
                         float dy = (float)(p.y - m[j].y);
 
-                        if (dx * dx + dy * dy < minDist)
+                        if (dx * dx + dy * dy < minDist2)
                         {
                             good = false;
                             goto break_out;
@@ -292,4 +317,1056 @@ void cv::gpu::HoughCirclesDownload(const GpuMat& d_circles, cv::OutputArray h_ci
     d_circles.download(h_circles);
 }
 
+//////////////////////////////////////////////////////////
+// GeneralizedHough
+
+namespace cv { namespace gpu { namespace device
+{
+    namespace hough
+    {
+        template <typename T>
+        int buildEdgePointList_gpu(PtrStepSzb edges, PtrStepSzb dx, PtrStepSzb dy, unsigned int* coordList, float* thetaList);
+        void buildRTable_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                             PtrStepSz<short2> r_table, int* r_sizes,
+                             short2 templCenter, int levels);
+
+        void GHT_Ballard_Pos_calcHist_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                          PtrStepSz<short2> r_table, const int* r_sizes,
+                                          PtrStepSzi hist,
+                                          float dp, int levels);
+        int GHT_Ballard_Pos_findPosInHist_gpu(PtrStepSzi hist, float4* out, int3* votes, int maxSize, float dp, int threshold);
+
+        void GHT_Ballard_PosScale_calcHist_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                               PtrStepSz<short2> r_table, const int* r_sizes,
+                                               PtrStepi hist, int rows, int cols,
+                                               float minScale, float scaleStep, int scaleRange,
+                                               float dp, int levels);
+        int GHT_Ballard_PosScale_findPosInHist_gpu(PtrStepi hist, int rows, int cols, int scaleRange, float4* out, int3* votes, int maxSize,
+                                                   float minScale, float scaleStep, float dp, int threshold);
+
+        void GHT_Ballard_PosRotation_calcHist_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                                  PtrStepSz<short2> r_table, const int* r_sizes,
+                                                  PtrStepi hist, int rows, int cols,
+                                                  float minAngle, float angleStep, int angleRange,
+                                                  float dp, int levels);
+        int GHT_Ballard_PosRotation_findPosInHist_gpu(PtrStepi hist, int rows, int cols, int angleRange, float4* out, int3* votes, int maxSize,
+                                                      float minAngle, float angleStep, float dp, int threshold);
+
+        void GHT_Guil_Full_setTemplFeatures(PtrStepb p1_pos, PtrStepb p1_theta, PtrStepb p2_pos, PtrStepb d12, PtrStepb r1, PtrStepb r2);
+        void GHT_Guil_Full_setImageFeatures(PtrStepb p1_pos, PtrStepb p1_theta, PtrStepb p2_pos, PtrStepb d12, PtrStepb r1, PtrStepb r2);
+        void GHT_Guil_Full_buildTemplFeatureList_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                                     int* sizes, int maxSize,
+                                                     float xi, float angleEpsilon, int levels,
+                                                     float2 center, float maxDist);
+        void GHT_Guil_Full_buildImageFeatureList_gpu(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                                     int* sizes, int maxSize,
+                                                     float xi, float angleEpsilon, int levels,
+                                                     float2 center, float maxDist);
+        void GHT_Guil_Full_calcOHist_gpu(const int* templSizes, const int* imageSizes, int* OHist,
+                                         float minAngle, float maxAngle, float angleStep, int angleRange,
+                                         int levels, int tMaxSize);
+        void GHT_Guil_Full_calcSHist_gpu(const int* templSizes, const int* imageSizes, int* SHist,
+                                         float angle, float angleEpsilon,
+                                         float minScale, float maxScale, float iScaleStep, int scaleRange,
+                                         int levels, int tMaxSize);
+        void GHT_Guil_Full_calcPHist_gpu(const int* templSizes, const int* imageSizes, PtrStepSzi PHist,
+                                         float angle, float angleEpsilon, float scale,
+                                         float dp,
+                                         int levels, int tMaxSize);
+        int GHT_Guil_Full_findPosInHist_gpu(PtrStepSzi hist, float4* out, int3* votes, int curSize, int maxSize,
+                                             float angle, int angleVotes, float scale, int scaleVotes,
+                                             float dp, int threshold);
+    }
+}}}
+
+namespace
+{
+    /////////////////////////////////////
+    // Common
+
+    template <typename T, class A> void releaseVector(vector<T, A>& v)
+    {
+        vector<T, A> empty;
+        empty.swap(v);
+    }
+
+    class GHT_Pos : public GeneralizedHough_GPU
+    {
+    public:
+        GHT_Pos();
+
+    protected:
+        void setTemplateImpl(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, Point templCenter);
+        void detectImpl(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, GpuMat& positions);
+        void releaseImpl();
+
+        virtual void processTempl() = 0;
+        virtual void processImage() = 0;
+
+        void buildEdgePointList(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy);
+        void filterMinDist();
+        void convertTo(GpuMat& positions);
+
+        int maxSize;
+        double minDist;
+
+        Size templSize;
+        Point templCenter;
+        GpuMat templEdges;
+        GpuMat templDx;
+        GpuMat templDy;
+
+        Size imageSize;
+        GpuMat imageEdges;
+        GpuMat imageDx;
+        GpuMat imageDy;
+
+        GpuMat edgePointList;
+
+        GpuMat outBuf;
+        int posCount;
+
+        vector<float4> oldPosBuf;
+        vector<int3> oldVoteBuf;
+        vector<float4> newPosBuf;
+        vector<int3> newVoteBuf;
+        vector<int> indexies;
+    };
+
+    GHT_Pos::GHT_Pos()
+    {
+        maxSize = 10000;
+        minDist = 1.0;
+    }
+
+    void GHT_Pos::setTemplateImpl(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, Point templCenter_)
+    {
+        templSize = edges.size();
+        templCenter = templCenter_;
+
+        ensureSizeIsEnough(templSize, edges.type(), templEdges);
+        ensureSizeIsEnough(templSize, dx.type(), templDx);
+        ensureSizeIsEnough(templSize, dy.type(), templDy);
+
+        edges.copyTo(templEdges);
+        dx.copyTo(templDx);
+        dy.copyTo(templDy);
+
+        processTempl();
+    }
+
+    void GHT_Pos::detectImpl(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, GpuMat& positions)
+    {
+        imageSize = edges.size();
+
+        ensureSizeIsEnough(imageSize, edges.type(), imageEdges);
+        ensureSizeIsEnough(imageSize, dx.type(), imageDx);
+        ensureSizeIsEnough(imageSize, dy.type(), imageDy);
+
+        edges.copyTo(imageEdges);
+        dx.copyTo(imageDx);
+        dy.copyTo(imageDy);
+
+        posCount = 0;
+
+        processImage();
+
+        if (posCount == 0)
+            positions.release();
+        else
+        {
+            if (minDist > 1)
+                filterMinDist();
+            convertTo(positions);
+        }
+    }
+
+    void GHT_Pos::releaseImpl()
+    {
+        templSize = Size();
+        templCenter = Point(-1, -1);
+        templEdges.release();
+        templDx.release();
+        templDy.release();
+
+        imageSize = Size();
+        imageEdges.release();
+        imageDx.release();
+        imageDy.release();
+
+        edgePointList.release();
+
+        outBuf.release();
+        posCount = 0;
+
+        releaseVector(oldPosBuf);
+        releaseVector(oldVoteBuf);
+        releaseVector(newPosBuf);
+        releaseVector(newVoteBuf);
+        releaseVector(indexies);
+    }
+
+    void GHT_Pos::buildEdgePointList(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy)
+    {
+        using namespace cv::gpu::device::hough;
+
+        typedef int (*func_t)(PtrStepSzb edges, PtrStepSzb dx, PtrStepSzb dy, unsigned int* coordList, float* thetaList);
+        static const func_t funcs[] =
+        {
+            0,
+            0,
+            0,
+            buildEdgePointList_gpu<short>,
+            buildEdgePointList_gpu<int>,
+            buildEdgePointList_gpu<float>,
+            0
+        };
+
+        CV_Assert(edges.type() == CV_8UC1);
+        CV_Assert(dx.size() == edges.size());
+        CV_Assert(dy.type() == dx.type() && dy.size() == edges.size());
+
+        const func_t func = funcs[dx.depth()];
+        CV_Assert(func != 0);
+
+        edgePointList.cols = edgePointList.step / sizeof(int);
+        ensureSizeIsEnough(2, edges.size().area(), CV_32SC1, edgePointList);
+
+        edgePointList.cols = func(edges, dx, dy, edgePointList.ptr<unsigned int>(0), edgePointList.ptr<float>(1));
+    }
+
+    #define votes_cmp_gt(l1, l2) (aux[l1].x > aux[l2].x)
+    static CV_IMPLEMENT_QSORT_EX( sortIndexies, int, votes_cmp_gt, const int3* )
+
+    void GHT_Pos::filterMinDist()
+    {
+        oldPosBuf.resize(posCount);
+        oldVoteBuf.resize(posCount);
+
+        cudaSafeCall( cudaMemcpy(&oldPosBuf[0], outBuf.ptr(0), posCount * sizeof(float4), cudaMemcpyDeviceToHost) );
+        cudaSafeCall( cudaMemcpy(&oldVoteBuf[0], outBuf.ptr(1), posCount * sizeof(int3), cudaMemcpyDeviceToHost) );
+
+        indexies.resize(posCount);
+        for (int i = 0; i < posCount; ++i)
+            indexies[i] = i;
+        sortIndexies(&indexies[0], posCount, &oldVoteBuf[0]);
+
+        newPosBuf.clear();
+        newVoteBuf.clear();
+        newPosBuf.reserve(posCount);
+        newVoteBuf.reserve(posCount);
+
+        const int cellSize = cvRound(minDist);
+        const int gridWidth = (imageSize.width + cellSize - 1) / cellSize;
+        const int gridHeight = (imageSize.height + cellSize - 1) / cellSize;
+
+        vector< vector<Point2f> > grid(gridWidth * gridHeight);
+
+        const double minDist2 = minDist * minDist;
+
+        for (int i = 0; i < posCount; ++i)
+        {
+            const int ind = indexies[i];
+
+            Point2f p(oldPosBuf[ind].x, oldPosBuf[ind].y);
+
+            bool good = true;
+
+            const int xCell = static_cast<int>(p.x / cellSize);
+            const int yCell = static_cast<int>(p.y / cellSize);
+
+            int x1 = xCell - 1;
+            int y1 = yCell - 1;
+            int x2 = xCell + 1;
+            int y2 = yCell + 1;
+
+            // boundary check
+            x1 = std::max(0, x1);
+            y1 = std::max(0, y1);
+            x2 = std::min(gridWidth - 1, x2);
+            y2 = std::min(gridHeight - 1, y2);
+
+            for (int yy = y1; yy <= y2; ++yy)
+            {
+                for (int xx = x1; xx <= x2; ++xx)
+                {
+                    const vector<Point2f>& m = grid[yy * gridWidth + xx];
+
+                    for(size_t j = 0; j < m.size(); ++j)
+                    {
+                        const Point2f d = p - m[j];
+
+                        if (d.ddot(d) < minDist2)
+                        {
+                            good = false;
+                            goto break_out;
+                        }
+                    }
+                }
+            }
+
+            break_out:
+
+            if(good)
+            {
+                grid[yCell * gridWidth + xCell].push_back(p);
+
+                newPosBuf.push_back(oldPosBuf[ind]);
+                newVoteBuf.push_back(oldVoteBuf[ind]);
+            }
+        }
+
+        posCount = static_cast<int>(newPosBuf.size());
+        cudaSafeCall( cudaMemcpy(outBuf.ptr(0), &newPosBuf[0], posCount * sizeof(float4), cudaMemcpyHostToDevice) );
+        cudaSafeCall( cudaMemcpy(outBuf.ptr(1), &newVoteBuf[0], posCount * sizeof(int3), cudaMemcpyHostToDevice) );
+    }
+
+    void GHT_Pos::convertTo(GpuMat& positions)
+    {
+        ensureSizeIsEnough(2, posCount, CV_32FC4, positions);
+        GpuMat(2, posCount, CV_32FC4, outBuf.data, outBuf.step).copyTo(positions);
+    }
+
+    /////////////////////////////////////
+    // POSITION Ballard
+
+    class GHT_Ballard_Pos : public GHT_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Ballard_Pos();
+
+    protected:
+        void releaseImpl();
+
+        void processTempl();
+        void processImage();
+
+        virtual void calcHist();
+        virtual void findPosInHist();
+
+        int levels;
+        int votesThreshold;
+        double dp;
+
+        GpuMat r_table;
+        GpuMat r_sizes;
+
+        GpuMat hist;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Ballard_Pos, "GeneralizedHough_GPU.POSITION",
+                      obj.info()->addParam(obj, "maxSize", obj.maxSize, false, 0, 0,
+                                           "Maximal size of inner buffers.");
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "R-Table levels.");
+                      obj.info()->addParam(obj, "votesThreshold", obj.votesThreshold, false, 0, 0,
+                                           "The accumulator threshold for the template centers at the detection stage. The smaller it is, the more false positions may be detected.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution."));
+
+    GHT_Ballard_Pos::GHT_Ballard_Pos()
+    {
+        levels = 360;
+        votesThreshold = 100;
+        dp = 1.0;
+    }
+
+    void GHT_Ballard_Pos::releaseImpl()
+    {
+        GHT_Pos::releaseImpl();
+
+        r_table.release();
+        r_sizes.release();
+
+        hist.release();
+    }
+
+    void GHT_Ballard_Pos::processTempl()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(levels > 0);
+
+        buildEdgePointList(templEdges, templDx, templDy);
+
+        ensureSizeIsEnough(levels + 1, maxSize, CV_16SC2, r_table);
+        ensureSizeIsEnough(1, levels + 1, CV_32SC1, r_sizes);
+        r_sizes.setTo(Scalar::all(0));
+
+        if (edgePointList.cols > 0)
+        {
+            buildRTable_gpu(edgePointList.ptr<unsigned int>(0), edgePointList.ptr<float>(1), edgePointList.cols,
+                            r_table, r_sizes.ptr<int>(), make_short2(templCenter.x, templCenter.y), levels);
+            min(r_sizes, maxSize, r_sizes);
+        }
+    }
+
+    void GHT_Ballard_Pos::processImage()
+    {
+        calcHist();
+        findPosInHist();
+    }
+
+    void GHT_Ballard_Pos::calcHist()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(levels > 0 && r_table.rows == (levels + 1) && r_sizes.cols == (levels + 1));
+        CV_Assert(dp > 0.0);
+
+        const double idp = 1.0 / dp;
+
+        buildEdgePointList(imageEdges, imageDx, imageDy);
+
+        ensureSizeIsEnough(cvCeil(imageSize.height * idp) + 2, cvCeil(imageSize.width * idp) + 2, CV_32SC1, hist);
+        hist.setTo(Scalar::all(0));
+
+        if (edgePointList.cols > 0)
+        {
+            GHT_Ballard_Pos_calcHist_gpu(edgePointList.ptr<unsigned int>(0), edgePointList.ptr<float>(1), edgePointList.cols,
+                                         r_table, r_sizes.ptr<int>(),
+                                         hist,
+                                         dp, levels);
+        }
+    }
+
+    void GHT_Ballard_Pos::findPosInHist()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(votesThreshold > 0);
+
+        ensureSizeIsEnough(2, maxSize, CV_32FC4, outBuf);
+
+        posCount = GHT_Ballard_Pos_findPosInHist_gpu(hist, outBuf.ptr<float4>(0), outBuf.ptr<int3>(1), maxSize, dp, votesThreshold);
+    }
+
+    /////////////////////////////////////
+    // POSITION & SCALE
+
+    class GHT_Ballard_PosScale : public GHT_Ballard_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Ballard_PosScale();
+
+    protected:
+        void calcHist();
+        void findPosInHist();
+
+        double minScale;
+        double maxScale;
+        double scaleStep;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Ballard_PosScale, "GeneralizedHough_GPU.POSITION_SCALE",
+                      obj.info()->addParam(obj, "maxSize", obj.maxSize, false, 0, 0,
+                                           "Maximal size of inner buffers.");
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "R-Table levels.");
+                      obj.info()->addParam(obj, "votesThreshold", obj.votesThreshold, false, 0, 0,
+                                           "The accumulator threshold for the template centers at the detection stage. The smaller it is, the more false positions may be detected.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution.");
+                      obj.info()->addParam(obj, "minScale", obj.minScale, false, 0, 0,
+                                           "Minimal scale to detect.");
+                      obj.info()->addParam(obj, "maxScale", obj.maxScale, false, 0, 0,
+                                           "Maximal scale to detect.");
+                      obj.info()->addParam(obj, "scaleStep", obj.scaleStep, false, 0, 0,
+                                           "Scale step."));
+
+    GHT_Ballard_PosScale::GHT_Ballard_PosScale()
+    {
+        minScale = 0.5;
+        maxScale = 2.0;
+        scaleStep = 0.05;
+    }
+
+    void GHT_Ballard_PosScale::calcHist()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(levels > 0 && r_table.rows == (levels + 1) && r_sizes.cols == (levels + 1));
+        CV_Assert(dp > 0.0);
+        CV_Assert(minScale > 0.0 && minScale < maxScale);
+        CV_Assert(scaleStep > 0.0);
+
+        const double idp = 1.0 / dp;
+        const int scaleRange = cvCeil((maxScale - minScale) / scaleStep);
+        const int rows = cvCeil(imageSize.height * idp);
+        const int cols = cvCeil(imageSize.width * idp);
+
+        buildEdgePointList(imageEdges, imageDx, imageDy);
+
+        ensureSizeIsEnough((scaleRange + 2) * (rows + 2), cols + 2, CV_32SC1, hist);
+        hist.setTo(Scalar::all(0));
+
+        if (edgePointList.cols > 0)
+        {
+            GHT_Ballard_PosScale_calcHist_gpu(edgePointList.ptr<unsigned int>(0), edgePointList.ptr<float>(1), edgePointList.cols,
+                                              r_table, r_sizes.ptr<int>(),
+                                              hist, rows, cols,
+                                              minScale, scaleStep, scaleRange, dp, levels);
+        }
+    }
+
+    void GHT_Ballard_PosScale::findPosInHist()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(votesThreshold > 0);
+
+        const double idp = 1.0 / dp;
+        const int scaleRange = cvCeil((maxScale - minScale) / scaleStep);
+        const int rows = cvCeil(imageSize.height * idp);
+        const int cols = cvCeil(imageSize.width * idp);
+
+        ensureSizeIsEnough(2, maxSize, CV_32FC4, outBuf);
+
+        posCount =  GHT_Ballard_PosScale_findPosInHist_gpu(hist, rows, cols, scaleRange, outBuf.ptr<float4>(0), outBuf.ptr<int3>(1), maxSize, minScale, scaleStep, dp, votesThreshold);
+    }
+
+    /////////////////////////////////////
+    // POSITION & Rotation
+
+    class GHT_Ballard_PosRotation : public GHT_Ballard_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Ballard_PosRotation();
+
+    protected:
+        void calcHist();
+        void findPosInHist();
+
+        double minAngle;
+        double maxAngle;
+        double angleStep;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Ballard_PosRotation, "GeneralizedHough_GPU.POSITION_ROTATION",
+                      obj.info()->addParam(obj, "maxSize", obj.maxSize, false, 0, 0,
+                                           "Maximal size of inner buffers.");
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "R-Table levels.");
+                      obj.info()->addParam(obj, "votesThreshold", obj.votesThreshold, false, 0, 0,
+                                           "The accumulator threshold for the template centers at the detection stage. The smaller it is, the more false positions may be detected.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution.");
+                      obj.info()->addParam(obj, "minAngle", obj.minAngle, false, 0, 0,
+                                           "Minimal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "maxAngle", obj.maxAngle, false, 0, 0,
+                                           "Maximal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "angleStep", obj.angleStep, false, 0, 0,
+                                           "Angle step in degrees."));
+
+    GHT_Ballard_PosRotation::GHT_Ballard_PosRotation()
+    {
+        minAngle = 0.0;
+        maxAngle = 360.0;
+        angleStep = 1.0;
+    }
+
+    void GHT_Ballard_PosRotation::calcHist()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(levels > 0 && r_table.rows == (levels + 1) && r_sizes.cols == (levels + 1));
+        CV_Assert(dp > 0.0);
+        CV_Assert(minAngle >= 0.0 && minAngle < maxAngle && maxAngle <= 360.0);
+        CV_Assert(angleStep > 0.0 && angleStep < 360.0);
+
+        const double idp = 1.0 / dp;
+        const int angleRange = cvCeil((maxAngle - minAngle) / angleStep);
+        const int rows = cvCeil(imageSize.height * idp);
+        const int cols = cvCeil(imageSize.width * idp);
+
+        buildEdgePointList(imageEdges, imageDx, imageDy);
+
+        ensureSizeIsEnough((angleRange + 2) * (rows + 2), cols + 2, CV_32SC1, hist);
+        hist.setTo(Scalar::all(0));
+
+        if (edgePointList.cols > 0)
+        {
+            GHT_Ballard_PosRotation_calcHist_gpu(edgePointList.ptr<unsigned int>(0), edgePointList.ptr<float>(1), edgePointList.cols,
+                                                 r_table, r_sizes.ptr<int>(),
+                                                 hist, rows, cols,
+                                                 minAngle, angleStep, angleRange, dp, levels);
+        }
+    }
+
+    void GHT_Ballard_PosRotation::findPosInHist()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(votesThreshold > 0);
+
+        const double idp = 1.0 / dp;
+        const int angleRange = cvCeil((maxAngle - minAngle) / angleStep);
+        const int rows = cvCeil(imageSize.height * idp);
+        const int cols = cvCeil(imageSize.width * idp);
+
+        ensureSizeIsEnough(2, maxSize, CV_32FC4, outBuf);
+
+        posCount = GHT_Ballard_PosRotation_findPosInHist_gpu(hist, rows, cols, angleRange, outBuf.ptr<float4>(0), outBuf.ptr<int3>(1), maxSize, minAngle, angleStep, dp, votesThreshold);
+    }
+
+    /////////////////////////////////////////
+    // POSITION & SCALE & ROTATION
+
+    double toRad(double a)
+    {
+        return a * CV_PI / 180.0;
+    }
+
+    double clampAngle(double a)
+    {
+        double res = a;
+
+        while (res > 360.0)
+            res -= 360.0;
+        while (res < 0)
+            res += 360.0;
+
+        return res;
+    }
+
+    bool angleEq(double a, double b, double eps = 1.0)
+    {
+        return (fabs(clampAngle(a - b)) <= eps);
+    }
+
+    class GHT_Guil_Full : public GHT_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Guil_Full();
+
+    protected:
+        void releaseImpl();
+
+        void processTempl();
+        void processImage();
+
+        struct Feature
+        {
+            GpuMat p1_pos;
+            GpuMat p1_theta;
+            GpuMat p2_pos;
+
+            GpuMat d12;
+
+            GpuMat r1;
+            GpuMat r2;
+
+            GpuMat sizes;
+            int maxSize;
+
+            void create(int levels, int maxCapacity, bool isTempl);
+            void release();
+        };
+
+        typedef void (*set_func_t)(PtrStepb p1_pos, PtrStepb p1_theta, PtrStepb p2_pos, PtrStepb d12, PtrStepb r1, PtrStepb r2);
+        typedef void (*build_func_t)(const unsigned int* coordList, const float* thetaList, int pointsCount,
+                                     int* sizes, int maxSize,
+                                     float xi, float angleEpsilon, int levels,
+                                     float2 center, float maxDist);
+
+        void buildFeatureList(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, Feature& features,
+                              set_func_t set_func, build_func_t build_func, bool isTempl, Point2d center = Point2d());
+
+        void calcOrientation();
+        void calcScale(double angle);
+        void calcPosition(double angle, int angleVotes, double scale, int scaleVotes);
+
+        double xi;
+        int levels;
+        double angleEpsilon;
+
+        double minAngle;
+        double maxAngle;
+        double angleStep;
+        int angleThresh;
+
+        double minScale;
+        double maxScale;
+        double scaleStep;
+        int scaleThresh;
+
+        double dp;
+        int posThresh;
+
+        Feature templFeatures;
+        Feature imageFeatures;
+
+        vector< pair<double, int> > angles;
+        vector< pair<double, int> > scales;
+
+        GpuMat hist;
+        vector<int> h_buf;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Guil_Full, "GeneralizedHough_GPU.POSITION_SCALE_ROTATION",
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "maxSize", obj.maxSize, false, 0, 0,
+                                           "Maximal size of inner buffers.");
+                      obj.info()->addParam(obj, "xi", obj.xi, false, 0, 0,
+                                           "Angle difference in degrees between two points in feature.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "Feature table levels.");
+                      obj.info()->addParam(obj, "angleEpsilon", obj.angleEpsilon, false, 0, 0,
+                                           "Maximal difference between angles that treated as equal.");
+                      obj.info()->addParam(obj, "minAngle", obj.minAngle, false, 0, 0,
+                                           "Minimal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "maxAngle", obj.maxAngle, false, 0, 0,
+                                           "Maximal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "angleStep", obj.angleStep, false, 0, 0,
+                                           "Angle step in degrees.");
+                      obj.info()->addParam(obj, "angleThresh", obj.angleThresh, false, 0, 0,
+                                           "Angle threshold.");
+                      obj.info()->addParam(obj, "minScale", obj.minScale, false, 0, 0,
+                                           "Minimal scale to detect.");
+                      obj.info()->addParam(obj, "maxScale", obj.maxScale, false, 0, 0,
+                                           "Maximal scale to detect.");
+                      obj.info()->addParam(obj, "scaleStep", obj.scaleStep, false, 0, 0,
+                                           "Scale step.");
+                      obj.info()->addParam(obj, "scaleThresh", obj.scaleThresh, false, 0, 0,
+                                           "Scale threshold.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution.");
+                      obj.info()->addParam(obj, "posThresh", obj.posThresh, false, 0, 0,
+                                           "Position threshold."));
+
+    GHT_Guil_Full::GHT_Guil_Full()
+    {
+        maxSize = 1000;
+        xi = 90.0;
+        levels = 360;
+        angleEpsilon = 1.0;
+
+        minAngle = 0.0;
+        maxAngle = 360.0;
+        angleStep = 1.0;
+        angleThresh = 15000;
+
+        minScale = 0.5;
+        maxScale = 2.0;
+        scaleStep = 0.05;
+        scaleThresh = 1000;
+
+        dp = 1.0;
+        posThresh = 100;
+    }
+
+    void GHT_Guil_Full::releaseImpl()
+    {
+        GHT_Pos::releaseImpl();
+
+        templFeatures.release();
+        imageFeatures.release();
+
+        releaseVector(angles);
+        releaseVector(scales);
+
+        hist.release();
+        releaseVector(h_buf);
+    }
+
+    void GHT_Guil_Full::processTempl()
+    {
+        using namespace cv::gpu::device::hough;
+
+        buildFeatureList(templEdges, templDx, templDy, templFeatures,
+            GHT_Guil_Full_setTemplFeatures, GHT_Guil_Full_buildTemplFeatureList_gpu,
+            true, templCenter);
+
+        h_buf.resize(templFeatures.sizes.cols);
+        cudaSafeCall( cudaMemcpy(&h_buf[0], templFeatures.sizes.data, h_buf.size() * sizeof(int), cudaMemcpyDeviceToHost) );
+        templFeatures.maxSize = *max_element(h_buf.begin(), h_buf.end());
+    }
+
+    void GHT_Guil_Full::processImage()
+    {
+        using namespace cv::gpu::device::hough;
+
+        CV_Assert(levels > 0);
+        CV_Assert(templFeatures.sizes.cols == levels + 1);
+        CV_Assert(minAngle >= 0.0 && minAngle < maxAngle && maxAngle <= 360.0);
+        CV_Assert(angleStep > 0.0 && angleStep < 360.0);
+        CV_Assert(angleThresh > 0);
+        CV_Assert(minScale > 0.0 && minScale < maxScale);
+        CV_Assert(scaleStep > 0.0);
+        CV_Assert(scaleThresh > 0);
+        CV_Assert(dp > 0.0);
+        CV_Assert(posThresh > 0);
+
+        const double iAngleStep = 1.0 / angleStep;
+        const int angleRange = cvCeil((maxAngle - minAngle) * iAngleStep);
+
+        const double iScaleStep = 1.0 / scaleStep;
+        const int scaleRange = cvCeil((maxScale - minScale) * iScaleStep);
+
+        const double idp = 1.0 / dp;
+        const int histRows = cvCeil(imageSize.height * idp);
+        const int histCols = cvCeil(imageSize.width * idp);
+
+        ensureSizeIsEnough(histRows + 2, std::max(angleRange + 1, std::max(scaleRange + 1, histCols + 2)), CV_32SC1, hist);
+        h_buf.resize(std::max(angleRange + 1, scaleRange + 1));
+
+        ensureSizeIsEnough(2, maxSize, CV_32FC4, outBuf);
+
+        buildFeatureList(imageEdges, imageDx, imageDy, imageFeatures,
+            GHT_Guil_Full_setImageFeatures, GHT_Guil_Full_buildImageFeatureList_gpu,
+            false);
+
+        calcOrientation();
+
+        for (size_t i = 0; i < angles.size(); ++i)
+        {
+            const double angle = angles[i].first;
+            const int angleVotes = angles[i].second;
+
+            calcScale(angle);
+
+            for (size_t j = 0; j < scales.size(); ++j)
+            {
+                const double scale = scales[j].first;
+                const int scaleVotes = scales[j].second;
+
+                calcPosition(angle, angleVotes, scale, scaleVotes);
+            }
+        }
+    }
+
+    void GHT_Guil_Full::Feature::create(int levels, int maxCapacity, bool isTempl)
+    {
+        if (!isTempl)
+        {
+            ensureSizeIsEnough(levels + 1, maxCapacity, CV_32FC2, p1_pos);
+            ensureSizeIsEnough(levels + 1, maxCapacity, CV_32FC2, p2_pos);
+        }
+
+        ensureSizeIsEnough(levels + 1, maxCapacity, CV_32FC1, p1_theta);
+
+        ensureSizeIsEnough(levels + 1, maxCapacity, CV_32FC1, d12);
+
+        if (isTempl)
+        {
+            ensureSizeIsEnough(levels + 1, maxCapacity, CV_32FC2, r1);
+            ensureSizeIsEnough(levels + 1, maxCapacity, CV_32FC2, r2);
+        }
+
+        ensureSizeIsEnough(1, levels + 1, CV_32SC1, sizes);
+        sizes.setTo(Scalar::all(0));
+
+        maxSize = 0;
+    }
+
+    void GHT_Guil_Full::Feature::release()
+    {
+        p1_pos.release();
+        p1_theta.release();
+        p2_pos.release();
+
+        d12.release();
+
+        r1.release();
+        r2.release();
+
+        sizes.release();
+
+        maxSize = 0;
+    }
+
+    void GHT_Guil_Full::buildFeatureList(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, Feature& features,
+                                         set_func_t set_func, build_func_t build_func, bool isTempl, Point2d center)
+    {
+        CV_Assert(levels > 0);
+
+        const double maxDist = sqrt((double) templSize.width * templSize.width + templSize.height * templSize.height) * maxScale;
+
+        features.create(levels, maxSize, isTempl);
+        set_func(features.p1_pos, features.p1_theta, features.p2_pos, features.d12, features.r1, features.r2);
+
+        buildEdgePointList(edges, dx, dy);
+
+        if (edgePointList.cols > 0)
+        {
+            build_func(edgePointList.ptr<unsigned int>(0), edgePointList.ptr<float>(1), edgePointList.cols,
+                features.sizes.ptr<int>(), maxSize, xi, angleEpsilon, levels, make_float2(center.x, center.y), maxDist);
+        }
+    }
+
+    void GHT_Guil_Full::calcOrientation()
+    {
+        using namespace cv::gpu::device::hough;
+
+        const double iAngleStep = 1.0 / angleStep;
+        const int angleRange = cvCeil((maxAngle - minAngle) * iAngleStep);
+
+        hist.setTo(Scalar::all(0));
+        GHT_Guil_Full_calcOHist_gpu(templFeatures.sizes.ptr<int>(), imageFeatures.sizes.ptr<int>(0),
+            hist.ptr<int>(), minAngle, maxAngle, angleStep, angleRange, levels, templFeatures.maxSize);
+        cudaSafeCall( cudaMemcpy(&h_buf[0], hist.data, h_buf.size() * sizeof(int), cudaMemcpyDeviceToHost) );
+
+        angles.clear();
+
+        for (int n = 0; n < angleRange; ++n)
+        {
+            if (h_buf[n] >= angleThresh)
+            {
+                const double angle = minAngle + n * angleStep;
+                angles.push_back(make_pair(angle, h_buf[n]));
+            }
+        }
+    }
+
+    void GHT_Guil_Full::calcScale(double angle)
+    {
+        using namespace cv::gpu::device::hough;
+
+        const double iScaleStep = 1.0 / scaleStep;
+        const int scaleRange = cvCeil((maxScale - minScale) * iScaleStep);
+
+        hist.setTo(Scalar::all(0));
+        GHT_Guil_Full_calcSHist_gpu(templFeatures.sizes.ptr<int>(), imageFeatures.sizes.ptr<int>(0),
+            hist.ptr<int>(), angle, angleEpsilon, minScale, maxScale, iScaleStep, scaleRange, levels, templFeatures.maxSize);
+        cudaSafeCall( cudaMemcpy(&h_buf[0], hist.data, h_buf.size() * sizeof(int), cudaMemcpyDeviceToHost) );
+
+        scales.clear();
+
+        for (int s = 0; s < scaleRange; ++s)
+        {
+            if (h_buf[s] >= scaleThresh)
+            {
+                const double scale = minScale + s * scaleStep;
+                scales.push_back(make_pair(scale, h_buf[s]));
+            }
+        }
+    }
+
+    void GHT_Guil_Full::calcPosition(double angle, int angleVotes, double scale, int scaleVotes)
+    {
+        using namespace cv::gpu::device::hough;
+
+        hist.setTo(Scalar::all(0));
+        GHT_Guil_Full_calcPHist_gpu(templFeatures.sizes.ptr<int>(), imageFeatures.sizes.ptr<int>(0),
+            hist, angle, angleEpsilon, scale, dp, levels, templFeatures.maxSize);
+
+        posCount = GHT_Guil_Full_findPosInHist_gpu(hist, outBuf.ptr<float4>(0), outBuf.ptr<int3>(1),
+            posCount, maxSize, angle, angleVotes, scale, scaleVotes, dp, posThresh);
+    }
+}
+
+Ptr<GeneralizedHough_GPU> cv::gpu::GeneralizedHough_GPU::create(int method)
+{
+    switch (method)
+    {
+    case GHT_POSITION:
+        CV_Assert( !GHT_Ballard_Pos_info_auto.name().empty() );
+        return new GHT_Ballard_Pos();
+
+    case (GHT_POSITION | GHT_SCALE):
+        CV_Assert( !GHT_Ballard_PosScale_info_auto.name().empty() );
+        return new GHT_Ballard_PosScale();
+
+    case (GHT_POSITION | GHT_ROTATION):
+        CV_Assert( !GHT_Ballard_PosRotation_info_auto.name().empty() );
+        return new GHT_Ballard_PosRotation();
+
+    case (GHT_POSITION | GHT_SCALE | GHT_ROTATION):
+        CV_Assert( !GHT_Guil_Full_info_auto.name().empty() );
+        return new GHT_Guil_Full();
+    }
+
+    CV_Error(CV_StsBadArg, "Unsupported method");
+    return Ptr<GeneralizedHough_GPU>();
+}
+
+cv::gpu::GeneralizedHough_GPU::~GeneralizedHough_GPU()
+{
+}
+
+void cv::gpu::GeneralizedHough_GPU::setTemplate(const GpuMat& templ, int cannyThreshold, Point templCenter)
+{
+    CV_Assert(templ.type() == CV_8UC1);
+    CV_Assert(cannyThreshold > 0);
+
+    ensureSizeIsEnough(templ.size(), CV_8UC1, edges_);
+    Canny(templ, cannyBuf_, edges_, cannyThreshold / 2, cannyThreshold);
+
+    if (templCenter == Point(-1, -1))
+        templCenter = Point(templ.cols / 2, templ.rows / 2);
+
+    setTemplateImpl(edges_, cannyBuf_.dx, cannyBuf_.dy, templCenter);
+}
+
+void cv::gpu::GeneralizedHough_GPU::setTemplate(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, Point templCenter)
+{
+    if (templCenter == Point(-1, -1))
+        templCenter = Point(edges.cols / 2, edges.rows / 2);
+
+    setTemplateImpl(edges, dx, dy, templCenter);
+}
+
+void cv::gpu::GeneralizedHough_GPU::detect(const GpuMat& image, GpuMat& positions, int cannyThreshold)
+{
+    CV_Assert(image.type() == CV_8UC1);
+    CV_Assert(cannyThreshold > 0);
+
+    ensureSizeIsEnough(image.size(), CV_8UC1, edges_);
+    Canny(image, cannyBuf_, edges_, cannyThreshold / 2, cannyThreshold);
+
+    detectImpl(edges_, cannyBuf_.dx, cannyBuf_.dy, positions);
+}
+
+void cv::gpu::GeneralizedHough_GPU::detect(const GpuMat& edges, const GpuMat& dx, const GpuMat& dy, GpuMat& positions)
+{
+    detectImpl(edges, dx, dy, positions);
+}
+
+void cv::gpu::GeneralizedHough_GPU::download(const GpuMat& d_positions, OutputArray h_positions_, OutputArray h_votes_)
+{
+    if (d_positions.empty())
+    {
+        h_positions_.release();
+        if (h_votes_.needed())
+            h_votes_.release();
+        return;
+    }
+
+    CV_Assert(d_positions.rows == 2 && d_positions.type() == CV_32FC4);
+
+    h_positions_.create(1, d_positions.cols, CV_32FC4);
+    Mat h_positions = h_positions_.getMat();
+    d_positions.row(0).download(h_positions);
+
+    if (h_votes_.needed())
+    {
+        h_votes_.create(1, d_positions.cols, CV_32SC3);
+        Mat h_votes = h_votes_.getMat();
+        GpuMat d_votes(1, d_positions.cols, CV_32SC3, const_cast<int3*>(d_positions.ptr<int3>(1)));
+        d_votes.download(h_votes);
+    }
+}
+
+void cv::gpu::GeneralizedHough_GPU::release()
+{
+    edges_.release();
+    cannyBuf_.release();
+    releaseImpl();
+}
+
 #endif /* !defined (HAVE_CUDA) */

modules/gpu/test/test_hough.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.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * 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 "test_precomp.hpp"
+
+#ifdef HAVE_CUDA
+
+namespace {
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////
+// HoughLines
+
+PARAM_TEST_CASE(HoughLines, cv::gpu::DeviceInfo, cv::Size, UseRoi)
+{
+    static void generateLines(cv::Mat& img)
+    {
+        img.setTo(cv::Scalar::all(0));
+
+        cv::line(img, cv::Point(20, 0), cv::Point(20, img.rows), cv::Scalar::all(255));
+        cv::line(img, cv::Point(0, 50), cv::Point(img.cols, 50), cv::Scalar::all(255));
+        cv::line(img, cv::Point(0, 0), cv::Point(img.cols, img.rows), cv::Scalar::all(255));
+        cv::line(img, cv::Point(img.cols, 0), cv::Point(0, img.rows), cv::Scalar::all(255));
+    }
+
+    static void drawLines(cv::Mat& dst, const std::vector<cv::Vec2f>& lines)
+    {
+        dst.setTo(cv::Scalar::all(0));
+
+        for (size_t i = 0; i < lines.size(); ++i)
+        {
+            float rho = lines[i][0], theta = lines[i][1];
+            cv::Point pt1, pt2;
+            double a = std::cos(theta), b = std::sin(theta);
+            double x0 = a*rho, y0 = b*rho;
+            pt1.x = cvRound(x0 + 1000*(-b));
+            pt1.y = cvRound(y0 + 1000*(a));
+            pt2.x = cvRound(x0 - 1000*(-b));
+            pt2.y = cvRound(y0 - 1000*(a));
+            cv::line(dst, pt1, pt2, cv::Scalar::all(255));
+        }
+    }
+};
+
+TEST_P(HoughLines, Accuracy)
+{
+    const cv::gpu::DeviceInfo devInfo = GET_PARAM(0);
+    cv::gpu::setDevice(devInfo.deviceID());
+    const cv::Size size = GET_PARAM(1);
+    const bool useRoi = GET_PARAM(2);
+
+    const float rho = 1.0f;
+    const float theta = 1.5f * CV_PI / 180.0f;
+    const int threshold = 100;
+
+    cv::Mat src(size, CV_8UC1);
+    generateLines(src);
+
+    cv::gpu::GpuMat d_lines;
+    cv::gpu::HoughLines(loadMat(src, useRoi), d_lines, rho, theta, threshold);
+
+    std::vector<cv::Vec2f> lines;
+    cv::gpu::HoughLinesDownload(d_lines, lines);
+
+    cv::Mat dst(size, CV_8UC1);
+    drawLines(dst, lines);
+
+    ASSERT_MAT_NEAR(src, dst, 0.0);
+}
+
+INSTANTIATE_TEST_CASE_P(GPU_ImgProc, HoughLines, testing::Combine(
+    ALL_DEVICES,
+    DIFFERENT_SIZES,
+    WHOLE_SUBMAT));
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////
+// HoughCircles
+
+PARAM_TEST_CASE(HoughCircles, cv::gpu::DeviceInfo, cv::Size, UseRoi)
+{
+    static void drawCircles(cv::Mat& dst, const std::vector<cv::Vec3f>& circles, bool fill)
+    {
+        dst.setTo(cv::Scalar::all(0));
+
+        for (size_t i = 0; i < circles.size(); ++i)
+            cv::circle(dst, cv::Point2f(circles[i][0], circles[i][1]), (int)circles[i][2], cv::Scalar::all(255), fill ? -1 : 1);
+    }
+};
+
+TEST_P(HoughCircles, Accuracy)
+{
+    const cv::gpu::DeviceInfo devInfo = GET_PARAM(0);
+    cv::gpu::setDevice(devInfo.deviceID());
+    const cv::Size size = GET_PARAM(1);
+    const bool useRoi = GET_PARAM(2);
+
+    const float dp = 2.0f;
+    const float minDist = 10.0f;
+    const int minRadius = 10;
+    const int maxRadius = 20;
+    const int cannyThreshold = 100;
+    const int votesThreshold = 20;
+
+    std::vector<cv::Vec3f> circles_gold(4);
+    circles_gold[0] = cv::Vec3i(20, 20, minRadius);
+    circles_gold[1] = cv::Vec3i(90, 87, minRadius + 3);
+    circles_gold[2] = cv::Vec3i(30, 70, minRadius + 8);
+    circles_gold[3] = cv::Vec3i(80, 10, maxRadius);
+
+    cv::Mat src(size, CV_8UC1);
+    drawCircles(src, circles_gold, true);
+
+    cv::gpu::GpuMat d_circles;
+    cv::gpu::HoughCircles(loadMat(src, useRoi), d_circles, CV_HOUGH_GRADIENT, dp, minDist, cannyThreshold, votesThreshold, minRadius, maxRadius);
+
+    std::vector<cv::Vec3f> circles;
+    cv::gpu::HoughCirclesDownload(d_circles, circles);
+
+    ASSERT_FALSE(circles.empty());
+
+    for (size_t i = 0; i < circles.size(); ++i)
+    {
+        cv::Vec3f cur = circles[i];
+
+        bool found = false;
+
+        for (size_t j = 0; j < circles_gold.size(); ++j)
+        {
+            cv::Vec3f gold = circles_gold[j];
+
+            if (std::fabs(cur[0] - gold[0]) < minDist && std::fabs(cur[1] - gold[1]) < minDist && std::fabs(cur[2] - gold[2]) < minDist)
+            {
+                found = true;
+                break;
+            }
+        }
+
+        ASSERT_TRUE(found);
+    }
+}
+
+INSTANTIATE_TEST_CASE_P(GPU_ImgProc, HoughCircles, testing::Combine(
+    ALL_DEVICES,
+    DIFFERENT_SIZES,
+    WHOLE_SUBMAT));
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////
+// GeneralizedHough
+
+PARAM_TEST_CASE(GeneralizedHough, cv::gpu::DeviceInfo, UseRoi)
+{
+};
+
+TEST_P(GeneralizedHough, POSITION)
+{
+    const cv::gpu::DeviceInfo devInfo = GET_PARAM(0);
+    cv::gpu::setDevice(devInfo.deviceID());
+    const bool useRoi = GET_PARAM(1);
+
+    cv::Mat templ = readImage("../cv/shared/templ.png", cv::IMREAD_GRAYSCALE);
+    ASSERT_FALSE(templ.empty());
+
+    cv::Point templCenter(templ.cols / 2, templ.rows / 2);
+
+    const size_t gold_count = 3;
+    cv::Point pos_gold[gold_count];
+    pos_gold[0] = cv::Point(templCenter.x + 10, templCenter.y + 10);
+    pos_gold[1] = cv::Point(2 * templCenter.x + 40, templCenter.y + 10);
+    pos_gold[2] = cv::Point(2 * templCenter.x + 40, 2 * templCenter.y + 40);
+
+    cv::Mat image(templ.rows * 3, templ.cols * 3, CV_8UC1, cv::Scalar::all(0));
+    for (size_t i = 0; i < gold_count; ++i)
+    {
+        cv::Rect rec(pos_gold[i].x - templCenter.x, pos_gold[i].y - templCenter.y, templ.cols, templ.rows);
+        cv::Mat imageROI = image(rec);
+        templ.copyTo(imageROI);
+    }
+
+    cv::Ptr<cv::gpu::GeneralizedHough_GPU> hough = cv::gpu::GeneralizedHough_GPU::create(cv::GHT_POSITION);
+    hough->set("votesThreshold", 200);
+
+    hough->setTemplate(loadMat(templ, useRoi));
+
+    cv::gpu::GpuMat d_pos;
+    hough->detect(loadMat(image, useRoi), d_pos);
+
+    std::vector<cv::Vec4f> pos;
+    hough->download(d_pos, pos);
+
+    ASSERT_EQ(gold_count, pos.size());
+
+    for (size_t i = 0; i < gold_count; ++i)
+    {
+        cv::Point gold = pos_gold[i];
+
+        bool found = false;
+
+        for (size_t j = 0; j < pos.size(); ++j)
+        {
+            cv::Point2f p(pos[j][0], pos[j][1]);
+
+            if (::fabs(p.x - gold.x) < 2 && ::fabs(p.y - gold.y) < 2)
+            {
+                found = true;
+                break;
+            }
+        }
+
+        ASSERT_TRUE(found);
+    }
+}
+
+INSTANTIATE_TEST_CASE_P(GPU_ImgProc, GeneralizedHough, testing::Combine(
+    ALL_DEVICES,
+    WHOLE_SUBMAT));
+
+} // namespace
+
+#endif // HAVE_CUDA

modules/gpu/test/test_imgproc.cpp

@@ -1126,142 +1126,6 @@ INSTANTIATE_TEST_CASE_P(GPU_ImgProc, CornerMinEigen, testing::Combine(
     testing::Values(BlockSize(3), BlockSize(5), BlockSize(7)),
     testing::Values(ApertureSize(0), ApertureSize(3), ApertureSize(5), ApertureSize(7))));
 
-///////////////////////////////////////////////////////////////////////////////////////////////////////
-// HoughLines
-
-PARAM_TEST_CASE(HoughLines, cv::gpu::DeviceInfo, cv::Size, UseRoi)
-{
-    static void generateLines(cv::Mat& img)
-    {
-        img.setTo(cv::Scalar::all(0));
-
-        cv::line(img, cv::Point(20, 0), cv::Point(20, img.rows), cv::Scalar::all(255));
-        cv::line(img, cv::Point(0, 50), cv::Point(img.cols, 50), cv::Scalar::all(255));
-        cv::line(img, cv::Point(0, 0), cv::Point(img.cols, img.rows), cv::Scalar::all(255));
-        cv::line(img, cv::Point(img.cols, 0), cv::Point(0, img.rows), cv::Scalar::all(255));
-    }
-
-    static void drawLines(cv::Mat& dst, const std::vector<cv::Vec2f>& lines)
-    {
-        dst.setTo(cv::Scalar::all(0));
-
-        for (size_t i = 0; i < lines.size(); ++i)
-        {
-            float rho = lines[i][0], theta = lines[i][1];
-            cv::Point pt1, pt2;
-            double a = std::cos(theta), b = std::sin(theta);
-            double x0 = a*rho, y0 = b*rho;
-            pt1.x = cvRound(x0 + 1000*(-b));
-            pt1.y = cvRound(y0 + 1000*(a));
-            pt2.x = cvRound(x0 - 1000*(-b));
-            pt2.y = cvRound(y0 - 1000*(a));
-            cv::line(dst, pt1, pt2, cv::Scalar::all(255));
-        }
-    }
-};
-
-TEST_P(HoughLines, Accuracy)
-{
-    const cv::gpu::DeviceInfo devInfo = GET_PARAM(0);
-    cv::gpu::setDevice(devInfo.deviceID());
-    const cv::Size size = GET_PARAM(1);
-    const bool useRoi = GET_PARAM(2);
-
-    const float rho = 1.0f;
-    const float theta = 1.5f * CV_PI / 180.0f;
-    const int threshold = 100;
-
-    cv::Mat src(size, CV_8UC1);
-    generateLines(src);
-
-    cv::gpu::GpuMat d_lines;
-    cv::gpu::HoughLines(loadMat(src, useRoi), d_lines, rho, theta, threshold);
-
-    std::vector<cv::Vec2f> lines;
-    cv::gpu::HoughLinesDownload(d_lines, lines);
-
-    cv::Mat dst(size, CV_8UC1);
-    drawLines(dst, lines);
-
-    ASSERT_MAT_NEAR(src, dst, 0.0);
-}
-
-INSTANTIATE_TEST_CASE_P(GPU_ImgProc, HoughLines, testing::Combine(
-    ALL_DEVICES,
-    DIFFERENT_SIZES,
-    WHOLE_SUBMAT));
-
-///////////////////////////////////////////////////////////////////////////////////////////////////////
-// HoughCircles
-
-PARAM_TEST_CASE(HoughCircles, cv::gpu::DeviceInfo, cv::Size, UseRoi)
-{
-    static void drawCircles(cv::Mat& dst, const std::vector<cv::Vec3f>& circles, bool fill)
-    {
-        dst.setTo(cv::Scalar::all(0));
-
-        for (size_t i = 0; i < circles.size(); ++i)
-            cv::circle(dst, cv::Point2f(circles[i][0], circles[i][1]), (int)circles[i][2], cv::Scalar::all(255), fill ? -1 : 1);
-    }
-};
-
-TEST_P(HoughCircles, Accuracy)
-{
-    const cv::gpu::DeviceInfo devInfo = GET_PARAM(0);
-    cv::gpu::setDevice(devInfo.deviceID());
-    const cv::Size size = GET_PARAM(1);
-    const bool useRoi = GET_PARAM(2);
-
-    const float dp = 2.0f;
-    const float minDist = 10.0f;
-    const int minRadius = 10;
-    const int maxRadius = 20;
-    const int cannyThreshold = 100;
-    const int votesThreshold = 20;
-
-    std::vector<cv::Vec3f> circles_gold(4);
-    circles_gold[0] = cv::Vec3i(20, 20, minRadius);
-    circles_gold[1] = cv::Vec3i(90, 87, minRadius + 3);
-    circles_gold[2] = cv::Vec3i(30, 70, minRadius + 8);
-    circles_gold[3] = cv::Vec3i(80, 10, maxRadius);
-
-    cv::Mat src(size, CV_8UC1);
-    drawCircles(src, circles_gold, true);
-
-    cv::gpu::GpuMat d_circles;
-    cv::gpu::HoughCircles(loadMat(src, useRoi), d_circles, CV_HOUGH_GRADIENT, dp, minDist, cannyThreshold, votesThreshold, minRadius, maxRadius);
-
-    std::vector<cv::Vec3f> circles;
-    cv::gpu::HoughCirclesDownload(d_circles, circles);
-
-    ASSERT_FALSE(circles.empty());
-
-    for (size_t i = 0; i < circles.size(); ++i)
-    {
-        cv::Vec3f cur = circles[i];
-
-        bool found = false;
-
-        for (size_t j = 0; j < circles_gold.size(); ++j)
-        {
-            cv::Vec3f gold = circles_gold[j];
-
-            if (std::fabs(cur[0] - gold[0]) < minDist && std::fabs(cur[1] - gold[1]) < minDist && std::fabs(cur[2] - gold[2]) < minDist)
-            {
-                found = true;
-                break;
-            }
-        }
-
-        ASSERT_TRUE(found);
-    }
-}
-
-INSTANTIATE_TEST_CASE_P(GPU_ImgProc, HoughCircles, testing::Combine(
-    ALL_DEVICES,
-    DIFFERENT_SIZES,
-    WHOLE_SUBMAT));
-
 } // namespace
 
 #endif // HAVE_CUDA

modules/imgproc/include/opencv2/imgproc/imgproc.hpp

@@ -489,6 +489,42 @@ CV_EXPORTS_W void HoughCircles( InputArray image, OutputArray circles,
                                double param1=100, double param2=100,
                                int minRadius=0, int maxRadius=0 );
 
+enum
+{
+    GHT_POSITION = 0,
+    GHT_SCALE = 1,
+    GHT_ROTATION = 2
+};
+
+//! finds arbitrary template in the grayscale image using Generalized Hough Transform
+//! Ballard, D.H. (1981). Generalizing the Hough transform to detect arbitrary shapes. Pattern Recognition 13 (2): 111-122.
+//! Guil, N., González-Linares, J.M. and Zapata, E.L. (1999). Bidimensional shape detection using an invariant approach. Pattern Recognition 32 (6): 1025-1038.
+class CV_EXPORTS GeneralizedHough : public Algorithm
+{
+public:
+    static Ptr<GeneralizedHough> create(int method);
+
+    virtual ~GeneralizedHough();
+
+    //! set template to search
+    void setTemplate(InputArray templ, int cannyThreshold = 100, Point templCenter = Point(-1, -1));
+    void setTemplate(InputArray edges, InputArray dx, InputArray dy, Point templCenter = Point(-1, -1));
+
+    //! find template on image
+    void detect(InputArray image, OutputArray positions, OutputArray votes = cv::noArray(), int cannyThreshold = 100);
+    void detect(InputArray edges, InputArray dx, InputArray dy, OutputArray positions, OutputArray votes = cv::noArray());
+
+    void release();
+
+protected:
+    virtual void setTemplateImpl(const Mat& edges, const Mat& dx, const Mat& dy, Point templCenter) = 0;
+    virtual void detectImpl(const Mat& edges, const Mat& dx, const Mat& dy, OutputArray positions, OutputArray votes) = 0;
+    virtual void releaseImpl() = 0;
+
+private:
+    Mat edges_, dx_, dy_;
+};
+
 //! erodes the image (applies the local minimum operator)
 CV_EXPORTS_W void erode( InputArray src, OutputArray dst, InputArray kernel,
                          Point anchor=Point(-1,-1), int iterations=1,

modules/imgproc/src/generalized_hough.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.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * 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 "precomp.hpp"
+
+using namespace std;
+using namespace cv;
+
+namespace
+{
+    /////////////////////////////////////
+    // Common
+
+    template <typename T, class A> void releaseVector(vector<T, A>& v)
+    {
+        vector<T, A> empty;
+        empty.swap(v);
+    }
+
+    double toRad(double a)
+    {
+        return a * CV_PI / 180.0;
+    }
+
+    bool notNull(float v)
+    {
+        return fabs(v) > numeric_limits<float>::epsilon();
+    }
+    bool notNull(double v)
+    {
+        return fabs(v) > numeric_limits<double>::epsilon();
+    }
+
+    class GHT_Pos : public GeneralizedHough
+    {
+    public:
+        GHT_Pos();
+
+    protected:
+        void setTemplateImpl(const Mat& edges, const Mat& dx, const Mat& dy, Point templCenter);
+        void detectImpl(const Mat& edges, const Mat& dx, const Mat& dy, OutputArray positions, OutputArray votes);
+        void releaseImpl();
+
+        virtual void processTempl() = 0;
+        virtual void processImage() = 0;
+
+        void filterMinDist();
+        void convertTo(OutputArray positions, OutputArray votes);
+
+        double minDist;
+
+        Size templSize;
+        Point templCenter;
+        Mat templEdges;
+        Mat templDx;
+        Mat templDy;
+
+        Size imageSize;
+        Mat imageEdges;
+        Mat imageDx;
+        Mat imageDy;
+
+        vector<Vec4f> posOutBuf;
+        vector<Vec3i> voteOutBuf;
+    };
+
+    GHT_Pos::GHT_Pos()
+    {
+        minDist = 1.0;
+    }
+
+    void GHT_Pos::setTemplateImpl(const Mat& edges, const Mat& dx, const Mat& dy, Point templCenter_)
+    {
+        templSize = edges.size();
+        templCenter = templCenter_;
+        edges.copyTo(templEdges);
+        dx.copyTo(templDx);
+        dy.copyTo(templDy);
+
+        processTempl();
+    }
+
+    void GHT_Pos::detectImpl(const Mat& edges, const Mat& dx, const Mat& dy, OutputArray positions, OutputArray votes)
+    {
+        imageSize = edges.size();
+        edges.copyTo(imageEdges);
+        dx.copyTo(imageDx);
+        dy.copyTo(imageDy);
+
+        posOutBuf.clear();
+        voteOutBuf.clear();
+
+        processImage();
+
+        if (!posOutBuf.empty())
+        {
+            if (minDist > 1)
+                filterMinDist();
+            convertTo(positions, votes);
+        }
+        else
+        {
+            positions.release();
+            if (votes.needed())
+                votes.release();
+        }
+    }
+
+    void GHT_Pos::releaseImpl()
+    {
+        templSize = Size();
+        templCenter = Point(-1, -1);
+        templEdges.release();
+        templDx.release();
+        templDy.release();
+
+        imageSize = Size();
+        imageEdges.release();
+        imageDx.release();
+        imageDy.release();
+
+        releaseVector(posOutBuf);
+        releaseVector(voteOutBuf);
+    }
+
+    #define votes_cmp_gt(l1, l2) (aux[l1][0] > aux[l2][0])
+    static CV_IMPLEMENT_QSORT_EX( sortIndexies, size_t, votes_cmp_gt, const Vec3i* )
+
+    void GHT_Pos::filterMinDist()
+    {
+        size_t oldSize = posOutBuf.size();
+        const bool hasVotes = !voteOutBuf.empty();
+
+        CV_Assert(!hasVotes || voteOutBuf.size() == oldSize);
+
+        vector<Vec4f> oldPosBuf(posOutBuf);
+        vector<Vec3i> oldVoteBuf(voteOutBuf);
+
+        vector<size_t> indexies(oldSize);
+        for (size_t i = 0; i < oldSize; ++i)
+            indexies[i] = i;
+        sortIndexies(&indexies[0], oldSize, &oldVoteBuf[0]);
+
+        posOutBuf.clear();
+        voteOutBuf.clear();
+
+        const int cellSize = cvRound(minDist);
+        const int gridWidth = (imageSize.width + cellSize - 1) / cellSize;
+        const int gridHeight = (imageSize.height + cellSize - 1) / cellSize;
+
+        vector< vector<Point2f> > grid(gridWidth * gridHeight);
+
+        const double minDist2 = minDist * minDist;
+
+        for (size_t i = 0; i < oldSize; ++i)
+        {
+            const size_t ind = indexies[i];
+
+            Point2f p(oldPosBuf[ind][0], oldPosBuf[ind][1]);
+
+            bool good = true;
+
+            const int xCell = static_cast<int>(p.x / cellSize);
+            const int yCell = static_cast<int>(p.y / cellSize);
+
+            int x1 = xCell - 1;
+            int y1 = yCell - 1;
+            int x2 = xCell + 1;
+            int y2 = yCell + 1;
+
+            // boundary check
+            x1 = std::max(0, x1);
+            y1 = std::max(0, y1);
+            x2 = std::min(gridWidth - 1, x2);
+            y2 = std::min(gridHeight - 1, y2);
+
+            for (int yy = y1; yy <= y2; ++yy)
+            {
+                for (int xx = x1; xx <= x2; ++xx)
+                {
+                    const vector<Point2f>& m = grid[yy * gridWidth + xx];
+
+                    for(size_t j = 0; j < m.size(); ++j)
+                    {
+                        const Point2f d = p - m[j];
+
+                        if (d.ddot(d) < minDist2)
+                        {
+                            good = false;
+                            goto break_out;
+                        }
+                    }
+                }
+            }
+
+            break_out:
+
+            if(good)
+            {
+                grid[yCell * gridWidth + xCell].push_back(p);
+
+                posOutBuf.push_back(oldPosBuf[ind]);
+                if (hasVotes)
+                    voteOutBuf.push_back(oldVoteBuf[ind]);
+            }
+        }
+    }
+
+    void GHT_Pos::convertTo(OutputArray _positions, OutputArray _votes)
+    {
+        const int total = static_cast<int>(posOutBuf.size());
+        const bool hasVotes = !voteOutBuf.empty();
+
+        CV_Assert(!hasVotes || voteOutBuf.size() == posOutBuf.size());
+
+        _positions.create(1, total, CV_32FC4);
+        Mat positions = _positions.getMat();
+        Mat(1, total, CV_32FC4, &posOutBuf[0]).copyTo(positions);
+
+        if (_votes.needed())
+        {
+            if (!hasVotes)
+                _votes.release();
+            else
+            {
+                _votes.create(1, total, CV_32SC3);
+                Mat votes = _votes.getMat();
+                Mat(1, total, CV_32SC3, &voteOutBuf[0]).copyTo(votes);
+            }
+        }
+    }
+
+    /////////////////////////////////////
+    // POSITION Ballard
+
+    class GHT_Ballard_Pos : public GHT_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Ballard_Pos();
+
+    protected:
+        void releaseImpl();
+
+        void processTempl();
+        void processImage();
+
+        virtual void calcHist();
+        virtual void findPosInHist();
+
+        int levels;
+        int votesThreshold;
+        double dp;
+
+        vector< vector<Point> > r_table;
+        Mat hist;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Ballard_Pos, "GeneralizedHough.POSITION",
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "R-Table levels.");
+                      obj.info()->addParam(obj, "votesThreshold", obj.votesThreshold, false, 0, 0,
+                                           "The accumulator threshold for the template centers at the detection stage. The smaller it is, the more false positions may be detected.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution."));
+
+    GHT_Ballard_Pos::GHT_Ballard_Pos()
+    {
+        levels = 360;
+        votesThreshold = 100;
+        dp = 1.0;
+    }
+
+    void GHT_Ballard_Pos::releaseImpl()
+    {
+        GHT_Pos::releaseImpl();
+
+        releaseVector(r_table);
+        hist.release();
+    }
+
+    void GHT_Ballard_Pos::processTempl()
+    {
+        CV_Assert(templEdges.type() == CV_8UC1);
+        CV_Assert(templDx.type() == CV_32FC1 && templDx.size() == templSize);
+        CV_Assert(templDy.type() == templDx.type() && templDy.size() == templSize);
+        CV_Assert(levels > 0);
+
+        const double thetaScale = levels / 360.0;
+
+        r_table.resize(levels + 1);
+        for_each(r_table.begin(), r_table.end(), mem_fun_ref(&vector<Point>::clear));
+
+        for (int y = 0; y < templSize.height; ++y)
+        {
+            const uchar* edgesRow = templEdges.ptr(y);
+            const float* dxRow = templDx.ptr<float>(y);
+            const float* dyRow = templDy.ptr<float>(y);
+
+            for (int x = 0; x < templSize.width; ++x)
+            {
+                const Point p(x, y);
+
+                if (edgesRow[x] && (notNull(dyRow[x]) || notNull(dxRow[x])))
+                {
+                    const float theta = fastAtan2(dyRow[x], dxRow[x]);
+                    const int n = cvRound(theta * thetaScale);
+                    r_table[n].push_back(p - templCenter);
+                }
+            }
+        }
+    }
+
+    void GHT_Ballard_Pos::processImage()
+    {
+        calcHist();
+        findPosInHist();
+    }
+
+    void GHT_Ballard_Pos::calcHist()
+    {
+        CV_Assert(imageEdges.type() == CV_8UC1);
+        CV_Assert(imageDx.type() == CV_32FC1 && imageDx.size() == imageSize);
+        CV_Assert(imageDy.type() == imageDx.type() && imageDy.size() == imageSize);
+        CV_Assert(levels > 0 && r_table.size() == static_cast<size_t>(levels + 1));
+        CV_Assert(dp > 0.0);
+
+        const double thetaScale = levels / 360.0;
+        const double idp = 1.0 / dp;
+
+        hist.create(cvCeil(imageSize.height * idp) + 2, cvCeil(imageSize.width * idp) + 2, CV_32SC1);
+        hist.setTo(0);
+
+        const int rows = hist.rows - 2;
+        const int cols = hist.cols - 2;
+
+        for (int y = 0; y < imageSize.height; ++y)
+        {
+            const uchar* edgesRow = imageEdges.ptr(y);
+            const float* dxRow = imageDx.ptr<float>(y);
+            const float* dyRow = imageDy.ptr<float>(y);
+
+            for (int x = 0; x < imageSize.width; ++x)
+            {
+                const Point p(x, y);
+
+                if (edgesRow[x] && (notNull(dyRow[x]) || notNull(dxRow[x])))
+                {
+                    const float theta = fastAtan2(dyRow[x], dxRow[x]);
+                    const int n = cvRound(theta * thetaScale);
+
+                    const vector<Point>& r_row = r_table[n];
+
+                    for (size_t j = 0; j < r_row.size(); ++j)
+                    {
+                        Point c = p - r_row[j];
+
+                        c.x = cvRound(c.x * idp);
+                        c.y = cvRound(c.y * idp);
+
+                        if (c.x >= 0 && c.x < cols && c.y >= 0 && c.y < rows)
+                            ++hist.at<int>(c.y + 1, c.x + 1);
+                    }
+                }
+            }
+        }
+    }
+
+    void GHT_Ballard_Pos::findPosInHist()
+    {
+        CV_Assert(votesThreshold > 0);
+
+        const int histRows = hist.rows - 2;
+        const int histCols = hist.cols - 2;
+
+        for(int y = 0; y < histRows; ++y)
+        {
+            const int* prevRow = hist.ptr<int>(y);
+            const int* curRow = hist.ptr<int>(y + 1);
+            const int* nextRow = hist.ptr<int>(y + 2);
+
+            for(int x = 0; x < histCols; ++x)
+            {
+                const int votes = curRow[x + 1];
+
+                if (votes > votesThreshold && votes > curRow[x] && votes >= curRow[x + 2] && votes > prevRow[x + 1] && votes >= nextRow[x + 1])
+                {
+                    posOutBuf.push_back(Vec4f(static_cast<float>(x * dp), static_cast<float>(y * dp), 1.0f, 0.0f));
+                    voteOutBuf.push_back(Vec3i(votes, 0, 0));
+                }
+            }
+        }
+    }
+
+    /////////////////////////////////////
+    // POSITION & SCALE
+
+    class GHT_Ballard_PosScale : public GHT_Ballard_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Ballard_PosScale();
+
+    protected:
+        void calcHist();
+        void findPosInHist();
+
+        double minScale;
+        double maxScale;
+        double scaleStep;
+
+        class Worker;
+        friend class Worker;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Ballard_PosScale, "GeneralizedHough.POSITION_SCALE",
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "R-Table levels.");
+                      obj.info()->addParam(obj, "votesThreshold", obj.votesThreshold, false, 0, 0,
+                                           "The accumulator threshold for the template centers at the detection stage. The smaller it is, the more false positions may be detected.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution.");
+                      obj.info()->addParam(obj, "minScale", obj.minScale, false, 0, 0,
+                                           "Minimal scale to detect.");
+                      obj.info()->addParam(obj, "maxScale", obj.maxScale, false, 0, 0,
+                                           "Maximal scale to detect.");
+                      obj.info()->addParam(obj, "scaleStep", obj.scaleStep, false, 0, 0,
+                                           "Scale step."));
+
+    GHT_Ballard_PosScale::GHT_Ballard_PosScale()
+    {
+        minScale = 0.5;
+        maxScale = 2.0;
+        scaleStep = 0.05;
+    }
+
+    class GHT_Ballard_PosScale::Worker : public ParallelLoopBody
+    {
+    public:
+        explicit Worker(GHT_Ballard_PosScale* base_) : base(base_) {}
+
+        void operator ()(const Range& range) const;
+
+    private:
+        GHT_Ballard_PosScale* base;
+    };
+
+    void GHT_Ballard_PosScale::Worker::operator ()(const Range& range) const
+    {
+        const double thetaScale = base->levels / 360.0;
+        const double idp = 1.0 / base->dp;
+
+        for (int s = range.start; s < range.end; ++s)
+        {
+            const double scale = base->minScale + s * base->scaleStep;
+
+            Mat curHist(base->hist.size[1], base->hist.size[2], CV_32SC1, base->hist.ptr(s + 1), base->hist.step[1]);
+
+            for (int y = 0; y < base->imageSize.height; ++y)
+            {
+                const uchar* edgesRow = base->imageEdges.ptr(y);
+                const float* dxRow = base->imageDx.ptr<float>(y);
+                const float* dyRow = base->imageDy.ptr<float>(y);
+
+                for (int x = 0; x < base->imageSize.width; ++x)
+                {
+                    const Point2d p(x, y);
+
+                    if (edgesRow[x] && (notNull(dyRow[x]) || notNull(dxRow[x])))
+                    {
+                        const float theta = fastAtan2(dyRow[x], dxRow[x]);
+                        const int n = cvRound(theta * thetaScale);
+
+                        const vector<Point>& r_row = base->r_table[n];
+
+                        for (size_t j = 0; j < r_row.size(); ++j)
+                        {
+                            Point2d d = r_row[j];
+                            Point2d c = p - d * scale;
+
+                            c.x *= idp;
+                            c.y *= idp;
+
+                            if (c.x >= 0 && c.x < base->hist.size[2] - 2 && c.y >= 0 && c.y < base->hist.size[1] - 2)
+                                ++curHist.at<int>(cvRound(c.y + 1), cvRound(c.x + 1));
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    void GHT_Ballard_PosScale::calcHist()
+    {
+        CV_Assert(imageEdges.type() == CV_8UC1);
+        CV_Assert(imageDx.type() == CV_32FC1 && imageDx.size() == imageSize);
+        CV_Assert(imageDy.type() == imageDx.type() && imageDy.size() == imageSize);
+        CV_Assert(levels > 0 && r_table.size() == static_cast<size_t>(levels + 1));
+        CV_Assert(dp > 0.0);
+        CV_Assert(minScale > 0.0 && minScale < maxScale);
+        CV_Assert(scaleStep > 0.0);
+
+        const double idp = 1.0 / dp;
+        const int scaleRange = cvCeil((maxScale - minScale) / scaleStep);
+
+        const int sizes[] = {scaleRange + 2, cvCeil(imageSize.height * idp) + 2, cvCeil(imageSize.width * idp) + 2};
+        hist.create(3, sizes, CV_32SC1);
+        hist.setTo(0);
+
+        parallel_for_(Range(0, scaleRange), Worker(this));
+    }
+
+    void GHT_Ballard_PosScale::findPosInHist()
+    {
+        CV_Assert(votesThreshold > 0);
+
+        const int scaleRange = hist.size[0] - 2;
+        const int histRows = hist.size[1] - 2;
+        const int histCols = hist.size[2] - 2;
+
+        for (int s = 0; s < scaleRange; ++s)
+        {
+            const float scale = static_cast<float>(minScale + s * scaleStep);
+
+            const Mat prevHist(histRows + 2, histCols + 2, CV_32SC1, hist.ptr(s), hist.step[1]);
+            const Mat curHist(histRows + 2, histCols + 2, CV_32SC1, hist.ptr(s + 1), hist.step[1]);
+            const Mat nextHist(histRows + 2, histCols + 2, CV_32SC1, hist.ptr(s + 2), hist.step[1]);
+
+            for(int y = 0; y < histRows; ++y)
+            {
+                const int* prevHistRow = prevHist.ptr<int>(y + 1);
+                const int* prevRow = curHist.ptr<int>(y);
+                const int* curRow = curHist.ptr<int>(y + 1);
+                const int* nextRow = curHist.ptr<int>(y + 2);
+                const int* nextHistRow = nextHist.ptr<int>(y + 1);
+
+                for(int x = 0; x < histCols; ++x)
+                {
+                    const int votes = curRow[x + 1];
+
+                    if (votes > votesThreshold &&
+                        votes > curRow[x] &&
+                        votes >= curRow[x + 2] &&
+                        votes > prevRow[x + 1] &&
+                        votes >= nextRow[x + 1] &&
+                        votes > prevHistRow[x + 1] &&
+                        votes >= nextHistRow[x + 1])
+                    {
+                        posOutBuf.push_back(Vec4f(static_cast<float>(x * dp), static_cast<float>(y * dp), scale, 0.0f));
+                        voteOutBuf.push_back(Vec3i(votes, votes, 0));
+                    }
+                }
+            }
+        }
+    }
+
+    /////////////////////////////////////
+    // POSITION & ROTATION
+
+    class GHT_Ballard_PosRotation : public GHT_Ballard_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Ballard_PosRotation();
+
+    protected:
+        void calcHist();
+        void findPosInHist();
+
+        double minAngle;
+        double maxAngle;
+        double angleStep;
+
+        class Worker;
+        friend class Worker;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Ballard_PosRotation, "GeneralizedHough.POSITION_ROTATION",
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "R-Table levels.");
+                      obj.info()->addParam(obj, "votesThreshold", obj.votesThreshold, false, 0, 0,
+                                           "The accumulator threshold for the template centers at the detection stage. The smaller it is, the more false positions may be detected.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution.");
+                      obj.info()->addParam(obj, "minAngle", obj.minAngle, false, 0, 0,
+                                           "Minimal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "maxAngle", obj.maxAngle, false, 0, 0,
+                                           "Maximal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "angleStep", obj.angleStep, false, 0, 0,
+                                           "Angle step in degrees."));
+
+    GHT_Ballard_PosRotation::GHT_Ballard_PosRotation()
+    {
+        minAngle = 0.0;
+        maxAngle = 360.0;
+        angleStep = 1.0;
+    }
+
+    class GHT_Ballard_PosRotation::Worker : public ParallelLoopBody
+    {
+    public:
+        explicit Worker(GHT_Ballard_PosRotation* base_) : base(base_) {}
+
+        void operator ()(const Range& range) const;
+
+    private:
+        GHT_Ballard_PosRotation* base;
+    };
+
+    void GHT_Ballard_PosRotation::Worker::operator ()(const Range& range) const
+    {
+        const double thetaScale = base->levels / 360.0;
+        const double idp = 1.0 / base->dp;
+
+        for (int a = range.start; a < range.end; ++a)
+        {
+            const double angle = base->minAngle + a * base->angleStep;
+
+            const double sinA = ::sin(toRad(angle));
+            const double cosA = ::cos(toRad(angle));
+
+            Mat curHist(base->hist.size[1], base->hist.size[2], CV_32SC1, base->hist.ptr(a + 1), base->hist.step[1]);
+
+            for (int y = 0; y < base->imageSize.height; ++y)
+            {
+                const uchar* edgesRow = base->imageEdges.ptr(y);
+                const float* dxRow = base->imageDx.ptr<float>(y);
+                const float* dyRow = base->imageDy.ptr<float>(y);
+
+                for (int x = 0; x < base->imageSize.width; ++x)
+                {
+                    const Point2d p(x, y);
+
+                    if (edgesRow[x] && (notNull(dyRow[x]) || notNull(dxRow[x])))
+                    {
+                        double theta = fastAtan2(dyRow[x], dxRow[x]) - angle;
+                        if (theta < 0)
+                            theta += 360.0;
+                        const int n = cvRound(theta * thetaScale);
+
+                        const vector<Point>& r_row = base->r_table[n];
+
+                        for (size_t j = 0; j < r_row.size(); ++j)
+                        {
+                            Point2d d = r_row[j];
+                            Point2d c = p - Point2d(d.x * cosA - d.y * sinA, d.x * sinA + d.y * cosA);
+
+                            c.x *= idp;
+                            c.y *= idp;
+
+                            if (c.x >= 0 && c.x < base->hist.size[2] - 2 && c.y >= 0 && c.y < base->hist.size[1] - 2)
+                                ++curHist.at<int>(cvRound(c.y + 1), cvRound(c.x + 1));
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    void GHT_Ballard_PosRotation::calcHist()
+    {
+        CV_Assert(imageEdges.type() == CV_8UC1);
+        CV_Assert(imageDx.type() == CV_32FC1 && imageDx.size() == imageSize);
+        CV_Assert(imageDy.type() == imageDx.type() && imageDy.size() == imageSize);
+        CV_Assert(levels > 0 && r_table.size() == static_cast<size_t>(levels + 1));
+        CV_Assert(dp > 0.0);
+        CV_Assert(minAngle >= 0.0 && minAngle < maxAngle && maxAngle <= 360.0);
+        CV_Assert(angleStep > 0.0 && angleStep < 360.0);
+
+        const double idp = 1.0 / dp;
+        const int angleRange = cvCeil((maxAngle - minAngle) / angleStep);
+
+        const int sizes[] = {angleRange + 2, cvCeil(imageSize.height * idp) + 2, cvCeil(imageSize.width * idp) + 2};
+        hist.create(3, sizes, CV_32SC1);
+        hist.setTo(0);
+
+        parallel_for_(Range(0, angleRange), Worker(this));
+    }
+
+    void GHT_Ballard_PosRotation::findPosInHist()
+    {
+        CV_Assert(votesThreshold > 0);
+
+        const int angleRange = hist.size[0] - 2;
+        const int histRows = hist.size[1] - 2;
+        const int histCols = hist.size[2] - 2;
+
+        for (int a = 0; a < angleRange; ++a)
+        {
+            const float angle = static_cast<float>(minAngle + a * angleStep);
+
+            const Mat prevHist(histRows + 2, histCols + 2, CV_32SC1, hist.ptr(a), hist.step[1]);
+            const Mat curHist(histRows + 2, histCols + 2, CV_32SC1, hist.ptr(a + 1), hist.step[1]);
+            const Mat nextHist(histRows + 2, histCols + 2, CV_32SC1, hist.ptr(a + 2), hist.step[1]);
+
+            for(int y = 0; y < histRows; ++y)
+            {
+                const int* prevHistRow = prevHist.ptr<int>(y + 1);
+                const int* prevRow = curHist.ptr<int>(y);
+                const int* curRow = curHist.ptr<int>(y + 1);
+                const int* nextRow = curHist.ptr<int>(y + 2);
+                const int* nextHistRow = nextHist.ptr<int>(y + 1);
+
+                for(int x = 0; x < histCols; ++x)
+                {
+                    const int votes = curRow[x + 1];
+
+                    if (votes > votesThreshold &&
+                        votes > curRow[x] &&
+                        votes >= curRow[x + 2] &&
+                        votes > prevRow[x + 1] &&
+                        votes >= nextRow[x + 1] &&
+                        votes > prevHistRow[x + 1] &&
+                        votes >= nextHistRow[x + 1])
+                    {
+                        posOutBuf.push_back(Vec4f(static_cast<float>(x * dp), static_cast<float>(y * dp), 1.0f, angle));
+                        voteOutBuf.push_back(Vec3i(votes, 0, votes));
+                    }
+                }
+            }
+        }
+    }
+
+    /////////////////////////////////////////
+    // POSITION & SCALE & ROTATION
+
+    double clampAngle(double a)
+    {
+        double res = a;
+
+        while (res > 360.0)
+            res -= 360.0;
+        while (res < 0)
+            res += 360.0;
+
+        return res;
+    }
+
+    bool angleEq(double a, double b, double eps = 1.0)
+    {
+        return (fabs(clampAngle(a - b)) <= eps);
+    }
+
+    class GHT_Guil_Full : public GHT_Pos
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        GHT_Guil_Full();
+
+    protected:
+        void releaseImpl();
+
+        void processTempl();
+        void processImage();
+
+        struct ContourPoint
+        {
+            Point2d pos;
+            double theta;
+        };
+
+        struct Feature
+        {
+            ContourPoint p1;
+            ContourPoint p2;
+
+            double alpha12;
+            double d12;
+
+            Point2d r1;
+            Point2d r2;
+        };
+
+        void buildFeatureList(const Mat& edges, const Mat& dx, const Mat& dy, vector< vector<Feature> >& features, Point2d center = Point2d());
+        void getContourPoints(const Mat& edges, const Mat& dx, const Mat& dy, vector<ContourPoint>& points);
+
+        void calcOrientation();
+        void calcScale(double angle);
+        void calcPosition(double angle, int angleVotes, double scale, int scaleVotes);
+
+        int maxSize;
+        double xi;
+        int levels;
+        double angleEpsilon;
+
+        double minAngle;
+        double maxAngle;
+        double angleStep;
+        int angleThresh;
+
+        double minScale;
+        double maxScale;
+        double scaleStep;
+        int scaleThresh;
+
+        double dp;
+        int posThresh;
+
+        vector< vector<Feature> > templFeatures;
+        vector< vector<Feature> > imageFeatures;
+
+        vector< pair<double, int> > angles;
+        vector< pair<double, int> > scales;
+    };
+
+    CV_INIT_ALGORITHM(GHT_Guil_Full, "GeneralizedHough.POSITION_SCALE_ROTATION",
+                      obj.info()->addParam(obj, "minDist", obj.minDist, false, 0, 0,
+                                           "Minimum distance between the centers of the detected objects.");
+                      obj.info()->addParam(obj, "maxSize", obj.maxSize, false, 0, 0,
+                                           "Maximal size of inner buffers.");
+                      obj.info()->addParam(obj, "xi", obj.xi, false, 0, 0,
+                                           "Angle difference in degrees between two points in feature.");
+                      obj.info()->addParam(obj, "levels", obj.levels, false, 0, 0,
+                                           "Feature table levels.");
+                      obj.info()->addParam(obj, "angleEpsilon", obj.angleEpsilon, false, 0, 0,
+                                           "Maximal difference between angles that treated as equal.");
+                      obj.info()->addParam(obj, "minAngle", obj.minAngle, false, 0, 0,
+                                           "Minimal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "maxAngle", obj.maxAngle, false, 0, 0,
+                                           "Maximal rotation angle to detect in degrees.");
+                      obj.info()->addParam(obj, "angleStep", obj.angleStep, false, 0, 0,
+                                           "Angle step in degrees.");
+                      obj.info()->addParam(obj, "angleThresh", obj.angleThresh, false, 0, 0,
+                                           "Angle threshold.");
+                      obj.info()->addParam(obj, "minScale", obj.minScale, false, 0, 0,
+                                           "Minimal scale to detect.");
+                      obj.info()->addParam(obj, "maxScale", obj.maxScale, false, 0, 0,
+                                           "Maximal scale to detect.");
+                      obj.info()->addParam(obj, "scaleStep", obj.scaleStep, false, 0, 0,
+                                           "Scale step.");
+                      obj.info()->addParam(obj, "scaleThresh", obj.scaleThresh, false, 0, 0,
+                                           "Scale threshold.");
+                      obj.info()->addParam(obj, "dp", obj.dp, false, 0, 0,
+                                           "Inverse ratio of the accumulator resolution to the image resolution.");
+                      obj.info()->addParam(obj, "posThresh", obj.posThresh, false, 0, 0,
+                                           "Position threshold."));
+
+    GHT_Guil_Full::GHT_Guil_Full()
+    {
+        maxSize = 1000;
+        xi = 90.0;
+        levels = 360;
+        angleEpsilon = 1.0;
+
+        minAngle = 0.0;
+        maxAngle = 360.0;
+        angleStep = 1.0;
+        angleThresh = 15000;
+
+        minScale = 0.5;
+        maxScale = 2.0;
+        scaleStep = 0.05;
+        scaleThresh = 1000;
+
+        dp = 1.0;
+        posThresh = 100;
+    }
+
+    void GHT_Guil_Full::releaseImpl()
+    {
+        GHT_Pos::releaseImpl();
+
+        releaseVector(templFeatures);
+        releaseVector(imageFeatures);
+
+        releaseVector(angles);
+        releaseVector(scales);
+    }
+
+    void GHT_Guil_Full::processTempl()
+    {
+        buildFeatureList(templEdges, templDx, templDy, templFeatures, templCenter);
+    }
+
+    void GHT_Guil_Full::processImage()
+    {
+        buildFeatureList(imageEdges, imageDx, imageDy, imageFeatures);
+
+        calcOrientation();
+
+        for (size_t i = 0; i < angles.size(); ++i)
+        {
+            const double angle = angles[i].first;
+            const int angleVotes = angles[i].second;
+
+            calcScale(angle);
+
+            for (size_t j = 0; j < scales.size(); ++j)
+            {
+                const double scale = scales[j].first;
+                const int scaleVotes = scales[j].second;
+
+                calcPosition(angle, angleVotes, scale, scaleVotes);
+            }
+        }
+    }
+
+    void GHT_Guil_Full::buildFeatureList(const Mat& edges, const Mat& dx, const Mat& dy, vector< vector<Feature> >& features, Point2d center)
+    {
+        CV_Assert(levels > 0);
+
+        const double maxDist = sqrt((double) templSize.width * templSize.width + templSize.height * templSize.height) * maxScale;
+
+        const double alphaScale = levels / 360.0;
+
+        vector<ContourPoint> points;
+        getContourPoints(edges, dx, dy, points);
+
+        features.resize(levels + 1);
+        for_each(features.begin(), features.end(), mem_fun_ref(&vector<Feature>::clear));
+        for_each(features.begin(), features.end(), bind2nd(mem_fun_ref(&vector<Feature>::reserve), maxSize));
+
+        for (size_t i = 0; i < points.size(); ++i)
+        {
+            ContourPoint p1 = points[i];
+
+            for (size_t j = 0; j < points.size(); ++j)
+            {
+                ContourPoint p2 = points[j];
+
+                if (angleEq(p1.theta - p2.theta, xi, angleEpsilon))
+                {
+                    const Point2d d = p1.pos - p2.pos;
+
+                    Feature f;
+
+                    f.p1 = p1;
+                    f.p2 = p2;
+
+                    f.alpha12 = clampAngle(fastAtan2(d.y, d.x) - p1.theta);
+                    f.d12 = norm(d);
+
+                    if (f.d12 > maxDist)
+                        continue;
+
+                    f.r1 = p1.pos - center;
+                    f.r2 = p2.pos - center;
+
+                    const int n = cvRound(f.alpha12 * alphaScale);
+
+                    if (features[n].size() < static_cast<size_t>(maxSize))
+                        features[n].push_back(f);
+                }
+            }
+        }
+    }
+
+    void GHT_Guil_Full::getContourPoints(const Mat& edges, const Mat& dx, const Mat& dy, vector<ContourPoint>& points)
+    {
+        CV_Assert(edges.type() == CV_8UC1);
+        CV_Assert(dx.type() == CV_32FC1 && dx.size == edges.size);
+        CV_Assert(dy.type() == dx.type() && dy.size == edges.size);
+
+        points.clear();
+        points.reserve(edges.size().area());
+
+        for (int y = 0; y < edges.rows; ++y)
+        {
+            const uchar* edgesRow = edges.ptr(y);
+            const float* dxRow = dx.ptr<float>(y);
+            const float* dyRow = dy.ptr<float>(y);
+
+            for (int x = 0; x < edges.cols; ++x)
+            {
+                if (edgesRow[x] && (notNull(dyRow[x]) || notNull(dxRow[x])))
+                {
+                    ContourPoint p;
+
+                    p.pos = Point2d(x, y);
+                    p.theta = fastAtan2(dyRow[x], dxRow[x]);
+
+                    points.push_back(p);
+                }
+            }
+        }
+    }
+
+    void GHT_Guil_Full::calcOrientation()
+    {
+        CV_Assert(levels > 0);
+        CV_Assert(templFeatures.size() == static_cast<size_t>(levels + 1));
+        CV_Assert(imageFeatures.size() == templFeatures.size());
+        CV_Assert(minAngle >= 0.0 && minAngle < maxAngle && maxAngle <= 360.0);
+        CV_Assert(angleStep > 0.0 && angleStep < 360.0);
+        CV_Assert(angleThresh > 0);
+
+        const double iAngleStep = 1.0 / angleStep;
+        const int angleRange = cvCeil((maxAngle - minAngle) * iAngleStep);
+
+        vector<int> OHist(angleRange + 1, 0);
+        for (int i = 0; i <= levels; ++i)
+        {
+            const vector<Feature>& templRow = templFeatures[i];
+            const vector<Feature>& imageRow = imageFeatures[i];
+
+            for (size_t j = 0; j < templRow.size(); ++j)
+            {
+                Feature templF = templRow[j];
+
+                for (size_t k = 0; k < imageRow.size(); ++k)
+                {
+                    Feature imF = imageRow[k];
+
+                    const double angle = clampAngle(imF.p1.theta - templF.p1.theta);
+                    if (angle >= minAngle && angle <= maxAngle)
+                    {
+                        const int n = cvRound((angle - minAngle) * iAngleStep);
+                        ++OHist[n];
+                    }
+                }
+            }
+        }
+
+        angles.clear();
+
+        for (int n = 0; n < angleRange; ++n)
+        {
+            if (OHist[n] >= angleThresh)
+            {
+                const double angle = minAngle + n * angleStep;
+                angles.push_back(make_pair(angle, OHist[n]));
+            }
+        }
+    }
+
+    void GHT_Guil_Full::calcScale(double angle)
+    {
+        CV_Assert(levels > 0);
+        CV_Assert(templFeatures.size() == static_cast<size_t>(levels + 1));
+        CV_Assert(imageFeatures.size() == templFeatures.size());
+        CV_Assert(minScale > 0.0 && minScale < maxScale);
+        CV_Assert(scaleStep > 0.0);
+        CV_Assert(scaleThresh > 0);
+
+        const double iScaleStep = 1.0 / scaleStep;
+        const int scaleRange = cvCeil((maxScale - minScale) * iScaleStep);
+
+        vector<int> SHist(scaleRange + 1, 0);
+
+        for (int i = 0; i <= levels; ++i)
+        {
+            const vector<Feature>& templRow = templFeatures[i];
+            const vector<Feature>& imageRow = imageFeatures[i];
+
+            for (size_t j = 0; j < templRow.size(); ++j)
+            {
+                Feature templF = templRow[j];
+
+                templF.p1.theta += angle;
+
+                for (size_t k = 0; k < imageRow.size(); ++k)
+                {
+                    Feature imF = imageRow[k];
+
+                    if (angleEq(imF.p1.theta, templF.p1.theta, angleEpsilon))
+                    {
+                        const double scale = imF.d12 / templF.d12;
+                        if (scale >= minScale && scale <= maxScale)
+                        {
+                            const int s = cvRound((scale - minScale) * iScaleStep);
+                            ++SHist[s];
+                        }
+                    }
+                }
+            }
+        }
+
+        scales.clear();
+
+        for (int s = 0; s < scaleRange; ++s)
+        {
+            if (SHist[s] >= scaleThresh)
+            {
+                const double scale = minScale + s * scaleStep;
+                scales.push_back(make_pair(scale, SHist[s]));
+            }
+        }
+    }
+
+    void GHT_Guil_Full::calcPosition(double angle, int angleVotes, double scale, int scaleVotes)
+    {
+        CV_Assert(levels > 0);
+        CV_Assert(templFeatures.size() == static_cast<size_t>(levels + 1));
+        CV_Assert(imageFeatures.size() == templFeatures.size());
+        CV_Assert(dp > 0.0);
+        CV_Assert(posThresh > 0);
+
+        const double sinVal = sin(toRad(angle));
+        const double cosVal = cos(toRad(angle));
+        const double idp = 1.0 / dp;
+
+        const int histRows = cvCeil(imageSize.height * idp);
+        const int histCols = cvCeil(imageSize.width * idp);
+
+        Mat DHist(histRows + 2, histCols + 2, CV_32SC1, Scalar::all(0));
+
+        for (int i = 0; i <= levels; ++i)
+        {
+            const vector<Feature>& templRow = templFeatures[i];
+            const vector<Feature>& imageRow = imageFeatures[i];
+
+            for (size_t j = 0; j < templRow.size(); ++j)
+            {
+                Feature templF = templRow[j];
+
+                templF.p1.theta += angle;
+
+                templF.r1 *= scale;
+                templF.r2 *= scale;
+
+                templF.r1 = Point2d(cosVal * templF.r1.x - sinVal * templF.r1.y, sinVal * templF.r1.x + cosVal * templF.r1.y);
+                templF.r2 = Point2d(cosVal * templF.r2.x - sinVal * templF.r2.y, sinVal * templF.r2.x + cosVal * templF.r2.y);
+
+                for (size_t k = 0; k < imageRow.size(); ++k)
+                {
+                    Feature imF = imageRow[k];
+
+                    if (angleEq(imF.p1.theta, templF.p1.theta, angleEpsilon))
+                    {
+                        Point2d c1, c2;
+
+                        c1 = imF.p1.pos - templF.r1;
+                        c1 *= idp;
+
+                        c2 = imF.p2.pos - templF.r2;
+                        c2 *= idp;
+
+                        if (fabs(c1.x - c2.x) > 1 || fabs(c1.y - c2.y) > 1)
+                            continue;
+
+                        if (c1.y >= 0 && c1.y < histRows && c1.x >= 0 && c1.x < histCols)
+                            ++DHist.at<int>(cvRound(c1.y) + 1, cvRound(c1.x) + 1);
+                    }
+                }
+            }
+        }
+
+        for(int y = 0; y < histRows; ++y)
+        {
+            const int* prevRow = DHist.ptr<int>(y);
+            const int* curRow = DHist.ptr<int>(y + 1);
+            const int* nextRow = DHist.ptr<int>(y + 2);
+
+            for(int x = 0; x < histCols; ++x)
+            {
+                const int votes = curRow[x + 1];
+
+                if (votes > posThresh && votes > curRow[x] && votes >= curRow[x + 2] && votes > prevRow[x + 1] && votes >= nextRow[x + 1])
+                {
+                    posOutBuf.push_back(Vec4f(static_cast<float>(x * dp), static_cast<float>(y * dp), static_cast<float>(scale), static_cast<float>(angle)));
+                    voteOutBuf.push_back(Vec3i(votes, scaleVotes, angleVotes));
+                }
+            }
+        }
+    }
+}
+
+Ptr<GeneralizedHough> cv::GeneralizedHough::create(int method)
+{
+    switch (method)
+    {
+    case GHT_POSITION:
+        CV_Assert( !GHT_Ballard_Pos_info_auto.name().empty() );
+        return new GHT_Ballard_Pos();
+
+    case (GHT_POSITION | GHT_SCALE):
+        CV_Assert( !GHT_Ballard_PosScale_info_auto.name().empty() );
+        return new GHT_Ballard_PosScale();
+
+    case (GHT_POSITION | GHT_ROTATION):
+        CV_Assert( !GHT_Ballard_PosRotation_info_auto.name().empty() );
+        return new GHT_Ballard_PosRotation();
+
+    case (GHT_POSITION | GHT_SCALE | GHT_ROTATION):
+        CV_Assert( !GHT_Guil_Full_info_auto.name().empty() );
+        return new GHT_Guil_Full();
+    }
+
+    CV_Error(CV_StsBadArg, "Unsupported method");
+    return Ptr<GeneralizedHough>();
+}
+
+cv::GeneralizedHough::~GeneralizedHough()
+{
+}
+
+void cv::GeneralizedHough::setTemplate(InputArray _templ, int cannyThreshold, Point templCenter)
+{
+    Mat templ = _templ.getMat();
+
+    CV_Assert(templ.type() == CV_8UC1);
+    CV_Assert(cannyThreshold > 0);
+
+    Canny(templ, edges_, cannyThreshold / 2, cannyThreshold);
+    Sobel(templ, dx_, CV_32F, 1, 0);
+    Sobel(templ, dy_, CV_32F, 0, 1);
+
+    if (templCenter == Point(-1, -1))
+        templCenter = Point(templ.cols / 2, templ.rows / 2);
+
+    setTemplateImpl(edges_, dx_, dy_, templCenter);
+}
+
+void cv::GeneralizedHough::setTemplate(InputArray _edges, InputArray _dx, InputArray _dy, Point templCenter)
+{
+    Mat edges = _edges.getMat();
+    Mat dx = _dx.getMat();
+    Mat dy = _dy.getMat();
+
+    if (templCenter == Point(-1, -1))
+        templCenter = Point(edges.cols / 2, edges.rows / 2);
+
+    setTemplateImpl(edges, dx, dy, templCenter);
+}
+
+void cv::GeneralizedHough::detect(InputArray _image, OutputArray positions, OutputArray votes, int cannyThreshold)
+{
+    Mat image = _image.getMat();
+
+    CV_Assert(image.type() == CV_8UC1);
+    CV_Assert(cannyThreshold > 0);
+
+    Canny(image, edges_, cannyThreshold / 2, cannyThreshold);
+    Sobel(image, dx_, CV_32F, 1, 0);
+    Sobel(image, dy_, CV_32F, 0, 1);
+
+    detectImpl(edges_, dx_, dy_, positions, votes);
+}
+
+void cv::GeneralizedHough::detect(InputArray _edges, InputArray _dx, InputArray _dy, OutputArray positions, OutputArray votes)
+{
+    cv::Mat edges = _edges.getMat();
+    cv::Mat dx = _dx.getMat();
+    cv::Mat dy = _dy.getMat();
+
+    detectImpl(edges, dx, dy, positions, votes);
+}
+
+void cv::GeneralizedHough::release()
+{
+    edges_.release();
+    dx_.release();
+    dy_.release();
+    releaseImpl();
+}

samples/cpp/generalized_hough.cpp

+#include <vector>
+#include <iostream>
+#include <string>
+
+#include "opencv2/core/core.hpp"
+#include "opencv2/imgproc/imgproc.hpp"
+#include "opencv2/gpu/gpu.hpp"
+#include "opencv2/highgui/highgui.hpp"
+#include "opencv2/contrib/contrib.hpp"
+
+using namespace std;
+using namespace cv;
+using namespace cv::gpu;
+
+static Mat loadImage(const string& name)
+{
+    Mat image = imread(name, IMREAD_GRAYSCALE);
+    if (image.empty())
+    {
+        cerr << "Can't load image - " << name << endl;
+        exit(-1);
+    }
+    return image;
+}
+
+int main(int argc, const char* argv[])
+{
+    CommandLineParser cmd(argc, argv,
+        "{ image i        | pic1.png  | input image }"
+        "{ template t     | templ.png | template image }"
+        "{ scale s        |           | estimate scale }"
+        "{ rotation r     |           | estimate rotation }"
+        "{ gpu            |           | use gpu version }"
+        "{ minDist        | 100       | minimum distance between the centers of the detected objects }"
+        "{ levels         | 360       | R-Table levels }"
+        "{ votesThreshold | 30        | the accumulator threshold for the template centers at the detection stage. The smaller it is, the more false positions may be detected }"
+        "{ angleThresh    | 10000     | angle votes treshold }"
+        "{ scaleThresh    | 1000      | scale votes treshold }"
+        "{ posThresh      | 100       | position votes threshold }"
+        "{ dp             | 2         | inverse ratio of the accumulator resolution to the image resolution }"
+        "{ minScale       | 0.5       | minimal scale to detect }"
+        "{ maxScale       | 2         | maximal scale to detect }"
+        "{ scaleStep      | 0.05      | scale step }"
+        "{ minAngle       | 0         | minimal rotation angle to detect in degrees }"
+        "{ maxAngle       | 360       | maximal rotation angle to detect in degrees }"
+        "{ angleStep      | 1         | angle step in degrees }"
+        "{ maxSize        | 1000      | maximal size of inner buffers }"
+        "{ help h ?       |           | print help message }"
+    );
+
+    cmd.about("This program demonstrates arbitary object finding with the Generalized Hough transform.");
+
+    if (cmd.has("help"))
+    {
+        cmd.printMessage();
+        return 0;
+    }
+
+    const string templName = cmd.get<string>("template");
+    const string imageName = cmd.get<string>("image");
+    const bool estimateScale = cmd.has("scale");
+    const bool estimateRotation = cmd.has("rotation");
+    const bool useGpu = cmd.has("gpu");
+    const double minDist = cmd.get<double>("minDist");
+    const int levels = cmd.get<int>("levels");
+    const int votesThreshold = cmd.get<int>("votesThreshold");
+    const int angleThresh = cmd.get<int>("angleThresh");
+    const int scaleThresh = cmd.get<int>("scaleThresh");
+    const int posThresh = cmd.get<int>("posThresh");
+    const double dp = cmd.get<double>("dp");
+    const double minScale = cmd.get<double>("minScale");
+    const double maxScale = cmd.get<double>("maxScale");
+    const double scaleStep = cmd.get<double>("scaleStep");
+    const double minAngle = cmd.get<double>("minAngle");
+    const double maxAngle = cmd.get<double>("maxAngle");
+    const double angleStep = cmd.get<double>("angleStep");
+    const int maxSize = cmd.get<int>("maxSize");
+
+    if (!cmd.check())
+    {
+        cmd.printErrors();
+        return -1;
+    }
+
+    Mat templ = loadImage(templName);
+    Mat image = loadImage(imageName);
+
+    int method = GHT_POSITION;
+    if (estimateScale)
+        method += GHT_SCALE;
+    if (estimateRotation)
+        method += GHT_ROTATION;
+
+    vector<Vec4f> position;
+    cv::TickMeter tm;
+
+    if (useGpu)
+    {
+        GpuMat d_templ(templ);
+        GpuMat d_image(image);
+        GpuMat d_position;
+
+        Ptr<GeneralizedHough_GPU> d_hough = GeneralizedHough_GPU::create(method);
+        d_hough->set("minDist", minDist);
+        d_hough->set("levels", levels);
+        d_hough->set("dp", dp);
+        d_hough->set("maxSize", maxSize);
+        if (estimateScale && estimateRotation)
+        {
+            d_hough->set("angleThresh", angleThresh);
+            d_hough->set("scaleThresh", scaleThresh);
+            d_hough->set("posThresh", posThresh);
+        }
+        else
+        {
+            d_hough->set("votesThreshold", votesThreshold);
+        }
+        if (estimateScale)
+        {
+            d_hough->set("minScale", minScale);
+            d_hough->set("maxScale", maxScale);
+            d_hough->set("scaleStep", scaleStep);
+        }
+        if (estimateRotation)
+        {
+            d_hough->set("minAngle", minAngle);
+            d_hough->set("maxAngle", maxAngle);
+            d_hough->set("angleStep", angleStep);
+        }
+
+        d_hough->setTemplate(d_templ);
+
+        tm.start();
+
+        d_hough->detect(d_image, d_position);
+        d_hough->download(d_position, position);
+
+        tm.stop();
+    }
+    else
+    {
+        Ptr<GeneralizedHough> hough = GeneralizedHough::create(method);
+        hough->set("minDist", minDist);
+        hough->set("levels", levels);
+        hough->set("dp", dp);
+        if (estimateScale && estimateRotation)
+        {
+            hough->set("angleThresh", angleThresh);
+            hough->set("scaleThresh", scaleThresh);
+            hough->set("posThresh", posThresh);
+            hough->set("maxSize", maxSize);
+        }
+        else
+        {
+            hough->set("votesThreshold", votesThreshold);
+        }
+        if (estimateScale)
+        {
+            hough->set("minScale", minScale);
+            hough->set("maxScale", maxScale);
+            hough->set("scaleStep", scaleStep);
+        }
+        if (estimateRotation)
+        {
+            hough->set("minAngle", minAngle);
+            hough->set("maxAngle", maxAngle);
+            hough->set("angleStep", angleStep);
+        }
+
+        hough->setTemplate(templ);
+
+        tm.start();
+
+        hough->detect(image, position);
+
+        tm.stop();
+    }
+
+    cout << "Found : " << position.size() << " objects" << endl;
+    cout << "Detection time : " << tm.getTimeMilli() << " ms" << endl;
+
+    Mat out;
+    cvtColor(image, out, COLOR_GRAY2BGR);
+
+    for (size_t i = 0; i < position.size(); ++i)
+    {
+        Point2f pos(position[i][0], position[i][1]);
+        float scale = position[i][2];
+        float angle = position[i][3];
+
+        RotatedRect rect;
+        rect.center = pos;
+        rect.size = Size2f(templ.cols * scale, templ.rows * scale);
+        rect.angle = angle;
+
+        Point2f pts[4];
+        rect.points(pts);
+
+        line(out, pts[0], pts[1], Scalar(0, 0, 255), 3);
+        line(out, pts[1], pts[2], Scalar(0, 0, 255), 3);
+        line(out, pts[2], pts[3], Scalar(0, 0, 255), 3);
+        line(out, pts[3], pts[0], Scalar(0, 0, 255), 3);
+    }
+
+    imshow("out", out);
+    waitKey();
+
+    return 0;
+}