Commit 40d11588d501812254ef602402ec5143c4080da3

Authored by Tianshu Cheng
1 parent b5c39899

2D inseparable convolution

1 -#include <stim/cuda/gaussian_blur.cuh> 1 +#include <stim/cuda/arraymath.cuh>
2 2
3 -void blur(float* image, float sigma, unsigned int x, unsigned int y){ 3 +/*void blur(float* image, float sigma, unsigned int x, unsigned int y){
4 4
5 stim::cuda::cpu_gaussian_blur_2d<float>(image, sigma, x, y); 5 stim::cuda::cpu_gaussian_blur_2d<float>(image, sigma, x, y);
6 -}  
7 \ No newline at end of file 6 \ No newline at end of file
  7 +}*/
  8 +
  9 +void array_multiply(float* lhs, float rhs, unsigned int N){
  10 +
  11 + stim::cuda::cpu_multiply(lhs, rhs, N);
  12 +}
  13 +
  14 +void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N){
  15 +
  16 + stim::cuda::cpu_add(ptr1, ptr2, sum, N);
  17 +
  18 +}
  19 +
  20 +void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M){
  21 +
  22 + stim::cuda::cpu_conv2(img, mask, cpu_copy, w, h, M);
  23 +
  24 +}
  25 +
  26 +void array_abs(float* img, unsigned int N){
  27 +
  28 + stim::cuda::cpu_abs(img, N);
  29 +
  30 +}
  31 +
@@ -4,6 +4,9 @@ @@ -4,6 +4,9 @@
4 #include <stim/visualization/colormap.h> 4 #include <stim/visualization/colormap.h>
5 #include <stim/image/image_contour_detection.h> 5 #include <stim/image/image_contour_detection.h>
6 6
  7 +void array_multiply(float* lhs, float rhs, unsigned int N);
  8 +void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
  9 +
7 stim::image<float> func_mPb_theta(stim::image<float> lab, float theta, unsigned int w, unsigned int h){ 10 stim::image<float> func_mPb_theta(stim::image<float> lab, float theta, unsigned int w, unsigned int h){
8 11
9 stim::image<float> mPb_theta(w, h, 1); 12 stim::image<float> mPb_theta(w, h, 1);
@@ -26,16 +29,25 @@ stim::image&lt;float&gt; func_mPb_theta(stim::image&lt;float&gt; lab, float theta, unsigned @@ -26,16 +29,25 @@ stim::image&lt;float&gt; func_mPb_theta(stim::image&lt;float&gt; lab, float theta, unsigned
26 stim::image<float> l1,l2,l3,a1,a2,a3,b1,b2,b3; 29 stim::image<float> l1,l2,l3,a1,a2,a3,b1,b2,b3;
27 30
28 l1 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r1 * 2, theta, w, h); 31 l1 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r1 * 2, theta, w, h);
  32 + stim::cpu2image(l1.data(), "data_output/l1_tex.bmp", w, h, stim::cmBrewer);
29 l2 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r2 * 2, theta, w, h); 33 l2 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r2 * 2, theta, w, h);
  34 + stim::cpu2image(l2.data(), "data_output/l2_tex.bmp", w, h, stim::cmBrewer);
30 l3 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r3 * 2, theta, w, h); 35 l3 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r3 * 2, theta, w, h);
  36 + stim::cpu2image(l3.data(), "data_output/l3_tex.bmp", w, h, stim::cmBrewer);
31 a1 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r2 * 2, theta, w, h); 37 a1 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r2 * 2, theta, w, h);
  38 + stim::cpu2image(a1.data(), "data_output/a1_tex.bmp", w, h, stim::cmBrewer);
32 a2 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r3 * 2, theta, w, h); 39 a2 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r3 * 2, theta, w, h);
  40 + stim::cpu2image(a2.data(), "data_output/a2_tex.bmp", w, h, stim::cmBrewer);
33 a3 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r4 * 2, theta, w, h); 41 a3 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r4 * 2, theta, w, h);
  42 + stim::cpu2image(a3.data(), "data_output/a3_tex.bmp", w, h, stim::cmBrewer);
34 b1 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r2 * 2, theta, w, h); 43 b1 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r2 * 2, theta, w, h);
  44 + stim::cpu2image(b1.data(), "data_output/b1_tex.bmp", w, h, stim::cmBrewer);
35 b2 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r3 * 2, theta, w, h); 45 b2 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r3 * 2, theta, w, h);
  46 + stim::cpu2image(b2.data(), "data_output/b2_tex.bmp", w, h, stim::cmBrewer);
36 b3 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r4 * 2, theta, w, h); 47 b3 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r4 * 2, theta, w, h);
  48 + stim::cpu2image(b3.data(), "data_output/b3_tex.bmp", w, h, stim::cmBrewer);
37 49
38 - for (unsigned i = 0; i<N; i++){ 50 + /*for (unsigned i = 0; i<N; i++){
39 51
40 mPb_theta.data()[i] = l1.data()[i] * alpha[0] + 52 mPb_theta.data()[i] = l1.data()[i] * alpha[0] +
41 l2.data()[i] * alpha[1] + 53 l2.data()[i] * alpha[1] +
@@ -47,9 +59,32 @@ stim::image&lt;float&gt; func_mPb_theta(stim::image&lt;float&gt; lab, float theta, unsigned @@ -47,9 +59,32 @@ stim::image&lt;float&gt; func_mPb_theta(stim::image&lt;float&gt; lab, float theta, unsigned
47 b2.data()[i] * alpha[7] + 59 b2.data()[i] * alpha[7] +
48 b3.data()[i] * alpha[8] ; 60 b3.data()[i] * alpha[8] ;
49 61
50 - } 62 + }*/
  63 +
  64 +
  65 + array_multiply(l1.data(), alpha[0], N);
  66 + //stim::cpu2image(l1.data(), "data_output/array_add_l1.bmp", w, h, stim::cmBrewer);
  67 + array_multiply(l2.data(), alpha[1], N);
  68 + //stim::cpu2image(l2.data(), "data_output/array_add_l2.bmp", w, h, stim::cmBrewer);
  69 + array_multiply(l3.data(), alpha[2], N);
  70 + array_multiply(a1.data(), alpha[3], N);
  71 + array_multiply(a2.data(), alpha[4], N);
  72 + array_multiply(a3.data(), alpha[5], N);
  73 + array_multiply(b1.data(), alpha[6], N);
  74 + array_multiply(b2.data(), alpha[7], N);
  75 + array_multiply(b3.data(), alpha[8], N);
  76 +
  77 + array_add(l1.data(), l2.data(), mPb_theta.data(), N);
  78 + //stim::cpu2image(sum, "data_output/array_add_sum.bmp", w, h, stim::cmBrewer);
  79 + array_add(mPb_theta.data(), l3.data(), mPb_theta.data(), N);
  80 + array_add(mPb_theta.data(), a1.data(), mPb_theta.data(), N);
  81 + array_add(mPb_theta.data(), a2.data(), mPb_theta.data(), N);
  82 + array_add(mPb_theta.data(), a3.data(), mPb_theta.data(), N);
  83 + array_add(mPb_theta.data(), b1.data(), mPb_theta.data(), N);
  84 + array_add(mPb_theta.data(), b2.data(), mPb_theta.data(), N);
  85 + array_add(mPb_theta.data(), b3.data(), mPb_theta.data(), N);
51 86
52 - //stim::cpu2image(mPb_theta.data(), "data_output/cmap_mPb_theta0.bmp", w, h, stim::cmBrewer); 87 + //stim::cpu2image(mPb_theta.data(), "data_output/mPb_theta0_1.bmp", w, h, stim::cmBrewer);
53 88
54 89
55 //getch(); 90 //getch();
@@ -20,9 +20,9 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign @@ -20,9 +20,9 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign
20 float* ptr; 20 float* ptr;
21 ptr = (float*) malloc(size * sizeof(float) * theta_n); 21 ptr = (float*) malloc(size * sizeof(float) * theta_n);
22 22
23 - for (unsigned int n = 0; n < 1; n++){ 23 + for (unsigned int n = 0; n < theta_n; n++){
24 24
25 - ss << "data_output/mPb_theta"<< n << ".bmp"; 25 + ss << "data_output/mPb_theta"<< n << "_conv2.bmp";
26 float theta = 180 * ((float)n/theta_n); 26 float theta = 180 * ((float)n/theta_n);
27 27
28 mPb_theta = func_mPb_theta(lab, theta, w, h); 28 mPb_theta = func_mPb_theta(lab, theta, w, h);
@@ -37,7 +37,7 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign @@ -37,7 +37,7 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign
37 unsigned long idx = n * w * h * 1; //index for the nth slice 37 unsigned long idx = n * w * h * 1; //index for the nth slice
38 38
39 std::string sss = ss.str(); 39 std::string sss = ss.str();
40 - stim::cpu2image(&ptr[idx], sss, w, h, stim::cmBrewer); 40 + //stim::cpu2image(&ptr[idx], sss, w, h, stim::cmBrewer);
41 41
42 42
43 for(unsigned long i = 0; i < N; i++){ 43 for(unsigned long i = 0; i < N; i++){
@@ -57,7 +57,7 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign @@ -57,7 +57,7 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign
57 ss.str(""); 57 ss.str("");
58 } 58 }
59 59
60 - stim::cpu2image(mPb.data(), "data_output/mPb.bmp", w, h, stim::cmBrewer); 60 + //stim::cpu2image(mPb.data(), "data_output/mPb_conv2.bmp", w, h, stim::cmBrewer);
61 61
62 double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; 62 double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
63 std::cout<<"total time:"<< duration2 <<"s"<<'\n'; 63 std::cout<<"total time:"<< duration2 <<"s"<<'\n';
gauss_derivative_odd.cpp
1 #include <stim/image/image.h> 1 #include <stim/image/image.h>
2 #include <cmath> 2 #include <cmath>
3 #include <stim/visualization/colormap.h> 3 #include <stim/visualization/colormap.h>
4 -#include <iostream> 4 +//#include <iostream>
5 5
6 #define PI 3.1415926 6 #define PI 3.1415926
7 7
  8 +void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M);
  9 +void array_abs(float* img, unsigned int N);
  10 +void array_multiply(float* lhs, float rhs, unsigned int N);
  11 +
  12 +// winsize = 2 * r, side of mask = winsize + 1
8 stim::image<float> gaussian_derivative_filter_odd(stim::image<float> image, float sigma, unsigned int sigma_n, unsigned int winsize, float theta, unsigned int w, unsigned int h){ 13 stim::image<float> gaussian_derivative_filter_odd(stim::image<float> image, float sigma, unsigned int sigma_n, unsigned int winsize, float theta, unsigned int w, unsigned int h){
9 14
10 stim::image<float> mask_x(winsize+1, winsize+1), mask_y(winsize+1, winsize+1), mask_theta(winsize+1, winsize+1), derivative_x, derivative_y, derivative_theta(w, h); 15 stim::image<float> mask_x(winsize+1, winsize+1), mask_y(winsize+1, winsize+1), mask_theta(winsize+1, winsize+1), derivative_x, derivative_y, derivative_theta(w, h);
@@ -38,27 +43,37 @@ stim::image&lt;float&gt; gaussian_derivative_filter_odd(stim::image&lt;float&gt; image, floa @@ -38,27 +43,37 @@ stim::image&lt;float&gt; gaussian_derivative_filter_odd(stim::image&lt;float&gt; image, floa
38 43
39 //stim::cpu2image(mask_x.data(), "data_output/cmapgray_mask_x.bmp", winsize+1, winsize+1, stim::cmBrewer); 44 //stim::cpu2image(mask_x.data(), "data_output/cmapgray_mask_x.bmp", winsize+1, winsize+1, stim::cmBrewer);
40 45
41 - //stim::cpu2image(mask_y.data(), "data_output/cmapgray_mask_y.bmp", winsize+1, winsize+1, stim::cmBrewer); 46 + stim::cpu2image(image.data(), "data_output/image.bmp", w, h, stim::cmBrewer);
42 47
43 48
44 - //stim::cpu2image(mask_theta.data(), "data_output/cmapgray_mask_theta.bmp", winsize+1, winsize+1, stim::cmBrewer); 49 + stim::cpu2image(mask_theta.data(), "data_output/mask.bmp", winsize+1, winsize+1, stim::cmBrewer);
45 50
46 // 2D convolution 51 // 2D convolution
47 - derivative_theta = image.convolve2(mask_theta); 52 + //derivative_theta = image.convolve2(mask_theta);
  53 + //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta1.bmp", w, h, stim::cmBrewer);
  54 + conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize+1);
  55 + //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta_tex1.bmp", w, h, stim::cmBrewer);
48 56
49 - for (unsigned k = 0; k < w * h; k++){ 57 + //array_abs(derivative_theta.data(), N);
  58 +
  59 + /*for (unsigned k = 0; k < w * h; k++){
50 60
51 derivative_theta.data()[k] = abs(derivative_theta.data()[k]); 61 derivative_theta.data()[k] = abs(derivative_theta.data()[k]);
52 62
53 - } 63 + }*/
54 64
55 - float max = derivative_theta.max(); 65 + //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta2_abs.bmp", w, h, stim::cmBrewer);
56 66
  67 + /*float max = derivative_theta.max();
  68 +
  69 + array_multiply(derivative_theta.data(), 1/max, N);*/
  70 +
  71 + /*(
57 for (unsigned k = 0; k < w * h; k++){ 72 for (unsigned k = 0; k < w * h; k++){
58 73
59 derivative_theta.data()[k] = derivative_theta.data()[k]/max; 74 derivative_theta.data()[k] = derivative_theta.data()[k]/max;
60 75
61 - } 76 + })*/
62 77
63 //float max2 = derivative_theta.max(); 78 //float max2 = derivative_theta.max();
64 79
@@ -9,11 +9,14 @@ void main() @@ -9,11 +9,14 @@ void main()
9 { 9 {
10 stim::image<float> rgb,gaussgradient; //generate an image object 10 stim::image<float> rgb,gaussgradient; //generate an image object
11 11
  12 + //unsigned int a = 5%5;
  13 + //unsigned int b = 5/5;
  14 +
12 rgb.load("101087.bmp"); //load the input image 15 rgb.load("101087.bmp"); //load the input image
13 unsigned int w = rgb.width(); //get the image size 16 unsigned int w = rgb.width(); //get the image size
14 unsigned int h = rgb.height(); 17 unsigned int h = rgb.height();
15 unsigned int s = rgb.size(); 18 unsigned int s = rgb.size();
16 - unsigned a = sizeof(float); 19 + //unsigned a = sizeof(float);
17 20
18 stim::image<float> lab; //create an image object for a single-channel (grayscale) image 21 stim::image<float> lab; //create an image object for a single-channel (grayscale) image
19 lab = rgb.srgb2lab(); //create the single-channel image 22 lab = rgb.srgb2lab(); //create the single-channel image