Commit 7fab7a98ad48d48af18a08cb808c9fe06b795121

Authored by Tianshu Cheng
1 parent abaf5630

a neat version of mPb code

@@ -48,3 +48,5 @@ target_link_libraries(bsds500 @@ -48,3 +48,5 @@ target_link_libraries(bsds500
48 #copy an image test case 48 #copy an image test case
49 configure_file(data/101085.bmp 101085.bmp COPYONLY) 49 configure_file(data/101085.bmp 101085.bmp COPYONLY)
50 configure_file(data/101087.bmp 101087.bmp COPYONLY) 50 configure_file(data/101087.bmp 101087.bmp COPYONLY)
  51 +configure_file(data/slice00.bmp slice00.bmp COPYONLY)
  52 +configure_file(data/slice00_500_500.bmp slice00_500_500.bmp COPYONLY)
1 #include <stim/image/image.h> 1 #include <stim/image/image.h>
2 #include <cmath> 2 #include <cmath>
3 -//#include <conio.h>  
4 #include <stim/visualization/colormap.h> 3 #include <stim/visualization/colormap.h>
5 #include <stim/image/image_contour_detection.h> 4 #include <stim/image/image_contour_detection.h>
  5 +#include <sstream>
6 6
7 void array_multiply(float* lhs, float rhs, unsigned int N); 7 void array_multiply(float* lhs, float rhs, unsigned int N);
8 void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N); 8 void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
9 9
10 -/// This function evaluates the theta-dependent gradient image given an input image 10 +/// This function evaluates the theta-dependent multicale Pb given an multi-channel image
11 11
12 -/// @param lab is the 3-channel image in the LAB color space 12 +/// @param img is the multi-channel image
13 /// @param theta is the angle used for computing the gradient 13 /// @param theta is the angle used for computing the gradient
  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> func_mPb_theta(stim::image<float> img, float theta, int* r, 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> mPb_theta(w, h, 1); // allocate space for theta-dependent multicale Pb
  26 + unsigned size = mPb_theta.size(); // get the size of theta-dependent multicale Pb
  27 + memset ( mPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of theta-dependent multicale Pb to 0
  28 +
  29 +
  30 + unsigned int N = w * h; // get the number of pixels
  31 + int sigma_n = 3; // set the number of standard deviations used to define the sigma
  32 +
  33 + std::ostringstream ss; // (optional) set the stream to designate the test result file
  34 +
  35 + stim::image<float> temp; // set the temporary image to store the addtion result
  36 +
  37 + for (int i = 0; i < c; i++){
  38 + for (int j = 0; j < s; j++){
  39 +
  40 + //ss << "data_output/mPb_theta_slice"<< i*s + j << ".bmp"; // set the name for test result file
  41 + //std::string sss = ss.str();
  42 +
  43 + // get the gaussian gradient by convolving each image slice with the mask
  44 + temp = gaussian_derivative_filter_odd(img.channel(i), r[i*s + j], sigma_n, theta);
  45 +
  46 + // (optional) output the test result of each slice
  47 + //stim::cpu2image(temp.data(), sss, w, h, stim::cmBrewer);
  48 +
  49 + // multiply each gaussian gradient with its weight
  50 + array_multiply(temp.data(), alpha[i*s + j], N);
  51 +
  52 + // add up all the weighted gaussian gradients
  53 + array_add(mPb_theta.data(), temp.data(), mPb_theta.data(), N);
  54 +
  55 + //ss.str(""); //(optional) clear the space for stream
  56 +
  57 + }
  58 + }
14 59
15 -stim::image<float> func_mPb_theta(stim::image<float> lab, float theta, unsigned int w, unsigned int h){  
16 -  
17 - //allocate space for the gradient image  
18 - stim::image<float> mPb_theta(w, h, 1);  
19 -  
20 - //allocate space for each individual channel  
21 - stim::image<float> pic_light, pic_colora, pic_colorb;  
22 - pic_light = lab.channel(0);  
23 - pic_colora = lab.channel(1);  
24 - pic_colorb = lab.channel(2);  
25 -  
26 - unsigned int N = w * h; //calculate the number of pixels in the image  
27 - float sigma = 2; //set the sigma value to \sigma = 2  
28 - unsigned int sigma_n = 3; //set the number of standard deviations used to define the kernel size  
29 - unsigned r1 = 3; //disk radii  
30 - unsigned r2 = 5;  
31 - unsigned r3 = 10;  
32 - unsigned r4 = 20;  
33 - float alpha[9] = {1,1,1,1,1,1,1,1,1}; //weighting for each channel  
34 -  
35 -  
36 - stim::image<float> l1,l2,l3,a1,a2,a3,b1,b2,b3;  
37 -  
38 - l1 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r3 * 2, theta, w, h);  
39 - stim::cpu2image(l1.data(), "data_output/testl_tex5.bmp", w, h, stim::cmBrewer);  
40 -  
41 - exit(0);  
42 - /*l2 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r2 * 2, theta, w, h);  
43 - stim::cpu2image(l2.data(), "data_output/l2_tex.bmp", w, h, stim::cmBrewer);  
44 - l3 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r3 * 2, theta, w, h);  
45 - stim::cpu2image(l3.data(), "data_output/l3_tex.bmp", w, h, stim::cmBrewer);  
46 - a1 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r2 * 2, theta, w, h);  
47 - stim::cpu2image(a1.data(), "data_output/a1_tex.bmp", w, h, stim::cmBrewer);  
48 - a2 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r3 * 2, theta, w, h);  
49 - stim::cpu2image(a2.data(), "data_output/a2_tex.bmp", w, h, stim::cmBrewer);  
50 - a3 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r4 * 2, theta, w, h);  
51 - stim::cpu2image(a3.data(), "data_output/a3_tex.bmp", w, h, stim::cmBrewer);  
52 - b1 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r2 * 2, theta, w, h);  
53 - stim::cpu2image(b1.data(), "data_output/b1_tex.bmp", w, h, stim::cmBrewer);  
54 - b2 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r3 * 2, theta, w, h);  
55 - stim::cpu2image(b2.data(), "data_output/b2_tex.bmp", w, h, stim::cmBrewer);  
56 - b3 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r4 * 2, theta, w, h);  
57 - stim::cpu2image(b3.data(), "data_output/b3_tex.bmp", w, h, stim::cmBrewer);*/  
58 -  
59 - /*for (unsigned i = 0; i<N; i++){  
60 -  
61 - mPb_theta.data()[i] = l1.data()[i] * alpha[0] +  
62 - l2.data()[i] * alpha[1] +  
63 - l3.data()[i] * alpha[2] +  
64 - a1.data()[i] * alpha[3] +  
65 - a2.data()[i] * alpha[4] +  
66 - a3.data()[i] * alpha[5] +  
67 - b1.data()[i] * alpha[6] +  
68 - b2.data()[i] * alpha[7] +  
69 - b3.data()[i] * alpha[8] ;  
70 -  
71 - }*/  
72 -  
73 -  
74 - array_multiply(l1.data(), alpha[0], N);  
75 - //stim::cpu2image(l1.data(), "data_output/array_add_l1.bmp", w, h, stim::cmBrewer);  
76 - array_multiply(l2.data(), alpha[1], N);  
77 - //stim::cpu2image(l2.data(), "data_output/array_add_l2.bmp", w, h, stim::cmBrewer);  
78 - array_multiply(l3.data(), alpha[2], N);  
79 - array_multiply(a1.data(), alpha[3], N);  
80 - array_multiply(a2.data(), alpha[4], N);  
81 - array_multiply(a3.data(), alpha[5], N);  
82 - array_multiply(b1.data(), alpha[6], N);  
83 - array_multiply(b2.data(), alpha[7], N);  
84 - array_multiply(b3.data(), alpha[8], N);  
85 -  
86 - array_add(l1.data(), l2.data(), mPb_theta.data(), N);  
87 - //stim::cpu2image(sum, "data_output/array_add_sum.bmp", w, h, stim::cmBrewer);  
88 - array_add(mPb_theta.data(), l3.data(), mPb_theta.data(), N);  
89 - array_add(mPb_theta.data(), a1.data(), mPb_theta.data(), N);  
90 - array_add(mPb_theta.data(), a2.data(), mPb_theta.data(), N);  
91 - array_add(mPb_theta.data(), a3.data(), mPb_theta.data(), N);  
92 - array_add(mPb_theta.data(), b1.data(), mPb_theta.data(), N);  
93 - array_add(mPb_theta.data(), b2.data(), mPb_theta.data(), N);  
94 - array_add(mPb_theta.data(), b3.data(), mPb_theta.data(), N);  
95 -  
96 - //stim::cpu2image(mPb_theta.data(), "data_output/mPb_theta0_1.bmp", w, h, stim::cmBrewer);  
97 -  
98 -  
99 - //getch();  
100 60
101 return mPb_theta; 61 return mPb_theta;
102 } 62 }
1 #include <stim/image/image.h> 1 #include <stim/image/image.h>
2 #include <cmath> 2 #include <cmath>
3 -//#include <conio.h>  
4 #include <stim/visualization/colormap.h> 3 #include <stim/visualization/colormap.h>
5 #include <stim/image/image_contour_detection.h> 4 #include <stim/image/image_contour_detection.h>
6 #include <sstream> 5 #include <sstream>
7 6
8 -stim::image<float> func_mPb(stim::image<float> lab, unsigned int theta_n, unsigned int w, unsigned int h){ 7 +/// This function evaluates the multicale Pb given an multi-channel image
9 8
10 - std::clock_t start;  
11 - start = std::clock(); 9 +/// @param img is the multi-channel image
  10 +/// @param theta_n is the number of angles used for computing the gradient
  11 +/// @param r is an array of radii for different scaled discs(filters)
  12 +/// @param alpha is an array of weights for different scaled discs(filters)
  13 +/// @param s is the number of scales
  14 +stim::image<float> func_mPb(stim::image<float> img, unsigned int theta_n, int* r, float* alpha, int s){
12 15
13 - //---------------pavel's suggesiton------------------------------------  
14 - std::ostringstream ss;  
15 - unsigned int N = w * h;  
16 - stim::image<float> mPb_theta(w,h), mPb(w,h);  
17 - unsigned size = mPb_theta.size();  
18 - memset ( mPb.data(), 0, size * sizeof(float)); 16 + std::clock_t start; // (optional) set the timer to calculate the total time
  17 + start = std::clock(); // (optional) set timer start point
19 18
20 - float* ptr;  
21 - ptr = (float*) malloc(size * sizeof(float) * theta_n); 19 +
  20 + std::ostringstream ss; // (optional) set the stream to designate the test result file
  21 +
  22 + unsigned int w = img.width(); // get the width of picture
  23 + unsigned int h = img.height(); // get the height of picture
  24 + unsigned int N = w * h; // get the number of pixels
  25 +
  26 + stim::image<float> mPb_theta(w,h); // allocate space for theta-dependent multicale Pb (mPb_theta)
  27 + stim::image<float> mPb(w,h); // allocate space for multicale Pb (mPb)
  28 +
  29 + unsigned size = mPb.size(); // get the size of mPb
  30 + memset ( mPb.data(), 0, size * sizeof(float)); // initialize all the pixels of mPb to 0
  31 +
  32 + float* ptr; // set a pointer
  33 + ptr = (float*) malloc(size * sizeof(float) * theta_n); // this pointer points to a continuous space allocated to store all the mPb_theta
22 34
23 for (unsigned int n = 0; n < theta_n; n++){ 35 for (unsigned int n = 0; n < theta_n; n++){
24 36
25 - ss << "data_output/mPb_theta"<< n << "_conv2.bmp";  
26 - float theta = 180 * ((float)n/theta_n); 37 + ss << "data_output/mPb_theta"<< n << "_0911.bmp"; // (optional) set the name for test result file
  38 + std::string sss = ss.str(); // (optional)
  39 + float theta = 180 * ((float)n/theta_n); // calculate the even-splited angle for each mPb_theta
27 40
28 - mPb_theta = func_mPb_theta(lab, theta, w, h);  
29 - //mPb_theta.load("101087.bmp");  
30 - float* ptr_n = &ptr[ n * w * h * 1 ];  
31 - mPb_theta.channel(0).data_noninterleaved(ptr_n); 41 + mPb_theta = func_mPb_theta(img, theta, r, alpha, s); // calculate the mPb_theta
32 42
33 - double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;  
34 - std::cout<<"mPb_theta_"<< theta <<" complished time:"<< duration1 <<"s"<<'\n'; 43 + float* ptr_n = &ptr[ n * w * h * 1 ]; // set a pointer which points to the space for each mPb_theta
  44 + mPb_theta.data_noninterleaved(ptr_n); // set this pointer to point to the each mPb_theta
35 45
  46 + double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the time for generating each mPb_theta
  47 + std::cout<<"mPb_theta_"<< theta <<" complished time:"<< duration1 <<"s"<<'\n'; // (optional) show this time
36 48
37 - unsigned long idx = n * w * h * 1; //index for the nth slice  
38 49
39 - std::string sss = ss.str();  
40 - //stim::cpu2image(&ptr[idx], sss, w, h, stim::cmBrewer);  
41 - 50 + unsigned long idx = n * w * h * 1; //index for the nth mPb_theta
  51 +
  52 +
  53 + stim::cpu2image(mPb_theta.data(), sss, w, h, stim::cmBrewer); // (optional) output the nth mPb_theta
42 54
43 for(unsigned long i = 0; i < N; i++){ 55 for(unsigned long i = 0; i < N; i++){
44 56
45 - float pixel = ptr[i+idx]; //get the ith pixel in nth slice 57 + float pixel = ptr[i+idx]; //get the ith pixel in nth mPb_theta
46 58
47 - if(pixel > mPb.data()[i]){ 59 + if(pixel > mPb.data()[i]){ //get the maximum value among all mPb_theta for ith pixel
48 mPb.data()[i] = pixel; 60 mPb.data()[i] = pixel;
49 } 61 }
50 62
@@ -53,83 +65,15 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign @@ -53,83 +65,15 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; lab, unsigned int theta_n, unsign
53 } 65 }
54 66
55 67
  68 + ss.str(""); //(optional) clear the space for stream
56 69
57 - ss.str("");  
58 } 70 }
59 71
60 - //stim::cpu2image(mPb.data(), "data_output/mPb_conv2.bmp", w, h, stim::cmBrewer); 72 + stim::cpu2image(mPb.data(), "data_output/mPb_500_0911_neat.bmp", w, h, stim::cmBrewer); // output the mPb
61 73
62 - double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;  
63 - std::cout<<"total time:"<< duration2 <<"s"<<'\n';  
64 -  
65 - //getch(); 74 + double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
  75 + std::cout<<"total time:"<< duration2 <<"s"<<'\n'; // (optional) show the total time
66 76
67 return mPb; 77 return mPb;
68 78
69 - //---------------my first method------------------------------------  
70 - /*  
71 - std::clock_t start;  
72 - start = std::clock();  
73 -  
74 - stim::image<float> mPb_stack(w,h,theta_n), mPb(w,h), mPb_theta(w,h), A, B, temp;  
75 - float* ptr[8];  
76 -  
77 - for (unsigned int n = 0; n < theta_n; n++){  
78 -  
79 - //int* x = new int(5);  
80 - //int* y = x;  
81 - //*y = 1;  
82 -  
83 - float theta = 180 * ((float)n/theta_n);  
84 - mPb_theta = func_mPb_theta(lab, theta, w, h);  
85 - mPb_stack.getslice(n) = mPb_theta;  
86 - float* ptr[n] = mPb_stack.getslice(n).data();  
87 -  
88 - double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;  
89 - std::cout<<"mPb_theta, theta = "<< theta <<" time:"<< duration1 <<"s"<<'\n';  
90 -  
91 -  
92 - for(unsigned long i = 0; i < N; i++){  
93 -  
94 - *(ptr[n]+i) = mPb_theta.data()[i];  
95 -  
96 -  
97 - //float a = mPb_theta.data()[i];  
98 - //float* B = ptr[n]+i;  
99 - //A.data()[i] = mPb_theta.data()[i];  
100 - //float* C = ptr[0]+1;  
101 - //*C = 1;  
102 -  
103 - //  
104 - }  
105 - stim::cpu2image(ptr[0], "data_output/mPb_theta.bmp", w, h, stim::cmBrewer);  
106 - }  
107 -  
108 - for (unsigned long i = 0; i < N; i++){  
109 -  
110 - mPb.data()[i] = 0;  
111 - for (unsigned int n = 0; n < theta_n; n++){  
112 -  
113 - float* ptr2 = ptr[i]+n;  
114 - float temp = *ptr2;  
115 -  
116 - if(temp > mPb.data()[i]){  
117 - mPb.data()[i] = temp;  
118 - }  
119 - else{  
120 - }  
121 - }  
122 - }  
123 -  
124 - stim::cpu2image(mPb.data(), "data_output/cmap_mPb.bmp", w, h, stim::cmBrewer);  
125 -  
126 - double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;  
127 - std::cout<<"total time:"<< duration2 <<"s"<<'\n';  
128 -  
129 - getch();  
130 -  
131 - return mPb; */  
132 -  
133 -  
134 -  
135 } 79 }
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>  
5 4
6 #define PI 3.1415926 5 #define PI 3.1415926
7 6
@@ -9,80 +8,56 @@ void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned in @@ -9,80 +8,56 @@ void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned in
9 void array_abs(float* img, unsigned int N); 8 void array_abs(float* img, unsigned int N);
10 void array_multiply(float* lhs, float rhs, unsigned int N); 9 void array_multiply(float* lhs, float rhs, unsigned int N);
11 10
12 -// winsize = 2 * r, side of mask = winsize + 1  
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){ 11 +/// This function evaluates the gaussian derivative gradient given an one-channel image
14 12
15 - stim::image<float> mask_x(winsize+1, winsize+1), mask_y(winsize+1, winsize+1), mask_theta(winsize+1, winsize+1), mask_delta(winsize+1, winsize+1, 1), derivative_x, derivative_y, derivative_theta(w, h);  
16 - //float* ptr = mask_x.data(); 13 +/// @param img is the one-channel image
  14 +/// @param r is an array of radii for different scaled discs(filters)
  15 +// @param sigma_n is the number of standard deviations used to define the sigma
  16 +/// @param theta is angle used for computing the gradient
17 17
  18 +stim::image<float> gaussian_derivative_filter_odd(stim::image<float> image, int r, unsigned int sigma_n, float theta){
18 19
19 - //DEBUG calculate a Dirac delta function kernel  
20 - memset ( mask_delta.data(), 0, mask_delta.size() * sizeof(float));  
21 - mask_delta.data()[winsize*(winsize+2)/2] = 1;  
22 - stim::cpu2image(mask_delta.data(), "data_output/mask_test.bmp", winsize+1, winsize+1, stim::cmBrewer); 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 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter)
  24 + float sigma = float(r)/float(sigma_n); // calculate the sigma used in gaussian function
23 25
  26 + stim::image<float> mask_x(winsize, winsize); // allocate space for x-axis-oriented filter
  27 + stim::image<float> mask_y(winsize, winsize); // allocate space for y-axis-oriented filter
  28 + stim::image<float> mask_theta(winsize, winsize);// allocate space for theta-oriented filter
  29 + stim::image<float> derivative_theta(w, h); // allocate space for theta-oriented gradient
24 30
25 - // set parameters  
26 - unsigned N = w * h;  
27 - float theta_r = (theta * PI)/180; 31 + float theta_r = (theta * PI)/180; //change angle unit from degree to rad
28 32
29 - float step = (2*sigma*sigma_n)/winsize; 33 + for (int j = 0; j < winsize; j++){
  34 + for (int i = 0; i< winsize; i++){
30 35
31 - for (unsigned j = 0; j <= winsize; j++){  
32 - for (unsigned i = 0; i<= winsize; i++){  
33 -  
34 - float x = (-1)*sigma*sigma_n + i * step; //range of x  
35 - float y = (-1)*sigma*sigma_n + j * step; //range of y 36 + int x = i - r; //range of x
  37 + int y = j - r; //range of y
36 38
37 // create the x-oriented gaussian derivative filter mask_x 39 // create the x-oriented gaussian derivative filter mask_x
38 - mask_x.data()[j*(winsize+1) + i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); 40 + mask_x.data()[j*winsize + i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
39 // create the y-oriented gaussian derivative filter mask_y 41 // create the y-oriented gaussian derivative filter mask_y
40 - mask_y.data()[j*(winsize+1) + i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); 42 + mask_y.data()[j*winsize + i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
41 // create the mask_theta 43 // create the mask_theta
42 - mask_theta.data()[j*(winsize+1) + i] = cos(theta_r) * mask_x.data()[j*(winsize+1) + i] + sin(theta_r) * mask_y.data()[j*(winsize+1) + i] ; 44 + mask_theta.data()[j*winsize + i] = cos(theta_r) * mask_x.data()[j*winsize + i] + sin(theta_r) * mask_y.data()[j*winsize + i] ;
43 45
44 } 46 }
45 } 47 }
46 48
47 - //stim::cpu2image(mask_x.data(), "data_output/cmapgray_mask_x.bmp", winsize+1, winsize+1, stim::cmBrewer);  
48 -  
49 - stim::cpu2image(image.data(), "data_output/image.bmp", w, h, stim::cmBrewer);  
50 -  
51 -  
52 - stim::cpu2image(mask_theta.data(), "data_output/mask.bmp", winsize+1, winsize+1, stim::cmBrewer);  
53 -  
54 - // 2D convolution  
55 - //derivative_theta = image.convolve2(mask_theta);  
56 - //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta1.bmp", w, h, stim::cmBrewer);  
57 - conv2(image.data(), mask_delta.data(), derivative_theta.data(), w, h, winsize+1);  
58 - //conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize+1);  
59 - stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta_tex1.bmp", w, h, stim::cmBrewer);  
60 -  
61 - //array_abs(derivative_theta.data(), N);  
62 -  
63 - /*for (unsigned k = 0; k < w * h; k++){  
64 -  
65 - derivative_theta.data()[k] = abs(derivative_theta.data()[k]);  
66 -  
67 - }*/  
68 -  
69 - //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta2_abs.bmp", w, h, stim::cmBrewer); 49 + //stim::cpu2image(mask_theta.data(), "data_output/mask_0911_2.bmp", winsize, winsize, stim::cmBrewer); // (optional) show the mask result
70 50
71 - /*float max = derivative_theta.max(); 51 + // do the 2D convolution with image and mask
  52 + conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize);
72 53
73 - array_multiply(derivative_theta.data(), 1/max, N);*/  
74 -  
75 - /*(  
76 - for (unsigned k = 0; k < w * h; k++){  
77 -  
78 - derivative_theta.data()[k] = derivative_theta.data()[k]/max; 54 + array_abs(derivative_theta.data(), N); // get the absolute value for each pixel (why slower than the "for loop" method sometimes?)
79 55
80 - })*/ 56 + float max = derivative_theta.max(); // get the maximum of gradient used for normalization
  57 + array_multiply(derivative_theta.data(), 1/max, N); // normalize the gradient
81 58
82 - //float max2 = derivative_theta.max();  
83 59
84 - //stim::cpu2image(derivative_theta.data(), "data_output/cmap_colorb_gradient_theta90_r5.bmp", w, h, stim::cmBrewer);  
85 - //derivative_x.save("data_output/gradient_x.bmp"); 60 + //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta_0911.bmp", w, h, stim::cmBrewer); // (optional) show the gradient result
86 61
87 return derivative_theta; 62 return derivative_theta;
88 63
image_contour_detection.h deleted
1 -#include <stim/image/image.h>  
2 -//#include <cmath>  
3 -//#include <stim/visualization/colormap.h>  
4 -  
5 -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);  
6 \ No newline at end of file 0 \ No newline at end of file
@@ -3,53 +3,26 @@ @@ -3,53 +3,26 @@
3 #include <stim/visualization/colormap.h> 3 #include <stim/visualization/colormap.h>
4 #include <stim/image/image_contour_detection.h> 4 #include <stim/image/image_contour_detection.h>
5 #include <iostream> 5 #include <iostream>
6 - 6 +/// calculate the mPb given a multi-channel image
7 7
8 int main() 8 int main()
9 { 9 {
10 - stim::image<float> rgb,gaussgradient; //generate an image object  
11 -  
12 - //unsigned int a = 5%5;  
13 - //unsigned int b = 5/5;  
14 -  
15 - rgb.load("101087.bmp"); //load the input image  
16 - unsigned int w = rgb.width(); //get the image size  
17 - unsigned int h = rgb.height();  
18 - unsigned int s = rgb.size();  
19 - //unsigned a = sizeof(float);  
20 -  
21 - stim::image<float> lab; //create an image object for a single-channel (grayscale) image  
22 - lab = rgb.srgb2lab(); //create the single-channel image  
23 -  
24 - /*  
25 - stim::image<float> pic_light, pic_colora, pic_colorb;  
26 - pic_light = lab.channel(0);  
27 - pic_light.save("pic_light.bmp");  
28 -  
29 - pic_colora = lab.channel(1);  
30 - pic_colorb = lab.channel(2);  
31 -  
32 - float sigma = 2;  
33 - unsigned int sigma_n = 3;  
34 - unsigned int r = 5;  
35 - unsigned int winsize = r * 2; //window size = winsize + 1  
36 - float theta = 90;  
37 -  
38 - gaussgradient = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, winsize, theta, w, h);  
39 - gaussgradient.save("data_output/pic_gray_gradient.bmp");  
40 - */  
41 -  
42 - //float theta = 0;  
43 - unsigned int theta_n = 8;  
44 -  
45 - //stim::image<float> mPb_stack(w,h,theta_n);  
46 -  
47 - //stim::image<float> mPb_theta;  
48 - //mPb_theta = func_mPb_theta(lab, theta, w, h);  
49 - //mPb_theta.save("data_output/pic_gray_gradient.bmp");  
50 -  
51 - stim::image<float> mPb;  
52 - mPb = func_mPb(lab, theta_n, w, h); 10 + stim::image<float> img; // generate an image object
  11 +
  12 + img.load("slice00_500_500.bmp"); // load the input image
  13 + img = img.channel(0); // get the first channel of black-and-white image
  14 +
  15 + unsigned int w = img.width(); // get the width of picture
  16 + unsigned int h = img.height(); // get the height of picture
  17 + int c = img.channels(); // get the number if channels of picture
  18 + int s = 3; // set the number of scales
  19 +
  20 + int r[3] = {3,5,10}; // set an array of radii for different scaled discs(filters)
  21 + float alpha[3] = {1,1,1}; // set an array of weights for different scaled discs(filters)
  22 + unsigned int theta_n = 8; // set the number of angles used for computing the gradient
  23 +
  24 + stim::image<float> mPb; // allocate the space for mPb
  25 + mPb = func_mPb(img, theta_n, r, alpha, s); // calculate the mPb
53 26
54 return 0; 27 return 0;
55 28