codebase / stimlib

Commit 85025dd1591246419a0643fd7ed7ce6803fb8366

Authored by David Mayerich
1 parent afb1883d

adapted bsds500 files to the stim library, tested cPb

Showing 14 changed files with 57 additions and 849 deletions   Show diff stats
Inline Side-by-side
stim/cuda/arraymath/array_multiply.cuh
  Show/Hide comments   View file @85025dd
... ... @@ -54,6 +54,63 @@ namespace stim{
54 54 //free allocated memory
55 55 cudaFree(gpuLHS);
56 56 }
  57 +
  58 +
  59 +// array .* array multiplication
  60 +
  61 + template<typename T>
  62 + __global__ void cuda_multiply(T* ptr1, T* ptr2, T* product, unsigned int N){
  63 +
  64 + //calculate the 1D index for this thread
  65 + int idx = blockIdx.x * blockDim.x + threadIdx.x;
  66 +
  67 + if(idx < N){
  68 + product[idx] = ptr1[idx] * ptr2[idx];
  69 + }
  70 +
  71 + }
  72 +
  73 + template<typename T>
  74 + void gpu_multiply(T* ptr1, T* ptr2, T* product, unsigned int N){
  75 +
  76 + //get the maximum number of threads per block for the CUDA device
  77 + int threads = stim::maxThreadsPerBlock();
  78 +
  79 + //calculate the number of blocks
  80 + int blocks = N / threads + 1;
  81 +
  82 + //call the kernel to do the multiplication
  83 + cuda_multiply <<< blocks, threads >>>(ptr1, ptr2, product, N);
  84 +
  85 + }
  86 +
  87 + template<typename T>
  88 + void cpu_multiply(T* ptr1, T* ptr2, T* cpu_product, unsigned int N){
  89 +
  90 + //allocate memory on the GPU for the array
  91 + T* gpu_ptr1;
  92 + T* gpu_ptr2;
  93 + T* gpu_product;
  94 + HANDLE_ERROR( cudaMalloc( &gpu_ptr1, N * sizeof(T) ) );
  95 + HANDLE_ERROR( cudaMalloc( &gpu_ptr2, N * sizeof(T) ) );
  96 + HANDLE_ERROR( cudaMalloc( &gpu_product, N * sizeof(T) ) );
  97 +
  98 + //copy the array to the GPU
  99 + HANDLE_ERROR( cudaMemcpy( gpu_ptr1, ptr1, N * sizeof(T), cudaMemcpyHostToDevice) );
  100 + HANDLE_ERROR( cudaMemcpy( gpu_ptr2, ptr2, N * sizeof(T), cudaMemcpyHostToDevice) );
  101 +
  102 + //call the GPU version of this function
  103 + gpu_multiply<T>(gpu_ptr1, gpu_ptr2 ,gpu_product, N);
  104 +
  105 + //copy the array back to the CPU
  106 + HANDLE_ERROR( cudaMemcpy( cpu_product, gpu_product, N * sizeof(T), cudaMemcpyDeviceToHost) );
  107 +
  108 + //free allocated memory
  109 + cudaFree(gpu_ptr1);
  110 + cudaFree(gpu_ptr2);
  111 + cudaFree(gpu_product);
  112 +
  113 + }
57 114  
58 115 }
59 116 }
... ...
stim/cuda/bsds500/Pb.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -//#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -
6   -#define SIGMA_N 3
7   -
8   -void array_abs(float* img, unsigned int N);
9   -void array_multiply(float* lhs, float rhs, unsigned int N);
10   -void array_cos(float* ptr1, float* cpu_out, unsigned int N);
11   -void array_sin(float* ptr1, float* cpu_out, unsigned int N);
12   -void array_atan(float* ptr1, float* cpu_out, unsigned int N);
13   -void array_divide(float* ptr1, float* ptr2,float* cpu_quotient, unsigned int N);
14   -void array_multiply(float* ptr1, float* ptr2, float* product, unsigned int N);
15   -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
16   -
17   -/// This function uses odd-symmetric gaussian derivative filter to evaluate
18   -/// the max probability of a contour on one scale, given an one-channel image
19   -
20   -/// @param img is an one-channel image
21   -/// @param r is an array of radii for different scaled discs(filters)
22   -/// @param sigma_n is the number of standard deviations used to define the sigma
23   -
24   -stim::image<float> Pb(stim::image<float> image, int sigma){
25   -
26   - unsigned int w = image.width(); // get the width of picture
27   - unsigned int h = image.height(); // get the height of picture
28   - unsigned N = w * h; // get the number of pixels of picture
29   -
30   - int r = SIGMA_N * sigma; // set the radius of filter
31   - int winsize = 2 * r + 1; // set the winsdow size of filter
32   -
33   - stim::image<float> I(w, h, 1, 2); // allocate space for return image of dG1_conv2
34   - stim::image<float> theta(w, h); // allocate space for theta matrix
35   - stim::image<float> cos(w, h); // allocate space for cos(theta)
36   - stim::image<float> sin(w, h); // allocate space for sin(theta)
37   - stim::image<float> temp(w, h); // allocate space for temp
38   - stim::image<float> Ix(w, h); // allocate space for Ix
39   - stim::image<float> Iy(w, h); // allocate space for Iy
40   - stim::image<float> Pb(w, h); // allocate space for Pb
41   -
42   - I = dG1_conv2(image, sigma); // calculate the Ix, Iy
43   - Ix = I.channel(0);
44   - array_abs(Ix.data(), N); //get |Ix|;
45   - //stim::cpu2image(Ix.data(), "data_output/Pb_Ix_0924.bmp", w, h, stim::cmBrewer);
46   - Iy = I.channel(1);
47   - array_abs(Iy.data(), N); //get |Iy|;
48   - //stim::cpu2image(Iy.data(), "data_output/Pb_Iy_0924.bmp", w, h, stim::cmBrewer);
49   -
50   - array_divide(Iy.data(), Ix.data(), temp.data(), N); //temp = Iy./Ix
51   - array_atan(temp.data(), theta.data(), N); //theta = atan(temp)
52   - array_cos(theta.data(), cos.data(), N); //cos = cos(theta)
53   - array_sin(theta.data(), sin.data(), N); //sin = sin(theta)
54   - array_multiply(Ix.data(), cos.data(), Ix.data(), N); //Ix = Ix.*cos
55   - array_multiply(Iy.data(), sin.data(), Iy.data(), N); //Iy = Iy.*sin
56   - array_add(Ix.data(), Iy.data(), Pb.data(), N); //Pb = Ix + Iy;
57   -
58   - float max = Pb.maxv(); // get the maximum of Pb used for normalization
59   - array_multiply(Pb.data(), 1/max, N); // normalize the Pb
60   -
61   - //stim::cpu2image(Pb.data(), "data_output/Pb_0924.bmp", w, h, stim::cmBrewer); show the Pb(optional)
62   -
63   - return Pb;
64   -
65   -}
66 0 \ No newline at end of file
stim/cuda/bsds500/cPb.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -#include <sstream>
6   -
7   -
8   -void array_multiply(float* lhs, float rhs, unsigned int N);
9   -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
10   -
11   -/// This function evaluates the cPb given an multi-channel image
12   -
13   -/// @param img is the multi-channel image
14   -/// @param r is an array of radii for different scaled discs(filters)
15   -/// @param alpha is is an array of weights for different scaled discs(filters)
16   -/// @param s is the number of scales
17   -
18   -stim::image<float> cPb(stim::image<float> img, int* sigma, float* alpha, int s){
19   -
20   - unsigned int w = img.width(); // get the width of picture
21   - unsigned int h = img.height(); // get the height of picture
22   - unsigned int c = img.channels(); // get the channels of picture
23   -
24   -
25   - stim::image<float> cPb(w, h, 1); // allocate space for cPb
26   - unsigned size = cPb.size(); // get the size of cPb
27   - memset ( cPb.data(), 0, size * sizeof(float)); // initialize all the pixels of cPb to 0
28   -
29   -
30   - unsigned int N = w * h; // get the number of pixels
31   -
32   - //std::ostringstream ss; // (optional) set the stream to designate the test result file
33   -
34   - stim::image<float> temp; // set the temporary image to store the addtion result
35   -
36   - for (int i = 0; i < c; i++){
37   - for (int j = 0; j < s; j++){
38   -
39   - //ss << "data_output/cPb_slice"<< i*s + j << ".bmp"; // set the name for test result file (optional)
40   - //std::string sss = ss.str();
41   -
42   - // get the gaussian gradient by convolving each image slice with the mask
43   - temp = Pb(img.channel(i), sigma[i*s + j]);
44   -
45   - // output the test result of each slice (optional)
46   - //stim::cpu2image(temp.data(), sss, w, h, stim::cmBrewer);
47   -
48   - // multiply each gaussian gradient with its weight
49   - array_multiply(temp.data(), alpha[i*s + j], N);
50   -
51   - // add up all the weighted gaussian gradients
52   - array_add(cPb.data(), temp.data(), cPb.data(), N);
53   -
54   - //ss.str(""); //(optional) clear the space for stream
55   -
56   - }
57   - }
58   -
59   - float max = cPb.maxv(); // get the maximum of cPb used for normalization
60   - array_multiply(cPb.data(), 1/max, N); // normalize the cPb
61   -
62   - // output the test result of cPb (optional)
63   - //stim::cpu2image(cPb.data(), "data_output/cPb_0916.bmp", w, h, stim::cmBrewer);
64   -
65   - return cPb;
66   -}
stim/cuda/bsds500/dG1_conv2.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -//#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -#define SIGMA_N 3
6   -
7   -/// This function generates the first-order gaussian derivative filter gx gy,
8   -/// convolves the image with gx gy,
9   -/// and returns an image class which channel(0) is Ix and channel(1) is Iy
10   -
11   -/// @param img is the one-channel image
12   -/// @param sigma is the parameter for gaussian function
13   -
14   -void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1);
15   -//void array_abs(float* img, unsigned int N);
16   -
17   -stim::image<float> dG1_conv2(stim::image<float> image, int sigma){
18   -
19   - unsigned int w = image.width(); // get the width of picture
20   - unsigned int h = image.height(); // get the height of picture
21   - unsigned N = w * h; // get the number of pixels of picture
22   -
23   - int r = SIGMA_N * sigma;
24   - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
25   -
26   - stim::image<float> I(w, h, 1, 2); // allocate space for return image class
27   - stim::image<float> Ix(w, h); // allocate space for Ix
28   - stim::image<float> Iy(w, h); // allocate space for Iy
29   - Ix = image; // initialize Ix
30   - Iy = image; // initialize Iy
31   -
32   - float* array_x1;
33   - array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gx
34   - float* array_y1;
35   - array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gx
36   - float* array_x2;
37   - array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gy
38   - float* array_y2;
39   - array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gy
40   -
41   -
42   - for (int i = 0; i < winsize; i++){
43   -
44   - int x = i - r; //range of x
45   - int y = i - r; //range of y
46   -
47   - // create the 1D x-oriented gaussian derivative filter array_x1 for gx
48   - array_x1[i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
49   - // create the 1D y-oriented gaussian derivative filter array_y1 for gx
50   - array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
51   - // create the 1D x-oriented gaussian derivative filter array_x2 for gy
52   - array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
53   - // create the 1D y-oriented gaussian derivative filter array_y2 for gy
54   - array_y2[i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
55   - }
56   -
57   - //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
58   - //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
59   - //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
60   - //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
61   -
62   - // get Ix by convolving the image with gx
63   - conv2_sep(Ix.data(), w, h, array_x1, winsize, array_y1, winsize);
64   -
65   - //stim::cpu2image(Ix.data(), "data_output/Ix_0915.bmp", w, h, stim::cmBrewer);
66   - // get Iy by convolving the image with gy
67   - conv2_sep(Iy.data(), w, h, array_x2, winsize, array_y2, winsize);
68   -
69   - //stim::cpu2image(Iy.data(), "data_output/Iy_0915.bmp", w, h, stim::cmBrewer);
70   -
71   - delete [] array_x1; //free the memory of array_x1
72   - delete [] array_y1; //free the memory of array_y1
73   - delete [] array_x2; //free the memory of array_x2
74   - delete [] array_y2; //free the memory of array_y2
75   -
76   - I.set_channel(0, Ix.data());
77   - I.set_channel(1, Iy.data());
78   -
79   - return I;
80   -
81   -}
82 0 \ No newline at end of file
stim/cuda/bsds500/dG1_theta_conv2.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -
6   -#define PI 3.1415926
7   -#define SIGMA_N 3
8   -
9   -void array_multiply(float* lhs, float rhs, unsigned int N);
10   -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
11   -void array_abs(float* img, unsigned int N);
12   -
13   -/// This function evaluates the theta-dependent odd symmetric gaussian derivative gradient of an one-channel image
14   -
15   -/// @param img is the one-channel image
16   -/// @param r is an array of radii for different scaled discs(filters)
17   -/// @param sigma_n is the number of standard deviations used to define the sigma
18   -/// @param theta is angle used for computing the gradient
19   -
20   -stim::image<float> dG1_theta_conv2(stim::image<float> image, int sigma, float theta){
21   -
22   - float theta_r = (theta * PI)/180; //change angle unit from degree to rad
23   -
24   - unsigned int w = image.width(); // get the width of picture
25   - unsigned int h = image.height(); // get the height of picture
26   -
27   - int r = SIGMA_N * sigma;
28   - unsigned N = w * h; // get the number of pixels of picture
29   -
30   - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
31   -
32   - stim::image<float> I(w, h, 1, 2); // allocate space for return image of dG1_conv2
33   - stim::image<float> Ix(w, h); // allocate space for Ix
34   - stim::image<float> Iy(w, h); // allocate space for Iy
35   - stim::image<float> dG1_theta(w, h); // allocate space for Pb
36   -
37   - I = dG1_conv2(image, sigma); // calculate the Ix, Iy
38   - Ix = I.channel(0);
39   - Iy = I.channel(1);
40   -
41   - array_multiply(Ix.data(), cos(theta_r), N); //Ix = Ix*cos(theta_r)
42   - array_multiply(Iy.data(), sin(theta_r), N); //Iy = Iy*sin(theta_r)
43   - array_add(Ix.data(), Iy.data(), dG1_theta.data(), N); //dG1_theta = Ix + Iy;
44   - array_abs(dG1_theta.data(), N);
45   -
46   - //stim::cpu2image(I.channel(0).data(), "data_output/dG1_theta_x_0919.bmp", w, h, stim::cmBrewer);
47   - //stim::cpu2image(I.channel(1).data(), "data_output/dG1_theta_y_0919.bmp", w, h, stim::cmBrewer);
48   - //stim::cpu2image(dG1_theta.data(), "data_output/dG1_theta_0919.bmp", w, h, stim::cmBrewer);
49   -
50   - return dG1_theta;
51   -
52   -}
stim/cuda/bsds500/dG2_conv2.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -//#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -#define SIGMA_N 3
6   -
7   -/// This function generates the second-order gaussian derivative filter gxx gyy,
8   -/// convolves the image with gxx gyy,
9   -/// and returns an image class which channel(0) is Ixx and channel(1) is Iyy
10   -
11   -/// @param img is the one-channel image
12   -/// @param sigma is the parameter for gaussian function
13   -
14   -void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1);
15   -//void array_abs(float* img, unsigned int N);
16   -
17   -stim::image<float> dG2_conv2(stim::image<float> image, int sigma){
18   -
19   - unsigned int w = image.width(); // get the width of picture
20   - unsigned int h = image.height(); // get the height of picture
21   - unsigned N = w * h; // get the number of pixels of picture
22   -
23   - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
24   - int r = SIGMA_N * sigma;
25   -
26   - stim::image<float> I(w, h, 1, 2); // allocate space for return image class
27   - stim::image<float> Ixx(w, h); // allocate space for Ixx
28   - stim::image<float> Iyy(w, h); // allocate space for Iyy
29   - Ixx = image; // initialize Ixx
30   - Iyy = image; // initialize Iyy
31   -
32   - float* array_x1;
33   - array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gxx
34   - float* array_y1;
35   - array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gxx
36   - float* array_x2;
37   - array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gyy
38   - float* array_y2;
39   - array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gyy
40   -
41   -
42   - for (int i = 0; i < winsize; i++){
43   -
44   - int x = i - r; //range of x
45   - int y = i - r; //range of y
46   -
47   - // create the 1D x-oriented gaussian derivative filter array_x1 for gxx
48   - array_x1[i] = (-1) * (1 - pow(x, 2)) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
49   - // create the 1D y-oriented gaussian derivative filter array_y1 for gxx
50   - array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
51   - // create the 1D x-oriented gaussian derivative filter array_x2 for gyy
52   - array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
53   - // create the 1D y-oriented gaussian derivative filter array_y2 for gyy
54   - array_y2[i] = (-1) * (1 - pow(y, 2)) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
55   - }
56   -
57   - //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
58   - //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
59   - //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
60   - //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
61   -
62   - // get Ixx by convolving the image with gxx
63   - conv2_sep(Ixx.data(), w, h, array_x1, winsize, array_y1, winsize);
64   -
65   - //stim::cpu2image(Ixx.data(), "data_output/Ixx_0915.bmp", w, h, stim::cmBrewer);
66   - // get Iyy by convolving the image with gyy
67   - conv2_sep(Iyy.data(), w, h, array_x2, winsize, array_y2, winsize);
68   -
69   - //stim::cpu2image(Iyy.data(), "data_output/Iyy_0915.bmp", w, h, stim::cmBrewer);
70   -
71   - delete [] array_x1; //free the memory of array_x1
72   - delete [] array_y1; //free the memory of array_y1
73   - delete [] array_x2; //free the memory of array_x2
74   - delete [] array_y2; //free the memory of array_y2
75   -
76   - I.set_channel(0, Ixx.data());
77   - I.set_channel(1, Iyy.data());
78   -
79   - return I;
80   -
81   -}
82 0 \ No newline at end of file
stim/cuda/bsds500/dG2_d2x_theta_conv2.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -
6   -#define SIGMA_N 3
7   -
8   -/// This function evaluates the theta-dependent even-symmetric gaussian derivative gradient of an one-channel image
9   -
10   -/// @param img is the one-channel image
11   -/// @param r is an array of radii for different scaled discs(filters)
12   -/// @param sigma_n is the number of standard deviations used to define the sigma
13   -/// @param theta is angle used for computing the gradient
14   -
15   -void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M);
16   -void array_abs(float* img, unsigned int N);
17   -
18   -stim::image<float> dG2_d2x_theta_conv2(stim::image<float> image, int sigma, float theta){
19   -
20   - unsigned int w = image.width(); // get the width of picture
21   - unsigned int h = image.height(); // get the height of picture
22   - unsigned N = w * h; // get the number of pixels of picture
23   -
24   - int r = SIGMA_N * sigma; // set the radius of filter
25   - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
26   -
27   - stim::image<float> I(w, h, 1, 2); // allocate space for return image class
28   - stim::image<float> dG2_d2x_theta(w, h); // allocate space for dG2_d2x_theta
29   - stim::image<float> mask_x(winsize, winsize); // allocate space for x-axis-oriented filter
30   - stim::image<float> mask_r(winsize, winsize); // allocate space for theta-oriented filter
31   -
32   - for (int j = 0; j < winsize; j++){
33   - for (int i = 0; i< winsize; i++){
34   -
35   - int x = i - r; //range of x
36   - int y = j - r; //range of y
37   -
38   - // create the x-oriented gaussian derivative filter mask_x
39   - mask_x.data()[j*winsize + i] = (-1) * (1 - pow(x, 2)) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
40   -
41   - }
42   - }
43   -
44   - mask_r = mask_x.rotate(theta, r, r);
45   - //mask_r = mask_x.rotate(45, r, r);
46   - //stim::cpu2image(mask_r.data(), "data_output/mask_r_0919.bmp", winsize, winsize, stim::cmBrewer);
47   -
48   - // do the 2D convolution with image and mask
49   - conv2(image.data(), mask_r.data(), dG2_d2x_theta.data(), w, h, winsize);
50   - array_abs(dG2_d2x_theta.data(), N);
51   -
52   - //stim::cpu2image(dG2_d2x_theta.data(), "data_output/dG2_d2x_theta_0919.bmp", w, h, stim::cmGrayscale);
53   -
54   - return dG2_d2x_theta;
55   -}
56 0 \ No newline at end of file
stim/cuda/bsds500/dG_d2x_theta_conv2.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -
6   -/// This function evaluates the theta-dependent even-symmetric gaussian derivative gradient of an one-channel image
7   -
8   -/// @param img is the one-channel image
9   -/// @param r is an array of radii for different scaled discs(filters)
10   -/// @param sigma_n is the number of standard deviations used to define the sigma
11   -/// @param theta is angle used for computing the gradient
12   -
13   -void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M);
14   -void array_abs(float* img, unsigned int N);
15   -
16   -stim::image<float> Gd_even(stim::image<float> image, int r, unsigned int sigma_n, float theta){
17   -
18   - unsigned int w = image.width(); // get the width of picture
19   - unsigned int h = image.height(); // get the height of picture
20   - unsigned N = w * h; // get the number of pixels of picture
21   - int winsize = 2 * r + 1; // set the winsdow size of disc(filter)
22   - float sigma = float(r)/float(sigma_n); // calculate the sigma used in gaussian function
23   -
24   - stim::image<float> I(w, h, 1, 2); // allocate space for return image class
25   - stim::image<float> Gd_even_theta(w, h); // allocate space for Gd_even_theta
26   - stim::image<float> mask_x(winsize, winsize); // allocate space for x-axis-oriented filter
27   - stim::image<float> mask_r(winsize, winsize); // allocate space for theta-oriented filter
28   -
29   - for (int j = 0; j < winsize; j++){
30   - for (int i = 0; i< winsize; i++){
31   -
32   - int x = i - r; //range of x
33   - int y = j - r; //range of y
34   -
35   - // create the x-oriented gaussian derivative filter mask_x
36   - mask_x.data()[j*winsize + i] = (-1) * (1 - pow(x, 2)) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
37   -
38   - }
39   - }
40   -
41   - mask_r = mask_x.rotate(theta, r, r);
42   - //mask_r = mask_x.rotate(45, r, r);
43   - //stim::cpu2image(mask_r.data(), "data_output/mask_r_0919.bmp", winsize, winsize, stim::cmBrewer);
44   -
45   - // do the 2D convolution with image and mask
46   - conv2(image.data(), mask_r.data(), Gd_even_theta.data(), w, h, winsize);
47   - array_abs(Gd_even_theta.data(), N);
48   -
49   - //stim::cpu2image(Gd_even_theta.data(), "data_output/Gd_even_theta_0919.bmp", w, h, stim::cmGrayscale);
50   -
51   - return Gd_even_theta;
52   -}
53 0 \ No newline at end of file
stim/cuda/bsds500/kmeans.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -//#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -#include <opencv2/opencv.hpp>
6   -#include <iostream>
7   -
8   -/// This function use cvkmeans to cluster given textons
9   -
10   -/// @param testons is a multi-channel image
11   -/// @param K is the number of clusters
12   -
13   -stim::image<float> kmeans(stim::image<float> textons, unsigned int K){
14   -
15   - unsigned int w = textons.width(); // get the width of picture
16   - unsigned int h = textons.height(); // get the height of picture
17   - unsigned int feature_n = textons.channels(); // get the spectrum of picture
18   - unsigned int N = w * h; // get the number of pixels
19   -
20   - float* sample1 = (float*) malloc(sizeof(float) * N * feature_n); //allocate the space for textons
21   -
22   - //reallocate a multi-channel texton image to a single-channel image
23   - for(unsigned int c = 0; c < feature_n; c++){
24   -
25   - stim::image<float> temp;
26   - temp = textons.channel(c);
27   -
28   - for(unsigned int j = 0; j < N; j++){
29   -
30   - sample1[c + j * feature_n] = temp.data()[j];
31   - }
32   - }
33   -
34   -
35   - cv::Mat sample2(N, feature_n, CV_32F, sample1); //copy image to cv::mat
36   -
37   - //(optional) show the test result
38   - //imshow("sample2", sample2);
39   -
40   -
41   - cv::TermCriteria criteria(CV_TERMCRIT_EPS+CV_TERMCRIT_ITER, 10, 0.1); // set stop-criteria for kmeans iteration
42   - cv::Mat labels(N, 1, CV_8U, cvScalarAll(0)); // allocate space for kmeans output
43   - cv::Mat centers; // allocate space for kmeans output
44   -
45   - unsigned int test_times = 2; // set the number of times of trying kmeans, it will return the best result
46   -
47   - cv::kmeans(sample2, K, labels, criteria, test_times, cv::KMEANS_PP_CENTERS, centers); // kmeans clustering
48   -
49   - //(optional) show the test result
50   - //imwrite( "data_output/labels_1D.bmp", labels);
51   -
52   - stim::image<float> texture(w, h, 1, 1); // allocate space for texture
53   -
54   - for(unsigned int i = 0; i < N; i++){ // reshape the labels from iD array to image
55   -
56   - texture.data()[i] = labels.at<int>(i);
57   -
58   - }
59   -
60   - //texture.save("data_output/kmeans_test0924_2.bmp");
61   -
62   - //(optional) show the test result
63   - //stim::cpu2image(texture.data(), "data_output/kmeans_test.bmp", w, h, stim::cmBrewer);
64   -
65   - return texture;
66   -
67   -}
68 0 \ No newline at end of file
stim/cuda/bsds500/laplacian_conv2.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -
6   -#define PI 3.1415926
7   -#define SIGMA_N 3
8   -
9   -void array_multiply(float* lhs, float rhs, unsigned int N);
10   -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
11   -void array_abs(float* img, unsigned int N);
12   -
13   -/// This function evaluates the center-surround(Laplacian of Gaussian) gaussian derivative gradient of an one-channel image
14   -
15   -/// @param img is the one-channel image
16   -/// @param r is an array of radii for different scaled discs(filters)
17   -/// @param sigma_n is the number of standard deviations used to define the sigma
18   -
19   -stim::image<float> laplacian_conv2(stim::image<float> image, int sigma){
20   -
21   - unsigned int w = image.width(); // get the width of picture
22   - unsigned int h = image.height(); // get the height of picture
23   - unsigned N = w * h; // get the number of pixels of picture
24   -
25   - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
26   -
27   - stim::image<float> I(w, h, 1, 2); // allocate space for return image of dG2_conv2
28   - stim::image<float> Ixx(w, h); // allocate space for Ixx
29   - stim::image<float> Iyy(w, h); // allocate space for Iyy
30   - stim::image<float> laplacian(w, h); // allocate space for Pb
31   -
32   - I = dG2_conv2(image, sigma); // calculate the Ixx, Iyy
33   - Ixx = I.channel(0);
34   - Iyy = I.channel(1);
35   -
36   - array_add(Ixx.data(), Iyy.data(), laplacian.data(), N); //laplacian = Ixx + Iyy;
37   - array_abs(laplacian.data(), N);
38   -
39   - //stim::cpu2image(laplacian.data(), "data_output/laplacian_0919.bmp", w, h, stim::cmBrewer);
40   -
41   - return laplacian;
42   -
43   -}
stim/cuda/bsds500/main.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -#include <iostream>
6   -#include <stim/visualization/colormap.h>
7   -#include <sstream>
8   -/// calculate the cPb, tPb and mPb given a multi-channel image
9   -
10   -void array_multiply(float* lhs, float rhs, unsigned int N);
11   -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
12   -
13   -int main()
14   -{
15   - stim::image<float> image; // generate an image object
16   -
17   - //std::ostringstream ss1,ss2; // (optional) set the stream to designate the test result file
18   -
19   - //for(unsigned int i = 27; i < 31; i++){
20   -
21   - //ss1 << "3d_sample/b0"<< i << ".bmp"; // (optional) set the name for test result file
22   - //std::string sss1 = ss1.str(); // (optional)
23   -
24   - //image.load("bone1001_3.bmp");
25   - //image.load(sss1);
26   - //image.load("101085.bmp"); // load the input image
27   - image.load("slice00_500_500.bmp"); // load the input image
28   - image = image.channel(0); // get the first channel of black-and-white image
29   -
30   -
31   - unsigned int w = image.width(); // get the width of picture
32   - unsigned int h = image.height(); // get the height of picture
33   - unsigned int N = w * h; // get the number of pixels
34   - int c = image.channels(); // get the number if channels of picture
35   - int s = 3; // set the number of scales
36   - int sigma[3] = {1,2,3}; // set an array of sigmas for different scaled gaussian filters for cPb, the length of array is c*s
37   - int r[3] = {5,10,20}; // set an array of radii for different scaled discs(filters) for tPb, the length of array is c*s
38   - float alpha[3] = {1,1,1}; // set an array of weights for different scaled discs(filters)
39   - unsigned int theta_n = 8; // set the number of angles used in filter orientation
40   - unsigned int bin_n = 16; // set the number of bins used in chi-distance
41   - unsigned int K = 16; // set the number of cludters, K should be multiple of bin_n
42   -
43   - stim::image<float> img_cPb(w, h); // allocate the space for cPb
44   - stim::image<float> img_tPb(w, h); // allocate the space for tPb
45   - stim::image<float> img_mPb(w, h); // allocate the space for tPb
46   -
47   - std::cout<<"imagesize: "<< w <<"*"<< h <<'\n';
48   -
49   - //*******************cPb********************
50   - std::clock_t start1; // (optional) set the timer to calculate the total time
51   - start1 = std::clock(); // (optional) set timer start point
52   -
53   - std::cout<<"begin cPb"<<'\n';
54   -
55   - img_cPb = cPb(image, sigma, alpha, s);
56   -
57   - // show the cPb (optional)
58   - stim::cpu2image(img_cPb.data(), "data_output/img_cPb_1008.bmp", w, h, stim::cmBrewer);
59   - //ss2 << "3d_sample/cPb0"<< i << ".bmp"; // (optional) set the name for test result file
60   - //std::string sss2 = ss2.str();
61   -
62   - //stim::cpu2image(img_cPb.data(), "data_output/bone_cPb3.1.bmp", w, h, stim::cmBrewer);
63   -
64   - double duration1 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
65   - std::cout<<"cPb time: "<< duration1 <<"s"<<'\n'; // (optional) show the total time
66   -
67   - //ss1.str(""); //(optional) clear the space for stream
68   - //ss2.str(""); //(optional) clear the space for stream
69   -
70   -
71   - //*******************tPb********************
72   -
73   - std::cout<<"begin tPb"<<'\n';
74   -
75   -
76   - std::clock_t start2; // (optional) set the timer to calculate the total time
77   - start2 = std::clock(); // (optional) set timer start point
78   -
79   - img_tPb = tPb(image, r, alpha, theta_n, bin_n, s, K);
80   -
81   - // show the tPb (optional)
82   - // stim::cpu2image(img_tPb.data(), "data_output/bone_tPb_3.1.bmp", w, h, stim::cmBrewer);
83   - stim::cpu2image(img_tPb.data(), "data_output/img_tPb_1008.bmp", w, h, stim::cmBrewer);
84   -
85   - double duration2 = ( std::clock() - start2 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
86   - std::cout<<"tPb time: "<< duration2 <<"s"<<'\n'; // (optional) show the total time
87   -
88   -
89   - //*******************************************
90   -
91   - double duration3 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
92   - std::cout<<"total running time: "<< duration3 <<"s"<<'\n'; // (optional) show the total time
93   - //}
94   - //******************mPb**********************
95   - // set parameters for linear combination
96   -
97   - float a = 1;
98   - float b = 0.5;
99   -
100   - // create mPb by linearly combined cPb and tPb
101   - array_multiply(img_cPb.data(), a, N);
102   - array_multiply(img_tPb.data(), b, N);
103   - array_add(img_cPb.data(), img_tPb.data(), img_mPb.data(), N);
104   -
105   - //ss2 << "3d_sample/mPb0"<< i << ".bmp"; // (optional) set the name for test result file
106   - //std::string sss2 = ss2.str();
107   - // show the mPb (optional)
108   - stim::cpu2image(img_mPb.data(), "data_output/img_mPb_1008.bmp", w, h, stim::cmBrewer);
109   -
110   - //stim::cpu2image(img_mPb.data(), sss2, w, h, stim::cmBrewer);
111   -
112   - //ss1.str(""); //(optional) clear the space for stream
113   - //ss2.str(""); //(optional) clear the space for stream
114   -
115   - //}
116   - return 0;
117   -
118   -}
stim/cuda/bsds500/tPb.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -#include <sstream>
6   -
7   -
8   -void array_multiply(float* lhs, float rhs, unsigned int N);
9   -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
10   -void chi_grad(float* img, float* cpu_copy, unsigned int w, unsigned int h, int r, unsigned int bin_n, unsigned int bin_size, float theta);
11   -
12   -/// This function evaluates the tPb given a grayscale image
13   -
14   -/// @param img is the multi-channel image
15   -/// @param theta_n is the number of angles used for computing oriented chi-gradient
16   -/// @param r is an array of radii for different scaled discs(filters)
17   -/// @param alpha is is an array of weights for different scaled discs(filters)
18   -/// @param s is the number of scales
19   -/// @param K is the number of clusters
20   -
21   -stim::image<float> tPb(stim::image<float> img, int* r, float* alpha, unsigned int theta_n, unsigned int bin_n, int s, unsigned K){
22   -
23   - unsigned int w = img.width(); // get the width of picture
24   - unsigned int h = img.height(); // get the height of picture
25   - unsigned int N = w * h; // get the number of pixels
26   -
27   - stim::image<float> img_textons(w, h, 1, theta_n*2+1); // allocate space for img_textons
28   - stim::image<float> img_texture(w, h, 1, 1); // allocate space for img_texture
29   - stim::image<float> tPb_theta(w, h, 1, 1); // allocate space for tPb_theta
30   - stim::image<float> tPb(w, h, 1, 1); // allocate space for tPb
31   - unsigned size = tPb_theta.size(); // get the size of tPb_theta
32   - memset (tPb.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb to 0
33   - stim::image<float> temp(w, h, 1, 1); // set the temporary image to store the addtion result
34   -
35   - std::ostringstream ss; // (optional) set the stream to designate the test result file
36   -
37   -
38   - img_textons = textons(img, theta_n);
39   -
40   - img_texture = kmeans(img_textons, K); // changing kmeans result into float type is required
41   -
42   - stim::cpu2image(img_texture.data(), "data_output/texture.bmp", w, h, stim::cmBrewer);
43   -
44   -
45   - unsigned int max1 = img_texture.maxv(); // get the maximum of Pb used for normalization
46   - unsigned int bin_size = (max1 + 1)/bin_n; // (whether"+1" or not depends on kmeans result)
47   -
48   - for (int i = 0; i < theta_n; i++){
49   -
50   - float theta = 180 * ((float)i/theta_n); // calculate the even-splited angle for each tPb_theta
51   -
52   - memset (tPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb_theta to 0
53   -
54   - //ss << "data_output/0922tPb_theta"<< theta << ".bmp"; // set the name for test result file (optional)
55   - //std::string sss = ss.str();
56   -
57   - for (int j = 0; j < s; j++){
58   -
59   - // get the chi-gradient by convolving each image slice with the mask
60   - chi_grad(img_texture.data(), temp.data(), w, h, r[j], bin_n, bin_size, theta);
61   -
62   - float max2 = temp.maxv(); // get the maximum of tPb_theta used for normalization
63   - array_multiply(temp.data(), 1/max2, N); // normalize the tPb_theta
64   -
65   - //output the test result of each slice (optional)
66   - //stim::cpu2image(temp.data(), "data_output/tPb_slice0924_2.bmp", w, h, stim::cmBrewer);
67   -
68   - // multiply each chi-gradient with its weight
69   - array_multiply(temp.data(), alpha[j], N);
70   -
71   - // add up all the weighted chi-gradients
72   - array_add(tPb_theta.data(), temp.data(), tPb_theta.data(), N);
73   -
74   -
75   - }
76   -
77   - //ss.str(""); //(optional) clear the space for stream
78   -
79   - for(unsigned long ti = 0; ti < N; ti++){
80   -
81   - if(tPb_theta.data()[ti] > tPb.data()[ti]){ //get the maximum value among all tPb_theta for ith pixel
82   - tPb.data()[ti] = tPb_theta.data()[ti];
83   - }
84   -
85   - else{
86   - }
87   - }
88   - }
89   -
90   - float max3 = tPb.maxv(); // get the maximum of tPb used for normalization
91   - array_multiply(tPb.data(), 1/max3, N); // normalize the tPb
92   -
93   - //output the test result of tPb (optional)
94   - //stim::cpu2image(tPb.data(), "data_output/tPb_0922.bmp", w, h, stim::cmBrewer);
95   -
96   - return tPb;
97   -}
stim/cuda/bsds500/textons.cpp deleted
View file @afb1883
1   -#include <stim/image/image.h>
2   -//#include <cmath>
3   -#include <stim/visualization/colormap.h>
4   -#include <stim/image/image_contour_detection.h>
5   -#include <sstream>
6   -
7   -/// This function convolve the grayscale image with a set of oriented Gaussian
8   -/// derivative filters, and return a texton image with (theta_n*2+1) channels
9   -
10   -/// @param image is an one-channel grayscale image
11   -/// @param theta_n is the number of angles used for computing the gradient
12   -
13   -stim::image<float> textons(stim::image<float> image, unsigned int theta_n){
14   -
15   - unsigned int w = image.width(); // get the width of picture
16   - unsigned int h = image.height(); // get the height of picture
17   - unsigned N = w * h; // get the number of pixels of picture
18   -
19   - stim::image<float> textons(w, h, 1, theta_n*2+1); // allocate space for textons
20   - stim::image<float> temp(w, h); // allocate space for temp
21   -
22   - int sigma = 1; // set sigma for odd-symmetric, even-symmetric and center-surround filter filter
23   -
24   - //std::ostringstream ss1, ss2; // (optional) set the stream to designate the test result file
25   -
26   - for (unsigned int i = 0; i < theta_n; i++){
27   -
28   - //ss1 << "data_output/textons_channel_"<< i << ".bmp"; // set the name for test result file (optional)
29   - //std::string sss1 = ss1.str();
30   - //ss2 << "data_output/textons_channel_"<< i+theta_n << ".bmp"; // set the name for test result file (optional)
31   - //std::string sss2 = ss2.str();
32   -
33   - float theta = 180 * ((float)i/theta_n); // calculate the even-splited angle for each oriented filter
34   -
35   - temp = dG1_theta_conv2(image, sigma, theta); // return dG1_theta_conv2 to temp
36   - //stim::cpu2image(temp.data(), sss1, w, h, stim::cmBrewer);
37   - textons.set_channel(i, temp.data()); // copy temp to ith channel of textons
38   -
39   - temp = dG2_d2x_theta_conv2(image, sigma, theta); // return dG2_d2x_theta_conv2 to temp
40   - //stim::cpu2image(temp.data(), sss2, w, h, stim::cmBrewer);
41   - textons.set_channel(i + theta_n, temp.data()); // copy temp to (i+theta_n)th channel of textons
42   -
43   - //ss1.str(""); //(optional) clear the space for stream
44   - //ss2.str(""); //(optional) clear the space for stream
45   -
46   - }
47   -
48   - temp = laplacian_conv2(image, sigma); // return laplacian_conv2 to temp
49   - //stim::cpu2image(temp.data(), "data_output/textons_channel_16.bmp", w, h, stim::cmBrewer);
50   - textons.set_channel(theta_n*2, temp.data()); // copy temp to (theta_n*2)th channel of textons
51   -
52   - return textons;
53   -
54   -}
55   -
56   -
57 0 \ No newline at end of file
stim/image/image_contour_detection.h deleted
View file @afb1883
1   -
2   -//stim::image<float> gaussian_derivative_filter_odd(stim::image<float> image, int r, unsigned int sigma_n, float theta);
3   -//stim::image<float> func_mPb_theta(stim::image<float> img, float theta, int* r, float* alpha, int s);
4   -//stim::image<float> func_mPb(stim::image<float> img, unsigned int theta_n, int* r, float* alpha, int s);
5   -
6   -stim::image<float> dG1_conv2(stim::image<float> image, int sigma);
7   -stim::image<float> dG2_conv2(stim::image<float> image, int sigma);
8   -stim::image<float> dG1_theta_conv2(stim::image<float> image, int sigma, float theta);
9   -stim::image<float> dG2_d2x_theta_conv2(stim::image<float> image, int sigma, float theta);
10   -stim::image<float> laplacian_conv2(stim::image<float> image, int sigma);
11   -
12   -stim::image<float> textons(stim::image<float> image, unsigned int theta_n);
13   -stim::image<float> kmeans(stim::image<float> textons, unsigned int K);
14   -stim::image<float> Pb(stim::image<float> image, int sigma);
15   -stim::image<float> cPb(stim::image<float> img, int* sigma, float* alpha, int s);
16   -stim::image<float> tPb(stim::image<float> img, int* r, float* alpha, unsigned int theta_n, unsigned int bin_n, int s, unsigned int K);
17 0 \ No newline at end of file