diff --git a/modules/ximgproc/include/fast_hough_transform.hpp b/modules/ximgproc/include/fast_hough_transform.hpp new file mode 100644 index 000000000..d808eb060 --- /dev/null +++ b/modules/ximgproc/include/fast_hough_transform.hpp @@ -0,0 +1,165 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2015, Smart Engines Ltd, all rights reserved. +// Copyright (C) 2015, Institute for Information Transmission Problems of the Russian Academy of Sciences (Kharkevich Institute), all rights reserved. +// Copyright (C) 2015, Dmitry Nikolaev, Simon Karpenko, Michail Aliev, Elena Kuznetsova, 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 the copyright holders 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*/ + +#ifndef __OPENCV_FAST_HOUGH_TRANSFORM_HPP__ +#define __OPENCV_FAST_HOUGH_TRANSFORM_HPP__ +#ifdef __cplusplus + +#include "opencv2/core.hpp" + + +namespace cv { namespace ximgproc { + +/** +* @brief Specifies the part of Hough space to calculate +* @details The enum specifies the part of Hough space to calculate. Each +* member specifies primarily direction of lines (horizontal or vertical) +* and the direction of angle changes. +* Direction of angle changes is from multiples of 90 to odd multiples of 45. +* The image considered to be written top-down and left-to-right. +* Angles are started from vertical line and go clockwise. +* Separate quarters and halves are written in orientation they should be in +* full Hough space. +*/ +enum AngleRangeOption +{ + ARO_0_45 = 0, //< Vertical primarily direction and clockwise angle changes + ARO_45_90 = 1, //< Horizontal primarily direction and counterclockwise angle changes + ARO_90_135 = 2, //< Horizontal primarily direction and clockwise angle changes + ARO_315_0 = 3, //< Vertical primarily direction and counterclockwise angle changes + ARO_315_45 = 4, //< Vertical primarily direction + ARO_45_135 = 5, //< Horizontal primarily direction + ARO_315_135 = 6, //< Full set of directions + ARO_CTR_HOR = 7, //< 90 +/- atan(0.5), interval approximately from 64.5 to 116.5 degrees. + //< It is used for calculating Fast Hough Transform for images skewed by atan(0.5). + ARO_CTR_VER = 8 //< +/- atan(0.5), interval approximately from 333.5(-26.5) to 26.5 degrees + //< It is used for calculating Fast Hough Transform for images skewed by atan(0.5). +}; + +/** + * @brief Specifies binary operations. + * @details The enum specifies binary operations, that is such ones which involve + * two operands. Formally, a binary operation @f$ f @f$ on a set @f$ S @f$ + * is a binary relation that maps elements of the Cartesian product + * @f$ S \times S @f$ to @f$ S @f$: +* @f[ f: S \times S \to S @f] + * @ingroup MinUtils_MathOper + */ +enum HoughOp +{ + FHT_MIN = 0, //< Binary minimum operation. The constant specifies the binary minimum operation + //< @f$ f @f$ that is defined as follows: @f[ f(x, y) = \min(x, y) @f] + FHT_MAX = 1, //< Binary maximum operation. The constant specifies the binary maximum operation + //< @f$ f @f$ that is defined as follows: @f[ f(x, y) = \max(x, y) @f] + FHT_ADD = 2, //< Binary addition operation. The constant specifies the binary addition operation + //< @f$ f @f$ that is defined as follows: @f[ f(x, y) = x + y @f] + FHT_AVE = 3 //< Binary average operation. The constant specifies the binary average operation + //< @f$ f @f$ that is defined as follows: @f[ f(x, y) = \frac{x + y}{2} @f] +}; + +/** +* @brief Specifies to do or not to do skewing of Hough transform image +* @details The enum specifies to do or not to do skewing of Hough transform image +* so it would be no cycling in Hough transform image through borders of image. +*/ +enum HoughDeskewOption +{ + HDO_RAW = 0, //< Use raw cyclic image + HDO_DESKEW = 1 //< Prepare deskewed image +}; + +/** + * @brief Specifies the degree of rules validation. + * @details The enum specifies the degree of rules validation. This can be used, + * for example, to choose a proper way of input arguments validation. + */ +typedef enum { + RO_STRICT = 0x00, ///< Validate each rule in a proper way. + RO_IGNORE_BORDERS = 0x01, ///< Skip validations of image borders. +} RulesOption; + +/** +* @brief Calculates 2D Fast Hough transform of an image. +* @param dst The destination image, result of transformation. +* @param src The source (input) image. +* @param dstMatDepth The depth of destination image +* @param operation The operation to be applied, see cv::HoughOp +* @param angleRange The part of Hough space to calculate, see cv::AngleRangeOption +* @param makeSkew Specifies to do or not to do image skewing, see cv::HoughDeskewOption +* +* The function calculates the fast Hough transform for full, half or quarter +* range of angles. +*/ +CV_EXPORTS_W void FastHoughTransform( InputArray src, + OutputArray dst, + int dstMatDepth, + int angleRange = ARO_315_135, + int op = FHT_ADD, + int makeSkew = HDO_DESKEW ); + +/** +* @brief Calculates coordinates of line segment corresponded by point in Hough space. +* @param line Coordinates of line segment corresponded by point in Hough space. +* @param houghPoint Point in Hough space. +* @param srcImgInfo The source (input) image of Hough transform. +* @param angleRange The part of Hough space where point is situated, see cv::AngleRangeOption +* @param makeSkew Specifies to do or not to do image skewing, see cv::HoughDeskewOption +* @param rules Specifies strictness of line segment calculating, see cv::RulesOption +* @remarks If rules parameter set to RO_STRICT + then returned line cut along the border of source image. +* @remarks If rules parameter set to RO_WEAK then in case of point, which belongs + the incorrect part of Hough image, returned line will not intersect source image. +* +* The function calculates coordinates of line segment corresponded by point in Hough space. +*/ +CV_EXPORTS_W void HoughPoint2Line( OutputArray line, + const Point &houghPoint, + const Mat &srcImgInfo, + int angleRange = ARO_315_135, + int makeSkew = HDO_DESKEW, + int rules = RO_IGNORE_BORDERS ); + +} }// namespace cv::ximgproc + +#endif //__cplusplus +#endif //__OPENCV_FAST_HOUGH_TRANSFORM_HPP__ diff --git a/modules/ximgproc/perf/perf_fast_hough_transform.cpp b/modules/ximgproc/perf/perf_fast_hough_transform.cpp new file mode 100644 index 000000000..63cd04bca --- /dev/null +++ b/modules/ximgproc/perf/perf_fast_hough_transform.cpp @@ -0,0 +1,91 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2015, Smart Engines Ltd, all rights reserved. +// Copyright (C) 2015, Institute for Information Transmission Problems of the Russian Academy of Sciences (Kharkevich Institute), all rights reserved. +// Copyright (C) 2015, Dmitry Nikolaev, Simon Karpenko, Michail Aliev, Elena Kuznetsova, 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 the copyright holders 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 "perf_precomp.hpp" +#include "fast_hough_transform.hpp" + +namespace cvtest { + +using namespace std; +using namespace cv; +using namespace cv::ximgproc; +using namespace perf; +using namespace testing; + +using std::tr1::make_tuple; +using std::tr1::get; + +typedef std::tr1::tuple srcSize_srcType_dstDepth_t; +typedef perf::TestBaseWithParam + srcSize_srcType_dstDepth; + +#define ALL_MAT_DEPHTS CV_8U, CV_8S, CV_16U, CV_16S, CV_32S, CV_32F, CV_64F + +PERF_TEST_P(srcSize_srcType_dstDepth, FastHoughTransform, + testing::Combine( + testing::Values(TYPICAL_MAT_SIZES), + testing::Values(TYPICAL_MAT_TYPES), + testing::Values(ALL_MAT_DEPHTS) + ) + ) +{ + Size srcSize = get<0>(GetParam()); + int srcType = get<1>(GetParam()); + int dstDepth = get<2>(GetParam()); + + Mat src(srcSize, srcType); + Mat fht; + + declare.in(src, WARMUP_RNG); + + TEST_CYCLE_N(3) + { + FastHoughTransform(src, fht, dstDepth); + } + + SANITY_CHECK_NOTHING(); +} + +#undef ALL_MAT_DEPHTS + +} // namespace cvtest diff --git a/modules/ximgproc/samples/fast_hough_transform.cpp b/modules/ximgproc/samples/fast_hough_transform.cpp new file mode 100644 index 000000000..d4eebdd6f --- /dev/null +++ b/modules/ximgproc/samples/fast_hough_transform.cpp @@ -0,0 +1,333 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2015, Smart Engines Ltd, all rights reserved. +// Copyright (C) 2015, Institute for Information Transmission Problems of the Russian Academy of Sciences (Kharkevich Institute), all rights reserved. +// Copyright (C) 2015, Dmitry Nikolaev, Simon Karpenko, Michail Aliev, Elena Kuznetsova, 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 the copyright holders 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 +#include +#include + +#include + +#include "fast_hough_transform.hpp" + +#include +#include +#include +#include +#include + +using namespace cv; +using namespace cv::ximgproc; +using namespace std; + +static void help() +{ + cout << "\nThis program demonstrates line finding with the Fast Hough transform.\n" + "Usage:\n" + "./fasthoughtransform\n" + ", default is '../../../samples/data/building.jpg'\n" + ", default is " << CV_32S << "\n" + ", default is " << 6 << " (@see cv::AngleRangeOption)\n" + ", default is " << 2 << " (@see cv::HoughOp)\n" + ", default is " << 1 << "(@see cv::HoughDeskewOption)" << endl; +} + +static bool parseArgs(int argc, const char **argv, + Mat &img, + int &houghDepth, + int &houghAngleRange, + int &houghOperator, + int &houghSkew) +{ + if (argc > 6) + { + cout << "Too many arguments" << endl; + return false; + } + + const char *filename = argc >= 2 ? argv[1] + : "../../../samples/data/building.jpg"; + img = imread(filename, 0); + if (img.empty()) + { + cout << "Failed to load image from '" << filename << "'" << endl; + return false; + } + + houghDepth = argc >= 3 ? atoi(argv[2]) : CV_32S; + houghAngleRange = argc >= 4 ? atoi(argv[3]) : 6;//ARO_315_135 + houghOperator = argc >= 5 ? atoi(argv[4]) : 2;//FHT_ADD + houghSkew = argc >= 6 ? atoi(argv[5]) : 1;//HDO_DESKEW + + return true; +} + +static bool getEdges(const Mat &src, Mat &dst) +{ + Mat ucharSingleSrc; + src.convertTo(ucharSingleSrc, CV_8UC1); + + Canny(ucharSingleSrc, dst, 50, 200, 3); + return true; +} + +static bool fht(const Mat &src, Mat &dst, + int dstDepth, int angleRange, int op, int skew) +{ + clock_t clocks = clock(); + //for (int i = 0; i < 100; ++i) + //{ + FastHoughTransform(src, dst, dstDepth, angleRange, op, skew); + //} + clocks = clock() - clocks; + double secs = (double)clocks / CLOCKS_PER_SEC; + cout << std::setprecision(2) << "FastHoughTransform finished in " << secs + << " seconds" << endl; + + //clocks = clock(); + //for (int i = 0; i < 100; ++i) + //{ + // vector s_lines; + // HoughLines(src, s_lines, 1, CV_PI/360, 100, 0, 0); + //} + //clocks = clock() - clocks; + //secs = (double)clocks / CLOCKS_PER_SEC; + //std::cout << std::setprecision(2) << "HoughLines finished in " << secs << " seconds" << std::endl; + + return true; +} + +template +bool rel(pair const &a, pair const &b) +{ + return a.first > b.first; +} + +template +bool incIfGreater(const T& a, const T& b, int *value) +{ + if (!value || a < b) + return false; + if (a > b) + ++(*value); + return true; +} + +static const int MAX_LEN = 10000; + +template +bool getLocalExtr(vector &lines, + const Mat &src, + const Mat &fht, + float minWeight, + int maxCount) +{ + vector > weightedPoints; + for (int y = 0; y < fht.rows; ++y) + { + if (weightedPoints.size() > MAX_LEN) + break; + + T const *pLine = (T *)fht.ptr(max(y - 1, 0)); + T const *cLine = (T *)fht.ptr(y); + T const *nLine = (T *)fht.ptr(min(y + 1, fht.rows - 1)); + + for (int x = 0; x < fht.cols; ++x) + { + if (weightedPoints.size() > MAX_LEN) + break; + + T const value = cLine[x]; + if (value >= minWeight) + { + int isLocalMax = 0; + for (int xx = max(x - 1, 0); + xx <= min(x + 1, fht.cols - 1); + ++xx) + { + if (!incIfGreater(value, pLine[xx], &isLocalMax) || + !incIfGreater(value, cLine[xx], &isLocalMax) || + !incIfGreater(value, nLine[xx], &isLocalMax)) + { + isLocalMax = 0; + break; + } + } + if (isLocalMax > 0) + weightedPoints.push_back(make_pair(value, Point(x, y))); + } + } + } + + if (weightedPoints.empty()) + return true; + + sort(weightedPoints.begin(), weightedPoints.end(), &rel); + weightedPoints.resize(min(static_cast(weightedPoints.size()), + maxCount)); + + for (size_t i = 0; i < weightedPoints.size(); ++i) + { + Vec4i houghLine(0, 0, 0, 0); + HoughPoint2Line(houghLine, weightedPoints[i].second, src); + lines.push_back(houghLine); + } + return true; +} + +static bool getLocalExtr(vector &lines, + const Mat &src, + const Mat &fht, + float minWeight, + int maxCount) +{ + int const depth = CV_MAT_DEPTH(fht.type()); + switch (depth) + { + case 0: + return getLocalExtr(lines, src, fht, minWeight, maxCount); + case 1: + return getLocalExtr(lines, src, fht, minWeight, maxCount); + case 2: + return getLocalExtr(lines, src, fht, minWeight, maxCount); + case 3: + return getLocalExtr(lines, src, fht, minWeight, maxCount); + case 4: + return getLocalExtr(lines, src, fht, minWeight, maxCount); + case 5: + return getLocalExtr(lines, src, fht, minWeight, maxCount); + case 6: + return getLocalExtr(lines, src, fht, minWeight, maxCount); + default: + return false; + } +} + +static void rescale(Mat const &src, Mat &dst, + int const maxHeight=500, + int const maxWidth = 1000) +{ + double scale = min(min(static_cast(maxWidth) / src.cols, + static_cast(maxHeight) / src.rows), 1.0); + resize(src, dst, Size(), scale, scale, INTER_LINEAR); +} + +static void showHumanReadableImg(string const &name, Mat const &img) +{ + Mat ucharImg; + img.convertTo(ucharImg, CV_MAKETYPE(CV_8U, img.channels())); + rescale(ucharImg, ucharImg); + imshow(name, ucharImg); +} + +static void showFht(Mat const &fht) +{ + double minv(0), maxv(0); + minMaxLoc(fht, &minv, &maxv); + Mat ucharFht; + fht.convertTo(ucharFht, CV_MAKETYPE(CV_8U, fht.channels()), + 255.0 / (maxv + minv), minv / (maxv + minv)); + rescale(ucharFht, ucharFht); + imshow("fast hough transform", ucharFht); +} + +static void showLines(Mat const &src, vector const &lines) +{ + Mat bgrSrc; + cvtColor(src, bgrSrc, COLOR_GRAY2BGR); + + for (size_t i = 0; i < lines.size(); ++i) + { + Vec4i const &l = lines[i]; + line(bgrSrc, Point(l[0], l[1]), Point(l[2], l[3]), + Scalar(0, 0, 255), 1, LINE_AA); + } + + rescale(bgrSrc, bgrSrc); + imshow("lines", bgrSrc); +} + +int main(int argc, const char **argv) +{ + Mat src; + int depth(0); + int angleRange(0); + int op(0); + int skew(0); + + if (!parseArgs(argc, argv, src, depth, angleRange, op, skew)) + { + help(); + return -1; + } + showHumanReadableImg("src", src); + + Mat canny; + if (!getEdges(src, canny)) + { + cout << "Failed to select canny edges"; + return -2; + } + showHumanReadableImg("canny", canny); + + Mat hough; + if (!fht(canny, hough, depth, angleRange, op, skew)) + { + cout << "Failed to compute Fast Hough Transform"; + return -2; + } + showFht(hough); + + vector lines; + if (!getLocalExtr(lines, canny, hough, + static_cast(255 * 0.3 * min(src.rows, src.cols)), + 50)) + { + cout << "Failed to find local maximums on FHT image"; + return -2; + } + showLines(canny, lines); + + waitKey(); + + return 0; +} diff --git a/modules/ximgproc/src/fast_hough_transform.cpp b/modules/ximgproc/src/fast_hough_transform.cpp new file mode 100644 index 000000000..52abb773c --- /dev/null +++ b/modules/ximgproc/src/fast_hough_transform.cpp @@ -0,0 +1,1031 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2015, Smart Engines Ltd, all rights reserved. +// Copyright (C) 2015, Institute for Information Transmission Problems of the Russian Academy of Sciences (Kharkevich Institute), all rights reserved. +// Copyright (C) 2015, Dmitry Nikolaev, Simon Karpenko, Michail Aliev, Elena Kuznetsova, 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 the copyright holders 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" +#include "fast_hough_transform.hpp" + +namespace cv { namespace ximgproc { + +#if defined(_WIN32) && !defined(int32_t) + typedef __int32 int32_t; +#endif + +template +struct HoughOperator { }; +#define SPECIALIZE_HOUGHOP(TOp, body) \ + template \ + struct HoughOperator { \ + static void operate(T *pDst, T *pSrc0, T* pSrc1, int len) { \ + Mat dst (Size(1, len), D, pDst); \ + Mat src0(Size(1, len), D, pSrc0); \ + Mat src1(Size(1, len), D, pSrc1); \ + body; \ + } \ + }; +SPECIALIZE_HOUGHOP(FHT_ADD, add(src0, src1, dst)); +SPECIALIZE_HOUGHOP(FHT_MIN, min(src0, src1, dst)); +SPECIALIZE_HOUGHOP(FHT_MAX, max(src0, src1, dst)); +SPECIALIZE_HOUGHOP(FHT_AVE, addWeighted(src0, 0.5, src1, 0.5, 0.0, dst)); +#undef SPECIALIZE_HOUGHOP + +//----------------------fht---------------------------------------------------- + +template +void fhtCore(Mat &img0, + Mat &img1, + int32_t y0, + int32_t h, + bool isPositiveShift, + int level, + double aspl) +{ + if (level <= 0) + return; + + CV_Assert(h > 0); + if (h == 1) + { + if ((aspl != 0.0) && (level == 1)) + { + int w = img0.cols; + uchar* pLine0 = img0.data + img0.step * y0; + uchar* pLine1 = img1.data + img1.step * y0; + int dLine = cvRound(y0 * aspl); + dLine = dLine % w; + dLine = dLine * (int)(img1.elemSize()); + int wLine = img0.cols * (int)(img0.elemSize()); + memcpy(pLine0, pLine1 + wLine - dLine, dLine); + memcpy(pLine0 + dLine, pLine1, wLine - dLine); + } + else + { + memcpy(img0.data + img0.step * y0, + img1.data + img1.step * y0, + img0.cols * (int)(img0.elemSize())); + } + return; + } + const int32_t k = h >> 1; + fhtCore(img1, img0, y0, k, + isPositiveShift, level - 1, aspl); + fhtCore(img1, img0, y0 + k, h - k, + isPositiveShift, level - 1, aspl); + + int au = 2 * k - 2; + int ad = 2 * h - 2 * k - 2; + int b = h - 1; + int d = 2 * h - 2; + int w = img0.cols; + int wm = (h / w + 1) * w; + + for (int32_t s = 0; s < h; s++) + { + int su = (s * au + b) / d; + int sd = (s * ad + b) / d; + int rd = isPositiveShift ? sd - s : s - sd; + rd = (rd + wm) % w; + uchar *pLine0 = img0.data + img0.step * (y0 + s); + uchar *pLineU = img1.data + img1.step * (y0 + su); + uchar *pLineD = img1.data + img1.step * (y0 + k + sd); + int w0 = img0.channels() * rd; + int w1 = img0.channels() * (w - rd); + + if ((aspl != 0.0) && (level == 1)) + { + int dU = cvRound((y0 + su) * aspl); + dU = dU % w; + dU *= img0.channels(); + int dD = cvRound((y0 + k + sd) * aspl); + dD = dD % w; + dD *= img0.channels(); + int wB = w * img0.channels(); + + int dX = dD - dU; + if (w0 >= dX) + { + if (w0 >= dD) + { + HoughOperator::operate((T *)pLine0 + dU, + (T *)pLineU, + (T *)pLineD + (w0 - dX), + w1 + dX); + HoughOperator::operate((T *)pLine0 + (w1 + dD), + (T *)pLineU + (w1 + dX), + (T *)pLineD, + w0 - dD); + HoughOperator::operate((T *)pLine0, + (T *)pLineU + (wB - dU), + (T *)pLineD + (w0 - dD), + dU); + } + else + { + HoughOperator::operate((T *)pLine0 + dU, + (T *)pLineU, + (T *)pLineD + (w0 - dX), + wB - dU); + HoughOperator::operate((T *)pLine0, + (T *)pLineU + (wB - dU), + (T *)pLineD + (w0 + wB - dD), + dD - w0); + HoughOperator::operate((T *)pLine0 + (dD - w0), + (T *)pLineU + (w1 + dX), + (T *)pLineD, + w0 - dX); + } + } + else + { + HoughOperator::operate((T *)pLine0 + dU, + (T *)pLineU, + (T *)pLineD + (wB - (dX - w0)), + dX - w0); + HoughOperator::operate((T *)pLine0 + (dD - w0), + (T *)pLineU + (dX - w0), + (T *)pLineD, + wB - (dX - w0) - dU); + HoughOperator::operate((T *)pLine0, + (T *)pLineU + (wB - dU), + (T *)pLineD + (wB - (dX - w0) - dU), + dU); + } + } + else + { + HoughOperator::operate((T *)pLine0, + (T *)pLineU, + (T *)pLineD + w0, + w1); + HoughOperator::operate((T *)pLine0 + w1, + (T *)pLineU + w1, + (T *)pLineD, + w0); + } + } +} + +template +void fhtVoT(Mat &img0, + Mat &img1, + bool isPositiveShift, + double aspl) +{ + int level = 0; + for (int thres = 1; img0.rows > thres; thres <<= 1) + level++; + + fhtCore(img0, img1, 0, img0.rows, isPositiveShift, level, aspl); +} + +template +void fhtVo(Mat &img0, + Mat &img1, + bool isPositiveShift, + int operation, + double aspl) +{ + switch (operation) + { + case FHT_ADD: + fhtVoT(img0, img1, isPositiveShift, aspl); + break; + case FHT_AVE: + fhtVoT(img0, img1, isPositiveShift, aspl); + break; + case FHT_MAX: + fhtVoT(img0, img1, isPositiveShift, aspl); + break; + case FHT_MIN: + fhtVoT(img0, img1, isPositiveShift, aspl); + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown operation %d", operation)); + break; + } +} + +static void fhtVo(Mat &img0, + Mat &img1, + bool isPositiveShift, + int operation, + double aspl) +{ + int const depth = img0.depth(); + switch (depth) + { + case CV_8U: + fhtVo(img0, img1, isPositiveShift, operation, aspl); + break; + case CV_8S: + fhtVo(img0, img1, isPositiveShift, operation, aspl); + break; + case CV_16U: + fhtVo(img0, img1, isPositiveShift, operation, aspl); + break; + case CV_16S: + fhtVo(img0, img1, isPositiveShift, operation, aspl); + break; + case CV_32S: + fhtVo(img0, img1, isPositiveShift, operation, aspl); + break; + case CV_32F: + fhtVo(img0, img1, isPositiveShift, operation, aspl); + break; + case CV_64F: + fhtVo(img0, img1, isPositiveShift, operation, aspl); + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown depth %d", depth)); + break; + } +} + +static void FHT(Mat &dst, + const Mat &src, + int operation, + bool isVertical, + bool isClockwise, + double aspl) +{ + CV_Assert(dst.cols > 0 && dst.rows > 0); + CV_Assert(src.channels() == dst.channels()); + if (isVertical) + CV_Assert(src.cols == dst.cols && src.rows == dst.rows); + else + CV_Assert(src.cols == dst.rows && src.rows == dst.cols); + + int level = 0; + for (int thres = 1; dst.rows > thres; thres <<= 1) + level++; + + Mat tmp; + src.convertTo(tmp, dst.type()); + if (!isVertical) + transpose(tmp, tmp); + tmp.copyTo(dst); + + fhtVo(dst, tmp, + isVertical ? isClockwise : !isClockwise, + operation, aspl); +} + +static void calculateFHTQuadrant(Mat &dst, + const Mat &src, + int operation, + int quadrant) +{ + bool bVert = true; + bool bClock = true; + double aspl = 0.0; + switch (quadrant) + { + case ARO_315_0: + bVert = true; + bClock = false; + break; + case ARO_0_45: + bVert = true; + bClock = true; + break; + case ARO_45_90: + bVert = false; + bClock = false; + break; + case ARO_90_135: + bVert = false; + bClock = true; + break; + case ARO_CTR_VER: + bVert = true; + bClock = false; + aspl = 0.5; + break; + case ARO_CTR_HOR: + bVert = false; + bClock = true; + aspl = 0.5; + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown quadrant %d", quadrant)); + } + + FHT(dst, src, operation, bVert, bClock, aspl); +} + +static void createDstFhtMat(OutputArray dst, + InputArray src, + int depth, + int angleRange) +{ + int const rows = src.size().height; + int const cols = src.size().width; + int const channels = src.channels(); + + int wd = cols + rows; + int ht = 0; + + switch (angleRange) + { + case ARO_315_0: + case ARO_0_45: + case ARO_CTR_VER: + ht = rows; + break; + case ARO_45_90: + case ARO_90_135: + case ARO_CTR_HOR: + ht = cols; + break; + case ARO_315_45: + ht = 2 * rows - 1; + break; + case ARO_45_135: + ht = 2 * cols - 1; + break; + case ARO_315_135: + ht = 2 * (cols + rows) - 3; + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown angleRange %d", angleRange)); + } + + dst.create(ht, wd, CV_MAKETYPE(depth, channels)); +} + +static void createFHTSrc(Mat &srcFull, + const Mat &src, + int angleRange) +{ + bool verticalTiling = false; + switch (angleRange) + { + case ARO_315_0: + case ARO_0_45: + case ARO_CTR_VER: + case ARO_315_45: + verticalTiling = false; + break; + case ARO_45_90: + case ARO_90_135: + case ARO_CTR_HOR: + case ARO_45_135: + verticalTiling = true; + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown angleRange %d", angleRange)); + } + + int wd = verticalTiling ? src.cols : src.cols + src.rows; + int ht = verticalTiling ? src.cols + src.rows : src.rows; + srcFull = Mat(ht, wd, src.type()); + + Mat imgReg; + if (verticalTiling) + imgReg = Mat(srcFull, Rect(0, src.rows, src.cols, src.cols)); + else + imgReg = Mat(srcFull, Rect(src.cols, 0, src.rows, src.rows)); + imgReg = Mat::zeros(imgReg.size(), imgReg.type()); + + imgReg = Mat(srcFull, Rect(0, 0, src.cols, src.rows)); + src.copyTo(imgReg); +} + +static void setFHTDstRegion(Mat &dstRegion, + const Mat &dst, + const Mat &src, + int quadrant, + int angleRange) +{ + int base = -1; + switch (angleRange) + { + case ARO_315_0: + case ARO_315_45: + case ARO_315_135: + base = 0; + break; + case ARO_0_45: + base = 1; + break; + case ARO_45_90: + case ARO_45_135: + base = 2; + break; + case ARO_90_135: + base = 3; + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown angleRange %d", angleRange)); + } + + int quad = -1; + switch (quadrant) + { + case ARO_315_0: + quad = 0; + break; + case ARO_0_45: + quad = 1; + break; + case ARO_45_90: + quad = 2; + break; + case ARO_90_135: + quad = 3; + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown quadrant %d", quadrant)); + } + + if (quad < base) + quad += 4; + + int shift = 0; + for (int i = base; i < quad; i++) + shift += (i & 2) ? src.cols - 1 : src.rows - 1; + const int ht = (quad & 2) ? src.cols : src.rows; + + dstRegion = Mat(dst, Rect(0, shift, src.rows + src.cols, ht)); +} + + + +static void rotateLineRightCyclic(uchar *pLine, + uchar *pBuf, + int len, + int shift) +{ + shift = shift % len; + shift = (shift + len) % len; + memcpy(pBuf, pLine, len); + memcpy(pLine + shift, pBuf, len - shift); + if (shift > 0) + memcpy(pLine, pBuf + len - shift, shift); +} + +static void skewQuadrant(Mat &quad, + const Mat &src, + uchar *pBuf, + int quadrant) +{ + CV_Assert(pBuf); + + const int wd = src.cols; + const int ht = src.rows; + + double start = 0.; + double step = .5; + switch (quadrant) + { + case ARO_315_0: + case ARO_CTR_VER: + step = -.5; + start = ht - 0.5; + break; + case ARO_0_45: + step = -.5; + start = ht * .5; + break; + case ARO_45_90: + break; + case ARO_90_135: + case ARO_CTR_HOR: + start = wd * .5 - 0.5; + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown quadrant %d", quadrant)); + } + + const int pixlen = static_cast(quad.elemSize()); + const int len = quad.cols * pixlen; + for (int y = 0; y < quad.rows; y++) + { + uchar *pLine = quad.ptr(y); + int shift = static_cast(start + step * y) * pixlen; + rotateLineRightCyclic(pLine, pBuf, len, shift); + } +} + +void FastHoughTransform(InputArray src, + OutputArray dst, + int dstMatDepth, + int angleRange, + int operation, + int makeSkew) +{ + Mat srcMat = src.getMat(); + if (!srcMat.isContinuous()) + srcMat = srcMat.clone(); + CV_Assert(srcMat.cols > 0 && srcMat.rows > 0); + + createDstFhtMat(dst, src, dstMatDepth, angleRange); + Mat dstMat = dst.getMat(); + + Mat imgRegDst; + const int len = dstMat.cols * static_cast(dstMat.elemSize()); + CV_Assert(len > 0); + std::vector buf_(len); + uchar *buf(&buf_[0]); + + if (angleRange == ARO_315_135) + { + { + Mat imgSrc; + createFHTSrc(imgSrc, srcMat, ARO_315_45); + + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_315_0, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_315_0); + flip(imgRegDst, imgRegDst, 0); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_315_0); + + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_0_45, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_0_45); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_0_45); + } + { + Mat imgSrc; + createFHTSrc(imgSrc, srcMat, ARO_45_135); + + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_45_90, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_45_90); + flip(imgRegDst, imgRegDst, 0); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_45_90); + + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_90_135, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_90_135); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_90_135); + } + return; + } + + Mat imgSrc; + createFHTSrc(imgSrc, srcMat, angleRange); + + switch (angleRange) + { + case ARO_315_0: + calculateFHTQuadrant(dstMat, imgSrc, operation, angleRange); + flip(dstMat, dstMat, 0); + if (HDO_DESKEW == makeSkew) + skewQuadrant(dstMat, imgSrc, buf, angleRange); + return; + case ARO_0_45: + calculateFHTQuadrant(dstMat, imgSrc, operation, angleRange); + if (HDO_DESKEW == makeSkew) + skewQuadrant(dstMat, imgSrc, buf, angleRange); + return; + case ARO_45_90: + calculateFHTQuadrant(dstMat, imgSrc, operation, angleRange); + flip(dstMat, dstMat, 0); + if (HDO_DESKEW == makeSkew) + skewQuadrant(dstMat, imgSrc, buf, angleRange); + return; + case ARO_90_135: + calculateFHTQuadrant(dstMat, imgSrc, operation, angleRange); + if (HDO_DESKEW == makeSkew) + skewQuadrant(dstMat, imgSrc, buf, angleRange); + return; + case ARO_315_45: + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_315_0, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_315_0); + flip(imgRegDst, imgRegDst, 0); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_315_0); + + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_0_45, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_0_45); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_0_45); + return; + case ARO_45_135: + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_45_90, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_45_90); + flip(imgRegDst, imgRegDst, 0); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_45_90); + + setFHTDstRegion(imgRegDst, dstMat, srcMat, ARO_90_135, angleRange); + calculateFHTQuadrant(imgRegDst, imgSrc, operation, ARO_90_135); + if (HDO_DESKEW == makeSkew) + skewQuadrant(imgRegDst, imgSrc, buf, ARO_90_135); + return; + case ARO_CTR_VER: + calculateFHTQuadrant(dstMat, imgSrc, operation, angleRange); + flip(dstMat, dstMat, 0); + if (HDO_DESKEW == makeSkew) + skewQuadrant(dstMat, imgSrc, buf, angleRange); + return; + case ARO_CTR_HOR: + calculateFHTQuadrant(dstMat, imgSrc, operation, angleRange); + if (HDO_DESKEW == makeSkew) + skewQuadrant(dstMat, imgSrc, buf, angleRange); + return; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown angleRange %d", angleRange)); + } +} + +//----------------------------------------------------------------------------- + +//----------------------fht point2line----------------------------------------- +struct LineSegment { + Point u, v; + LineSegment(const Point _u, const Point _v) : u(_u), v(_v) { } +}; + +static void getRawPoint(Point &rawHoughPoint, + int &quadRawPoint, + const Point &givenHoughPoint, + const Mat &srcImgInfo, + int angleRange, + int makeSkew) +{ + int base = -1; + switch (angleRange) + { + case ARO_315_0: + case ARO_315_45: + case ARO_CTR_VER: + case ARO_315_135: + base = 0; + break; + case ARO_0_45: + base = 1; + break; + case ARO_45_90: + case ARO_45_135: + base = 2; + break; + case ARO_90_135: + case ARO_CTR_HOR: + base = 3; + break; + default: + CV_Error_(CV_StsNotImplemented, ("Unknown angleRange %d", angleRange)); + } + + int const cols = srcImgInfo.cols; + int const rows = srcImgInfo.rows; + + rawHoughPoint.y = givenHoughPoint.y; + int quad = 0, qsize = 0; + for (quad = base; quad < 4; quad++) + { + qsize = (quad & 2) ? cols - 1 : rows - 1; + if (rawHoughPoint.y <= qsize) + break; + rawHoughPoint.y -= qsize; + } + if (quad >= 4) + CV_Error(CV_StsInternal, ""); + + quadRawPoint = quad; + + double skewShift = 0.0; + if (makeSkew == HDO_DESKEW) + { + switch (quad) + { + case 0: + skewShift = rows - (rawHoughPoint.y + 1) * 0.5; + break; + case 1: + skewShift = (rows - rawHoughPoint.y) * 0.5; + break; + case 2: + skewShift = rawHoughPoint.y * 0.5; + break; + default: + skewShift = 0.5 * (cols + rawHoughPoint.y - 1); + break; + } + } + + rawHoughPoint.x = givenHoughPoint.x - static_cast(skewShift); + if (rawHoughPoint.x < 0) + rawHoughPoint.x = rows + cols + rawHoughPoint.x; +} + +static bool checkRawPoint(const Point &rawHoughPoint, + int quadRawPoint, + const Mat &srcImgInfo) +{ + int const cols = srcImgInfo.cols; + int const rows = srcImgInfo.rows; + + switch (quadRawPoint) + { + case 0: + //down left triangle on FHT + if (rawHoughPoint.x - cols <= rawHoughPoint.y && + rawHoughPoint.x - cols >= 0) + return false; + break; + case 1: + //up right triangle on FHT + if (rawHoughPoint.x - cols >= rawHoughPoint.y) + return false; + break; + case 2: + //up right triangle on up-down FHT image + if (rawHoughPoint.x - rows >= cols - 1 - rawHoughPoint.y) + return false; + break; + default: + //down left triangle on up-down FHT image + if (rawHoughPoint.x - rows <= cols - 1 - rawHoughPoint.y && + rawHoughPoint.x - rows >= 0) + return false; + break; + } + return true; +} + +static void shiftLineSegment(LineSegment &segment, + const Point &shift) +{ + segment.u.x += shift.x; + segment.v.x += shift.x; + segment.u.y += shift.y; + segment.v.y += shift.y; +} + +static void lineFactors(double &a, + double &b, + double &c, + const Point &point1, + const Point &point2) +{ + CV_Assert(point1.x != point2.x || point1.y != point2.y); + + Point vectorSegment(point2.x - point1.x, point2.y - point1.y); + a = - vectorSegment.y; + b = vectorSegment.x; + c = - (a * point1.x + b * point1.y); +} + +static void crossSegments(Point &point, + const LineSegment &line1, + const LineSegment &line2) +{ + double a1 = 0.0, b1 = 0.0, c1 = 0.0; + double a2 = 0.0, b2 = 0.0, c2 = 0.0; + lineFactors(a1, b1, c1, line1.u, line1.v); + lineFactors(a2, b2, c2, line2.u, line2.v); + + double uLine1onLine2 = a2 * line1.u.x + b2 * line1.u.y + c2; + double vLine1onLine2 = a2 * line1.v.x + b2 * line1.v.y + c2; + + double ULine2onLine1 = a1 * line2.u.x + b1 * line2.u.y + c1; + double VLine2onLine1 = a1 * line2.v.x + b1 * line2.v.y + c1; + + CV_Assert(ULine2onLine1 != 0 || VLine2onLine1 != 0 || + uLine1onLine2 != 0 || vLine1onLine2 != 0); + CV_Assert(ULine2onLine1 * VLine2onLine1 <= 0 && + uLine1onLine2 * vLine1onLine2 <= 0); + + static const double double_eps = 1e-10; + CV_Assert(std::abs(uLine1onLine2 - vLine1onLine2) >= double_eps); + + double mul = uLine1onLine2 / (uLine1onLine2 - vLine1onLine2); + point.x = cvRound(line1.u.x + mul * (line1.v.x - line1.u.x)); + point.y = cvRound(line1.u.y + mul * (line1.v.y - line1.u.y)); +} + +void HoughPoint2Line(OutputArray line, + const Point &houghPoint, + const Mat &srcImgInfo, + int angleRange, + int makeSkew, + int rules) +{ + int const cols = srcImgInfo.cols; + int const rows = srcImgInfo.rows; + + CV_Assert(houghPoint.y >= 0); + CV_Assert(houghPoint.x < cols + rows); + + int quad = 0; + Point rawPoint(0, 0); + getRawPoint(rawPoint, quad, houghPoint, srcImgInfo, angleRange, makeSkew); + bool ret = checkRawPoint(rawPoint, quad, srcImgInfo); + if (!(rules & RO_IGNORE_BORDERS)) + { + CV_Assert(ret); + } + + LineSegment dstLine(Point(0, 0), Point(0, 0)); + switch (quad) + { + case 0: + dstLine.v.y = rows - 1; + dstLine.u.x = rawPoint.x; + dstLine.v.x = dstLine.u.x + rows - rawPoint.y - 1; + break; + case 1: + dstLine.v.y = rows - 1; + dstLine.u.x = rawPoint.x; + dstLine.v.x = dstLine.u.x - rawPoint.y; + break; + case 2: + dstLine.v.x = cols - 1; + dstLine.u.y = rawPoint.x; + dstLine.v.y = dstLine.u.y - cols + rawPoint.y + 1; + break; + default: + dstLine.v.x = cols - 1; + dstLine.u.y = rawPoint.x; + dstLine.v.y = dstLine.u.y + rawPoint.y; + break; + } + + if (angleRange == ARO_CTR_VER) + { + int shift = cvRound(0.5 * dstLine.u.y); + shift = shift % (cols + rows); + dstLine.u.x = dstLine.u.x - shift; + + shift = cvRound(0.5 * dstLine.v.y); + shift = shift % (cols + rows); + dstLine.v.x = dstLine.v.x - shift; + } + else if (angleRange == ARO_CTR_HOR) + { + int shift = cvRound(0.5 * dstLine.u.x); + shift = shift % (cols + rows); + dstLine.u.y = dstLine.u.y - shift; + + shift = cvRound(0.5 * dstLine.v.x); + shift = shift % (cols + rows); + dstLine.v.y = dstLine.v.y - shift; + } + + if (!ret) + { + Vec4i pts(dstLine.v.x, dstLine.v.y, dstLine.u.x, dstLine.u.y); + Mat(pts).copyTo(line); + return; + } + + if (rules & RO_IGNORE_BORDERS) + { + switch (quad) + { + case 0: + if (dstLine.v.x > cols + rows - 1) + { + Point shiftVector(-(cols + rows), 0); + shiftLineSegment(dstLine, shiftVector); + } + break; + case 3: + if (dstLine.v.y > rows + cols - 1) + { + Point shiftVector(0, -(cols + rows)); + shiftLineSegment(dstLine, shiftVector); + } + break; + default: + break; + } + + Vec4i pts(dstLine.v.x, dstLine.v.y, dstLine.u.x, dstLine.u.y); + Mat(pts).copyTo(line); + return; + } + + Point minIntersectPoint(0, 0); + + Point minLeftUpSrcPoint(0, 0); + Point minRightUpSrcPoint(cols - 1, 0); + Point minLeftDownSrcPoint(0, rows - 1); + Point minRightDownSrcPoint(cols - 1, rows - 1); + + switch (quad) + { + case 0: + if (dstLine.v.x > cols + rows - 1) + { + LineSegment minRightToDstLine(Point(cols + rows, 0), + Point(cols + rows, rows - 1) ); + crossSegments(minIntersectPoint, dstLine, minRightToDstLine); + dstLine.u.y = minIntersectPoint.y; + dstLine.u.x = 0; + dstLine.v.x = dstLine.v.x - (cols + rows); + } + if (dstLine.v.x > cols - 1) + { + LineSegment minRightSrcLine(minRightUpSrcPoint, + minRightDownSrcPoint); + crossSegments(minIntersectPoint, dstLine, minRightSrcLine); + dstLine.v.y = minIntersectPoint.y; + dstLine.v.x = cols - 1; + } + break; + case 1: + if (dstLine.v.x < 0) + { + LineSegment minLeftSrcLine(minLeftUpSrcPoint, minLeftDownSrcPoint); + crossSegments(minIntersectPoint, dstLine, minLeftSrcLine); + dstLine.v.y = minIntersectPoint.y; + dstLine.v.x = 0; + } + if (dstLine.u.x > cols - 1) + { + LineSegment minRightSrcLine(minRightUpSrcPoint, + minRightDownSrcPoint); + crossSegments(minIntersectPoint, dstLine, minRightSrcLine); + dstLine.u.y = minIntersectPoint.y; + dstLine.u.x = cols - 1; + } + break; + case 2: + if (dstLine.v.y < 0) + { + LineSegment minTopSrcLine(minLeftUpSrcPoint, minRightUpSrcPoint); + crossSegments(minIntersectPoint, dstLine, minTopSrcLine); + dstLine.v.x = minIntersectPoint.x; + dstLine.v.y = 0; + } + if (dstLine.u.y > rows - 1) + { + LineSegment minBottomSrcLine(minLeftDownSrcPoint, + minRightDownSrcPoint); + crossSegments(minIntersectPoint, dstLine, minBottomSrcLine ); + dstLine.u.x = minIntersectPoint.x; + dstLine.u.y = rows - 1; + } + break; + default: + if (dstLine.v.y > rows + cols - 1) + { + LineSegment minDownToDstLine(Point(0, rows + cols), + Point(cols - 1, rows + cols)); + crossSegments(minIntersectPoint, dstLine, minDownToDstLine); + dstLine.u.x = minIntersectPoint.x; + dstLine.u.y = 0; + dstLine.v.y = dstLine.v.y - (rows + cols); + } + if (dstLine.v.y > rows - 1) + { + LineSegment minBottomSrcLine(minLeftDownSrcPoint, + minRightDownSrcPoint); + crossSegments(minIntersectPoint, dstLine, minBottomSrcLine); + dstLine.v.x = minIntersectPoint.x; + dstLine.v.y = rows - 1; + } + break; + } + + Vec4i pts(dstLine.v.x, dstLine.v.y, dstLine.u.x, dstLine.u.y); + Mat(pts).copyTo(line); +} + +//----------------------------------------------------------------------------- + + +} } // namespace cv::ximgproc diff --git a/modules/ximgproc/test/test_fast_hough_transform.cpp b/modules/ximgproc/test/test_fast_hough_transform.cpp new file mode 100644 index 000000000..94819be4a --- /dev/null +++ b/modules/ximgproc/test/test_fast_hough_transform.cpp @@ -0,0 +1,472 @@ +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2015, Smart Engines Ltd, all rights reserved. +// Copyright (C) 2015, Institute for Information Transmission Problems of the Russian Academy of Sciences (Kharkevich Institute), all rights reserved. +// Copyright (C) 2015, Dmitry Nikolaev, Simon Karpenko, Michail Aliev, Elena Kuznetsova, 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 the copyright holders 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" +#include "fast_hough_transform.hpp" + +namespace cvtest +{ +using namespace cv; +using namespace cv::ximgproc; +using namespace std; +using namespace testing; + +using std::tr1::make_tuple; +using std::tr1::get; + +//----------------------utils--------------------------------------------------- + +template struct Eps +{ + static T get() { return 1; } +}; +template <> struct Eps { static float get() { return float(1e-3); } }; +template <> struct Eps { static double get() { return 1e-6; } }; + +template struct MinPos +{ + static T get() { return Eps::get(); } +}; + +template struct Max { static T get() +{ + return saturate_cast(numeric_limits::max()); } +}; + +template struct Rand +{ + static T get(T _min = MinPos::get(), T _max = Max::get()) + { + RNG& rng = TS::ptr()->get_rng(); + return saturate_cast(rng.uniform(int(std::max(MinPos::get(), + _min)), + int(std::min(Max::get(), + _max)))); + } +}; +template <> struct Rand +{ + static float get(float _min = MinPos::get(), + float _max = Max::get()) + { + RNG& rng = TS::ptr()->get_rng(); + return rng.uniform(std::max(MinPos::get(), _min), + std::min(Max::get(), _max)); + } +}; +template <> struct Rand +{ + static double get(double _min = MinPos::get(), + double _max = Max::get()) + { + RNG& rng = TS::ptr()->get_rng(); + return rng.uniform(std::max(MinPos::get(), _min), + std::min(Max::get(), _max)); + } +}; + +template struct Eq +{ + static bool get(T a, T b) + { + return a < b ? b - a < Eps::get() : a - b < Eps::get(); + } +}; + +//----------------------TestFHT------------------------------------------------- +class TestFHT +{ +public: + TestFHT() : ts(TS::ptr()) {} + + void run_n_tests(int depth, + int channels, + int pts_count, + int n_per_test); + +private: + template + int run_n_tests_t(int depth, + int channels, + int pts_count, + int n_per_test); + + template + int run_test(int depth, + int channels, + int pts_count); + + template + int put_random_points(Mat &img, + int count, + vector &pts); + + int run_func(Mat const&src, + Mat& fht); + + template + int validate_test_results(Mat const &fht, + Mat const &src, + vector const& pts); + + template int validate_sum(Mat const& src, Mat const& fht); + int validate_point(Mat const& fht, vector const &pts); + int validate_line(Mat const& fht, Mat const& src, vector const& pts); + +private: + TS *ts; +}; + +template +int TestFHT::put_random_points(Mat &img, int count, vector &pts) +{ + int code = TS::OK; + + pts.resize(count, Point(-1, -1)); + + for (int i = 0; i < count; ++i) + { + RNG rng = ts->get_rng(); + Point const pt(rng.uniform(0, img.cols), + rng.uniform(0, img.rows)); + pts[i] = pt; + + for (int c = 0; c < img.channels(); ++c) + { + T color = Rand::get(MinPos::get(), + T(Max::get() / count)); + + T *img_line = (T*)(img.data + img.step * pt.y); + img_line[pt.x * img.channels() + c] = color; + } + } + + return code; +} + +template +int TestFHT::validate_sum(Mat const& src, Mat const& fht) +{ + int const channels = src.channels(); + if (fht.channels() != channels) + return TS::FAIL_BAD_ARG_CHECK; + + + vector src_channels(channels); + split(src, src_channels); + vector fht_channels(channels); + split(fht, fht_channels); + + for (int c = 0; c < channels; ++c) + { + T const src_sum = saturate_cast(sum(src_channels[c]).val[0]); + for (int y = 0; y < fht.rows; ++y) + { + T const fht_sum = saturate_cast(sum(fht_channels[c].row(y)).val[0]); + if (!Eq::get(src_sum, fht_sum)) + { + ts->printf(TS::LOG, + "The sum of column #%d of channel #%d of the fast " + "hough transform result and the sum of source image" + " mismatch (=%g, should be =%g)\n", + y, c, (float)fht_sum, (float)src_sum); + return TS::FAIL_BAD_ACCURACY; + } + } + } + return TS::OK; +} + +int TestFHT::validate_point(Mat const& fht, + vector const &pts) +{ + if (pts.empty()) + return TS::OK; + + for (size_t i = 1; i < pts.size(); ++i) + { + if (pts[0] != pts[i]) + return TS::OK; + } + + int const channels = fht.channels(); + vector fht_channels(channels); + split(fht, fht_channels); + + for (int c = 0; c < channels; ++c) + { + for (int y = 0; y < fht.rows; ++y) + { + int cnt = countNonZero(fht_channels[c].row(y)); + if (cnt != 1) + { + ts->printf(TS::LOG, + "The incorrect count of non-zero values in column " + "#%d, channel #%d of FastHoughTransform result " + "image (=%d, should be %d)\n", + y, c, cnt, 1); + return TS::FAIL_BAD_ACCURACY; + } + } + } + return TS::OK; +} + +static const double MAX_LDIST = 2.0; +int TestFHT::validate_line(Mat const& fht, + Mat const& src, + vector const& pts) +{ + size_t const size = (int)pts.size(); + if (size < 2) + return TS::OK; + size_t first_pt_i = 0, second_pt_i = 1; + for (size_t i = first_pt_i + 1; i < size; ++i) + { + if (pts[i] != pts[first_pt_i]) + { + second_pt_i = first_pt_i; + break; + } + } + if (pts[second_pt_i] == pts[first_pt_i]) + return TS::OK; + for (size_t i = second_pt_i + 1; i < size; ++i) + { + if (pts[i] != pts[second_pt_i]) + return TS::OK; + } + const Point &f = pts[first_pt_i]; + const Point &s = pts[second_pt_i]; + + int const channels = fht.channels(); + vector fht_channels(channels); + split(fht, fht_channels); + + for (int ch = 0; ch < channels; ++ch) + { + Point fht_max(-1, -1); + minMaxLoc(fht_channels[ch], 0, 0, 0, &fht_max); + Vec4i src_line; + HoughPoint2Line(src_line, fht_max, src, + ARO_315_135, HDO_DESKEW, RO_STRICT); + + double const a = src_line[1] - src_line[3]; + double const b = src_line[2] - src_line[0]; + double const c = - (a * src_line[0] + b * src_line[1]); + + double const fd = abs(f.x * a + f.y * b + c) / sqrt(a * a + b * b); + double const sd = abs(s.x * a + s.y * b + c) / sqrt(a * a + b * b); + double const dist = std::max(fd, sd); + + if (dist > MAX_LDIST) + { + ts->printf(TS::LOG, + "Failed to detect max line in channels %d (distance " + "between point and line correspoinding of maximum in " + "FastHoughTransform space is #%g)\n", ch, dist); + return TS::FAIL_BAD_ACCURACY; + } + } + + return TS::OK; +} + +template +int TestFHT::validate_test_results(Mat const &fht, + Mat const &src, + vector const& pts) +{ + int code = validate_sum(src, fht); + if (code == TS::OK) + code = validate_point(fht, pts); + if (code == TS::OK) + code = validate_line(fht, src, pts); + return code; +} + +int TestFHT::run_func(Mat const&src, + Mat& fht) +{ + int code = TS::OK; + FastHoughTransform(src, fht, src.depth()); + return code; +} + +static Size random_size(int const max_size_log, + int const elem_size) +{ + RNG& rng = TS::ptr()->get_rng(); + return randomSize(rng, std::max(1, + max_size_log - cvRound(log(double(elem_size))))); +} + +static const int FHT_MAX_SIZE_LOG = 9; + +template +int TestFHT::run_test(int depth, + int channels, + int pts_count) +{ + int code = TS::OK; + + Size size = random_size(FHT_MAX_SIZE_LOG, + CV_ELEM_SIZE(CV_MAKE_TYPE(depth, channels))); + Mat src = Mat::zeros(size, CV_MAKETYPE(depth, channels)); + + vector pts; + code = put_random_points(src, pts_count, pts); + if (code != TS::OK) + return code; + + Mat fht; + code = run_func(src, fht); + if (code != TS::OK) + return code; + + code = validate_test_results(fht, src, pts); + return code; +} + +void TestFHT::run_n_tests(int depth, + int channels, + int pts_count, + int n) +{ + try + { + int code = TS::OK; + + switch (depth) + { + case CV_8U: + code = run_n_tests_t(depth, channels, pts_count, n); + break; + case CV_8S: + code = run_n_tests_t(depth, channels, pts_count, n); + break; + case CV_16U: + code = run_n_tests_t(depth, channels, pts_count, n); + break; + case CV_16S: + code = run_n_tests_t(depth, channels, pts_count, n); + break; + case CV_32S: + code = run_n_tests_t(depth, channels, pts_count, n); + break; + case CV_32F: + code = run_n_tests_t(depth, channels, pts_count, n); + break; + case CV_64F: + code = run_n_tests_t(depth, channels, pts_count, n); + break; + default: + code = TS::FAIL_BAD_ARG_CHECK; + ts->printf(TS::LOG, "Unknown depth %d\n", depth); + break; + } + if (code != TS::OK) + throw TS::FailureCode(code); + } + catch (const TS::FailureCode& fc) + { + std::string errorStr = TS::str_from_code(fc); + ts->printf(TS::LOG, + "General failure:\n\t%s (%d)\n", errorStr.c_str(), fc); + + ts->set_failed_test_info(fc); + } + catch(...) + { + ts->printf(TS::LOG, "Unknown failure\n"); + ts->set_failed_test_info(TS::FAIL_EXCEPTION); + } +} + +template +int TestFHT::run_n_tests_t(int depth, + int channels, + int pts_count, + int n) +{ + int code = TS::OK; + for (int iTest = 0; iTest < n; ++iTest) + { + code = run_test(depth, channels, pts_count); + if (code != TS::OK) + { + ts->printf(TS::LOG, "Test %d failed with code %d\n", iTest, code); + break; + } + } + return code; +} + +//----------------------TEST_P-------------------------------------------------- +typedef std::tr1::tuple Depth_Channels_PtsC_nPerTest; +typedef TestWithParam FastHoughTransformTest; + +TEST_P(FastHoughTransformTest, accuracy) +{ + int const depth = get<0>(GetParam()); + int const channels = get<1>(GetParam()); + int const pts_count = get<2>(GetParam()); + int const n_per_test = get<3>(GetParam()); + + TestFHT testFht; + testFht.run_n_tests(depth, channels, pts_count, n_per_test); +} + +#define FHT_ALL_DEPTHS CV_8U, CV_8S, CV_16U, CV_16S, CV_32S, CV_32F, CV_64F +#define FHT_ALL_CHANNELS 1, 2, 3, 4 + +INSTANTIATE_TEST_CASE_P(FullSet, FastHoughTransformTest, + Combine(Values(FHT_ALL_DEPTHS), + Values(FHT_ALL_CHANNELS), + Values(1, 2), + Values(5))); + +#undef FHT_ALL_DEPTHS +#undef FHT_ALL_CHANNELS + +} // namespace cvtest