Parallel implementation for Radon transform

Author: sturkmen72

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

@@ -10,24 +10,43 @@
 
 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.
-*
-*/
+ * @brief Computes the Radon transform of a given 2D image.
+ *
+ * The Radon transform is often used in image processing, particularly in applications
+ * like computed tomography, to analyze the structure of an image by projecting it along
+ * different angles. This function calculates the Radon Transform over a specified range
+ * of angles.
+ *
+ * The output type will vary depending on the input type:
+ * - 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 size of the output matrix depends on whether cropping is applied:
+ * - Without cropping, the output size will be `num_of_integral x src_diagonal_length`.
+ * - With cropping (circular), the output size will be `num_of_integral x min_edge`,
+ *   where `min_edge` is the smaller dimension of the cropped square.
+ *
+ * See https://engineering.purdue.edu/~malcolm/pct/CTI_Ch03.pdf for more details.
+ *
+ * @param src The input image on which the Radon transform is to be applied.
+ *            Must be a 2D single-channel array (e.g., grayscale image).
+ * @param dst The output array that will hold the result of the Radon transform.
+ *            The type of the output will depend on the input image type.
+ * @param theta The angle increment in degrees for the projection (resolution of the transform).
+ *              Default is 1 degree.
+ * @param start_angle The starting angle for the Radon transform in degrees.
+ *                    Default is 0 degrees.
+ * @param end_angle The ending angle for the Radon transform in degrees.
+ *                  Default is 180 degrees. The difference between end_angle and start_angle must
+ *                  be positive when multiplied by theta.
+ * @param crop A flag indicating whether to crop the input image to a square or circular shape
+ *             before the transformation. If enabled, the image is first cropped to a square
+ *             (smallest dimension) and then transformed into a circle.
+ * @param norm A flag indicating whether to normalize the output image to the range [0, 255] after
+ *             computation and convert the type to `CV_8U`. If normalization is not enabled,
+ *             the output will retain its original data range.
+ *
+ */
 CV_EXPORTS_W void RadonTransform(InputArray src,
                                       OutputArray dst,
                                       double theta = 1,

modules/ximgproc/src/radon_transform.cpp

@@ -5,77 +5,79 @@
 #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);
+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 col_num = cvRound((end_angle - start_angle) / theta);
+    int row_num, out_mat_type;
+    Point center;
+    Mat srcMat, masked_src;
 
-        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;
+    transpose(src, srcMat);
 
-        if (_srcMat.type() == CV_32FC1 || _srcMat.type() == CV_64FC1) {
-            _out_mat_type = CV_64FC1;
-        }
-        else {
-            _out_mat_type = CV_32SC1;
-        }
+    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)));
-        }
+    if (crop) {
+        // Crop the source into square
+        row_num = min(srcMat.rows, srcMat.cols);
+        Rect crop_ROI(
+            srcMat.cols / 2 - row_num / 2,
+            srcMat.rows / 2 - row_num / 2,
+            row_num, row_num);
+        srcMat = srcMat(crop_ROI);
 
-        double _t;
-        Mat _rotated_src;
-        Mat _radon(_row_num, _col_num, _out_mat_type);
+        // 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)));
+    }
 
-        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);
-        }
+    dst.create(row_num, col_num, out_mat_type);
+    Mat radon = dst.getMat();
+
+    // Define the parallel loop as a lambda function
+    parallel_for_(Range(0, col_num), [&](const Range& range) {
+        for (int col = range.start; col < range.end; col++) {
+            // Rotate the source by t
+            double t = (start_angle + col * theta);
+            Mat r_matrix = getRotationMatrix2D(center, t, 1);
+
+            Mat rotated_src;
+            warpAffine(masked_src, rotated_src, r_matrix, masked_src.size());
 
-        if (norm) {
-            normalize(_radon, _radon, 0, 255, NORM_MINMAX, CV_8UC1);
+            Mat col_mat = radon.col(col);
+            // Make projection
+            reduce(rotated_src, col_mat, 1, REDUCE_SUM, out_mat_type);
         }
+    });
 
-        _radon.copyTo(dst);
-        return;
+    if (norm) {
+        normalize(radon, dst.getMatRef(), 0, 255, NORM_MINMAX, CV_8UC1);
     }
+}
 } }

modules/ximgproc/test/test_radon_transform.cpp

@@ -9,37 +9,36 @@ 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);
 
+    ximgproc::RadonTransform(src, radon);
     EXPECT_EQ(363, radon.rows);
     EXPECT_EQ(180, radon.cols);
 
-    cv::ximgproc::RadonTransform(src, radon, 1, 0, 180, true);
-
+    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 src_float(Size(256, 256), CV_32FC1, Scalar(0));
+    Mat src_double(Size(256, 256), CV_64FC1, Scalar(0));
     Mat radon, radon_norm;
-    cv::ximgproc::RadonTransform(src_int, radon);
-    cv::ximgproc::RadonTransform(src_int, radon_norm, 1, 0, 180, false, true);
 
+    ximgproc::RadonTransform(src_int, radon);
+    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);
+    ximgproc::RadonTransform(src_float, radon);
+    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);
+
+    ximgproc::RadonTransform(src_double, radon);
+    ximgproc::RadonTransform(src_double, radon_norm, 1, 0, 180, false, true);
     EXPECT_EQ(CV_64FC1, radon.type());
     EXPECT_EQ(CV_8U, radon_norm.type());
 }
@@ -49,31 +48,29 @@ 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);
 
+    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);
+    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);
+    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]);
 }