Add Radon transform function to ximgproc

modules/ximgproc/include/opencv2/ximgproc.hpp

@@ -60,6 +60,7 @@
 #include "ximgproc/run_length_morphology.hpp"
 #include "ximgproc/edgepreserving_filter.hpp"
 #include "ximgproc/color_match.hpp"
+#include "ximgproc/radon_transform.hpp"
 
 
 /** @defgroup ximgproc Extended Image Processing

modules/ximgproc/include/opencv2/ximgproc/radon_transform.hpp

@@ -0,0 +1,40 @@
+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#ifndef __OPENCV_RADON_TRANSFORM_HPP__
+#define __OPENCV_RADON_TRANSFORM_HPP__
+
+#include "opencv2/core.hpp"
+#include "opencv2/imgproc.hpp"
+
+namespace cv { namespace ximgproc {
+/**
+* @brief   Calculate Radon Transform of an image.
+* @param   src         The source (input) image.
+* @param   dst         The destination image, result of transformation.
+* @param   theta       Angle resolution of the transform in degrees.
+* @param   start_angle Start angle of the transform in degrees.
+* @param   end_angle   End angle of the transform in degrees.
+* @param   crop        Crop the source image into a circle.
+* @param   norm        Normalize the output Mat to grayscale and convert type to CV_8U
+*
+* This function calculates the Radon Transform of a given image in any range.
+* See https://engineering.purdue.edu/~malcolm/pct/CTI_Ch03.pdf for detail.
+* If the input type is CV_8U, the output will be CV_32S.
+* If the input type is CV_32F or CV_64F, the output will be CV_64F
+* The output size will be num_of_integral x src_diagonal_length.
+* If crop is selected, the input image will be crop into square then circle,
+* and output size will be num_of_integral x min_edge.
+*
+*/
+CV_EXPORTS_W void RadonTransform(InputArray src,
+                                      OutputArray dst,
+                                      double theta = 1,
+                                      double start_angle = 0,
+                                      double end_angle = 180,
+                                      bool crop = false,
+                                      bool norm = false);
+} }
+
+#endif

modules/ximgproc/perf/perf_radon_transform.cpp

@@ -0,0 +1,35 @@
+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#include "perf_precomp.hpp"
+
+namespace opencv_test { namespace {
+
+typedef tuple<Size, MatType> RadonTransformPerfTestParam;
+typedef perf::TestBaseWithParam<RadonTransformPerfTestParam> RadonTransformPerfTest;
+
+PERF_TEST_P(RadonTransformPerfTest, perf,
+    testing::Combine(
+        testing::Values(TYPICAL_MAT_SIZES),
+        testing::Values(CV_8UC1, CV_32FC1, CV_64FC1)
+    )
+)
+{
+    Size srcSize = get<0>(GetParam());
+    int  srcType = get<1>(GetParam());
+
+    Mat src(srcSize, srcType);
+    Mat radon;
+
+    declare.in(src, WARMUP_RNG);
+
+    TEST_CYCLE()
+    {
+        RadonTransform(src, radon);
+    }
+
+    SANITY_CHECK_NOTHING();
+}
+
+} }

modules/ximgproc/samples/radon_transform_demo.cpp

@@ -0,0 +1,18 @@
+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#include <opencv2/highgui.hpp>
+#include <opencv2/ximgproc/radon_transform.hpp>
+
+using namespace cv;
+
+int main() {
+    Mat src = imread("peilin_plane.png", IMREAD_GRAYSCALE);
+    Mat radon;
+    ximgproc::RadonTransform(src, radon, 1, 0, 180, false, true);
+    imshow("src image", src);
+    imshow("Radon transform", radon);
+    waitKey();
+    return 0;
+}

modules/ximgproc/samples/radon_transform_demo.py

@@ -0,0 +1,13 @@
+# This file is part of OpenCV project.
+# It is subject to the license terms in the LICENSE file found in the top-level directory
+# of this distribution and at http://opencv.org/license.html.
+
+import numpy as np
+import cv2 as cv
+
+if __name__ == "__main__":
+    src = cv.imread("peilin_plane.png", cv.IMREAD_GRAYSCALE)
+    radon = cv.ximgproc.RadonTransform(src)
+    cv.imshow("src image", src)
+    cv.imshow("Radon transform", radon)
+    cv.waitKey()

modules/ximgproc/src/radon_transform.cpp

@@ -0,0 +1,81 @@
+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#include "precomp.hpp"
+
+namespace cv {namespace ximgproc {
+    void RadonTransform(InputArray src,
+                             OutputArray dst,
+                             double theta,
+                             double start_angle,
+                             double end_angle,
+                             bool crop,
+                             bool norm)
+    {
+        CV_Assert(src.dims() == 2);
+        CV_Assert(src.channels() == 1);
+        CV_Assert((end_angle - start_angle) * theta > 0);
+
+        Mat _srcMat = src.getMat();
+
+        int _row_num, _col_num, _out_mat_type;
+        _col_num = cvRound((end_angle - start_angle) / theta);
+        transpose(_srcMat, _srcMat);
+        Mat _masked_src;
+        cv::Point _center;
+
+        if (_srcMat.type() == CV_32FC1 || _srcMat.type() == CV_64FC1) {
+            _out_mat_type = CV_64FC1;
+        }
+        else {
+            _out_mat_type = CV_32SC1;
+        }
+
+        if (crop) {
+            // crop the source into square
+            _row_num = min(_srcMat.rows, _srcMat.cols);
+            cv::Rect _crop_ROI(
+                _srcMat.cols / 2 - _row_num / 2,
+                _srcMat.rows / 2 - _row_num / 2,
+                _row_num, _row_num);
+            _srcMat = _srcMat(_crop_ROI);
+            // crop the source into circle
+            Mat _mask(_srcMat.size(), CV_8UC1, Scalar(0));
+            _center = Point(_srcMat.cols / 2, _srcMat.rows / 2);
+            circle(_mask, _center, _srcMat.cols / 2, Scalar(255), FILLED);
+            _srcMat.copyTo(_masked_src, _mask);
+        }
+        else {
+            // avoid cropping corner when rotating
+            _row_num = cvCeil(sqrt(_srcMat.rows * _srcMat.rows + _srcMat.cols * _srcMat.cols));
+            _masked_src = Mat(Size(_row_num, _row_num), _srcMat.type(), Scalar(0));
+            _center = Point(_masked_src.cols / 2, _masked_src.rows / 2);
+            _srcMat.copyTo(_masked_src(Rect(
+                (_row_num - _srcMat.cols) / 2,
+                (_row_num - _srcMat.rows) / 2,
+                _srcMat.cols, _srcMat.rows)));
+        }
+
+        double _t;
+        Mat _rotated_src;
+        Mat _radon(_row_num, _col_num, _out_mat_type);
+
+        for (int _col = 0; _col < _col_num; _col++) {
+            // rotate the source by _t
+            _t = (start_angle + _col * theta);
+            cv::Mat _r_matrix = cv::getRotationMatrix2D(_center, _t, 1);
+            cv::warpAffine(_masked_src, _rotated_src, _r_matrix, _masked_src.size());
+            Mat _col_mat = _radon.col(_col);
+            // make projection
+            cv::reduce(_rotated_src, _col_mat, 1, REDUCE_SUM, _out_mat_type);
+        }
+
+        if (norm) {
+            normalize(_radon, _radon, 0, 255, NORM_MINMAX, CV_8UC1);
+        }
+
+        _radon.copyTo(dst);
+        return;
+    }
+} }

modules/ximgproc/test/test_radon_transform.cpp

@@ -0,0 +1,81 @@
+// This file is part of OpenCV project.
+// It is subject to the license terms in the LICENSE file found in the top-level directory
+// of this distribution and at http://opencv.org/license.html.
+
+#include "test_precomp.hpp"
+
+namespace opencv_test {namespace {
+
+TEST(RadonTransformTest, output_size)
+{
+    Mat src(Size(256, 256), CV_8U, Scalar(0));
+    circle(src, Point(128, 128), 64, Scalar(255), FILLED);
+    Mat radon;
+    cv::ximgproc::RadonTransform(src, radon);
+
+    EXPECT_EQ(363, radon.rows);
+    EXPECT_EQ(180, radon.cols);
+
+    cv::ximgproc::RadonTransform(src, radon, 1, 0, 180, true);
+
+    EXPECT_EQ(256, radon.rows);
+    EXPECT_EQ(180, radon.cols);
+}
+
+TEST(RadonTransformTest, output_type)
+{
+    Mat src_int(Size(256, 256), CV_8U, Scalar(0));
+    circle(src_int, Point(128, 128), 64, Scalar(255), FILLED);
+    Mat radon, radon_norm;
+    cv::ximgproc::RadonTransform(src_int, radon);
+    cv::ximgproc::RadonTransform(src_int, radon_norm, 1, 0, 180, false, true);
+
+    EXPECT_EQ(CV_32SC1, radon.type());
+    EXPECT_EQ(CV_8U, radon_norm.type());
+
+    Mat src_float(Size(256, 256), CV_32FC1, Scalar(0));
+    Mat src_double(Size(256, 256), CV_32FC1, Scalar(0));
+    cv::ximgproc::RadonTransform(src_float, radon);
+    cv::ximgproc::RadonTransform(src_float, radon_norm, 1, 0, 180, false, true);
+    EXPECT_EQ(CV_64FC1, radon.type());
+    EXPECT_EQ(CV_8U, radon_norm.type());
+    cv::ximgproc::RadonTransform(src_double, radon);
+    EXPECT_EQ(CV_64FC1, radon.type());
+    EXPECT_EQ(CV_8U, radon_norm.type());
+}
+
+TEST(RadonTransformTest, accuracy_by_pixel)
+{
+    Mat src(Size(256, 256), CV_8U, Scalar(0));
+    circle(src, Point(128, 128), 64, Scalar(255), FILLED);
+    Mat radon;
+    cv::ximgproc::RadonTransform(src, radon);
+
+    ASSERT_EQ(CV_32SC1, radon.type());
+
+    EXPECT_EQ(0, radon.at<int>(0, 0));
+
+    EXPECT_LT(18000, radon.at<int>(128, 128));
+    EXPECT_GT(19000, radon.at<int>(128, 128));
+}
+
+TEST(RadonTransformTest, accuracy_uchar)
+{
+    Mat src(Size(10, 10), CV_8UC1, Scalar(1));
+    cv::Mat radon;
+    ximgproc::RadonTransform(src, radon, 45, 0, 180, false, false);
+
+    EXPECT_EQ(100, sum(radon.col(0))[0]);
+}
+
+TEST(RadonTransformTest, accuracy_float)
+{
+    Mat src(Size(10, 10), CV_32FC1, Scalar(1.1));
+    cv::Mat radon;
+    ximgproc::RadonTransform(src, radon, 45, 0, 180, false, false);
+
+    EXPECT_LT(109, sum(radon.col(0))[0]);
+    EXPECT_GT(111, sum(radon.col(0))[0]);
+}
+
+} }