Julia

CMakeLists.txt

@@ -114,6 +114,7 @@ add_executable(test_algo ${PERF_SOURCES})
 target_link_libraries(test_algo
   algos
   "${binary_dir}/lib/libopencv_ts.a"
+  "${binary_dir}/lib/libopencv_imgcodecs.so"
   pthread
 )
 target_include_directories(test_algo PUBLIC

include/algos.hpp

@@ -24,6 +24,9 @@ void bgr2gray_opencv(const cv::Mat& src, cv::Mat& dst);
 void boxFilter_halide(uint16_t* src, uint16_t* dst, int height, int width);
 void boxFilter_opencv(const cv::Mat& src, cv::Mat& dst);
 
+void julia_ref(uint8_t* dst, int height, int width);
+void halide_julia(uint8_t* dst, int height, int width);
+
 #ifdef HAVE_OPENCV_DNN
 void convolution_nchw_halide(float* src, float* kernel, float* dst,
                              int inpChannels, int outChannels, int height, int width);

src/julia.cpp

@@ -0,0 +1,133 @@
+#include <opencv2/opencv.hpp>
+
+#ifdef __riscv
+#include <HalideBuffer.h>
+#include "bgr2gray_interleaved.h"
+#include "bgr2gray_planar.h"
+
+using namespace Halide::Runtime;
+
+#else
+
+#include <Halide.h>
+
+using namespace Halide;
+
+#endif
+
+static const float ci = -0.75;
+static const float cr = 0.0;
+
+// src :: Int -> fraction :: Float -> Float
+Expr halide_julia_norm(Expr src, Expr fraction) {
+    return (((cast<float>(src)) * 4.0f) / fraction) - 2.0f;
+}
+
+struct Complex {
+    Expr real, imag;
+
+    Complex(Tuple t)
+        : real(t[0]), imag(t[1]) {
+    }
+
+    Complex(Expr r, Expr i)
+        : real(r), imag(i) {
+    }
+
+    Complex(FuncRef t)
+        : Complex(Tuple(t)) {
+    }
+
+    operator Tuple() const {
+        return {real, imag};
+    }
+
+    Complex operator+(const Complex &other) const {
+        return {real + other.real, imag + other.imag};
+    }
+
+    Complex operator*(const Complex &other) const {
+        return {real * other.real - imag * other.imag,
+                real * other.imag + imag * other.real};
+    }
+
+    Expr magnitude_squared() const {
+        return real * real + imag * imag;
+    }
+};
+
+void halide_julia(uint8_t* dst, int height, int width) {
+    static const uint8_t step_bound = 255;
+    static const uint8_t upper_bound = 4;
+
+    Buffer output(dst, width, height);
+
+    Var x, y;
+
+    Complex c(cr, ci);
+
+    Expr x_ranged;
+    x_ranged = halide_julia_norm(x, height);
+
+    Expr y_ranged;
+    y_ranged = halide_julia_norm(y, width);
+
+    Func julia;
+    Var t;
+    julia(x, y, t) = Complex(x_ranged, y_ranged);
+
+    // loop
+    RDom index(1, step_bound);
+    Complex current = julia(x, y, index - 1);
+    julia(x, y, index) = current * current + c;
+
+    Expr esc_cond = Complex(julia(x, y, index)).magnitude_squared() < upper_bound;
+    Tuple first_escape = argmin(esc_cond);
+
+    // proj to result
+    Func result;
+    result(x, y) = cast<uint8_t>(first_escape[0]);
+
+    result.realize(output);
+}
+
+float julia_norm(int src, float fraction) {
+    return ((((float)src) * 4.0) / fraction) - 2.0;
+}
+
+uint8_t julia_step(float x_src, float y_src) {
+    static const uint8_t count_bound = 255;
+    static const uint8_t upper_bound = 4;
+
+    float x = x_src;
+    float y = y_src;
+
+    int count = 0;
+
+    for (;count < count_bound && x * x + y * y < upper_bound; count++) {
+        float new_x = x * x - y * y + cr;
+        float new_y = 2 * x * y + ci;
+
+        x = new_x;
+        y = new_y;
+    }
+
+    return count;
+}
+
+void julia_ref(uint8_t* dst, int height, int width) {
+    static const int channel_number = 1;
+
+    for (int y = 0; y < height; ++y) {
+        uint8_t* dstData = dst + y * width * channel_number;
+
+        for (int x = 0; x < width; ++x) {
+            const float ranged_x = julia_norm(x, (float)height);
+            const float ranged_y = julia_norm(y, (float)width);
+
+            const uint8_t result = julia_step(ranged_x, ranged_y);
+
+            dstData[x] = result;
+        }
+    }
+}

test/test_main.cpp

@@ -19,8 +19,27 @@ using namespace cv;
 const int width = 1920;
 const int height = 1080;
 
+static const int julia_width = 200;
+static const int julia_height = 200;
+
 CV_TEST_MAIN("")
 
+TEST(julia, ref) {
+    Mat dst(julia_height, julia_width, CV_8UC1, Scalar(0));
+
+    julia_ref(dst.ptr<uint8_t>(), dst.rows, dst.cols);
+
+    imwrite("julia_ref.png", dst);
+}
+
+TEST(julia, halide) {
+    Mat dst(julia_height, julia_width, CV_8UC1, Scalar(0));
+
+    halide_julia(dst.ptr<uint8_t>(), dst.rows, dst.cols);
+
+    imwrite("julia_halide.png", dst);
+}
+
 TEST(histogram, opencv) {
     Mat src(height, width, CV_8UC3), dst(3, 256, CV_32F), ref(3, 256, CV_32S);
     randu(src, 0, 256);