segmentation / bsds500

Commit 6dcc460e79106d5240be259023155d84ea07dae2

Authored by Tianshu Cheng
1 parent 7fab7a98

cPb+tPb

Showing 17 changed files with 844 additions and 13 deletions   Show diff stats
Inline Side-by-side
CMakeLists.txt
  Show/Hide comments   View file @6dcc460
... ... @@ -7,6 +7,9 @@ project(bsds500)
7 7 #set the module directory
8 8 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}")
9 9  
  10 +#find OpenCV
  11 +find_package(OpenCV REQUIRED )
  12 +
10 13 #set up CUDA
11 14 find_package(CUDA REQUIRED)
12 15  
... ... @@ -20,6 +23,7 @@ find_package(Threads)
20 23 find_package(X11)
21 24  
22 25 include_directories(
  26 + ${OpenCV_INCLUDE_DIRS}
23 27 ${STIM_INCLUDE_DIRS}
24 28 )
25 29  
... ... @@ -41,6 +45,7 @@ cuda_add_executable(bsds500
41 45 target_link_libraries(bsds500
42 46 #${CUDA_cufft_LIBRARY}
43 47 #${CUDA_cublas_LIBRARY}
  48 + ${OpenCV_LIBS}
44 49 ${CMAKE_THREAD_LIBS_INIT}
45 50 ${X11_LIBRARIES}
46 51 )
... ... @@ -48,5 +53,6 @@ target_link_libraries(bsds500
48 53 #copy an image test case
49 54 configure_file(data/101085.bmp 101085.bmp COPYONLY)
50 55 configure_file(data/101087.bmp 101087.bmp COPYONLY)
  56 +configure_file(data/101087_square.bmp 101087_square.bmp COPYONLY)
51 57 configure_file(data/slice00.bmp slice00.bmp COPYONLY)
52 58 configure_file(data/slice00_500_500.bmp slice00_500_500.bmp COPYONLY)
... ...
Pb.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 +void array_abs(float* img, unsigned int N);
  7 +void array_multiply(float* lhs, float rhs, unsigned int N);
  8 +void array_cos(float* ptr1, float* cpu_out, unsigned int N);
  9 +void array_sin(float* ptr1, float* cpu_out, unsigned int N);
  10 +void array_atan(float* ptr1, float* cpu_out, unsigned int N);
  11 +void array_divide(float* ptr1, float* ptr2,float* cpu_quotient, unsigned int N);
  12 +void array_multiply(float* ptr1, float* ptr2, float* product, unsigned int N);
  13 +void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
  14 +
  15 +/// This function uses odd-symmetric gaussian derivative filter to evaluate
  16 +/// the max probability of a contour on one scale, given an one-channel image
  17 +
  18 +/// @param img is an one-channel image
  19 +/// @param r is an array of radii for different scaled discs(filters)
  20 +/// @param sigma_n is the number of standard deviations used to define the sigma
  21 +
  22 +stim::image<float> Pb(stim::image<float> image, int r, unsigned int sigma_n){
  23 +
  24 + unsigned int w = image.width(); // get the width of picture
  25 + unsigned int h = image.height(); // get the height of picture
  26 + unsigned N = w * h; // get the number of pixels of picture
  27 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter)
  28 +
  29 + stim::image<float> I(w, h, 1, 2); // allocate space for return image of Gd1
  30 + stim::image<float> theta(w, h); // allocate space for theta matrix
  31 + stim::image<float> cos(w, h); // allocate space for cos(theta)
  32 + stim::image<float> sin(w, h); // allocate space for sin(theta)
  33 + stim::image<float> temp(w, h); // allocate space for temp
  34 + stim::image<float> Ix(w, h); // allocate space for Ix
  35 + stim::image<float> Iy(w, h); // allocate space for Iy
  36 + stim::image<float> Pb(w, h); // allocate space for Pb
  37 +
  38 + I = Gd1(image, r, sigma_n); // calculate the Ix, Iy
  39 + Ix = I.channel(0);
  40 + array_abs(Ix.data(), N); //get |Ix|;
  41 + //stim::cpu2image(Ix.data(), "data_output/Pb_Ix_0924.bmp", w, h, stim::cmBrewer);
  42 + Iy = I.channel(1);
  43 + array_abs(Iy.data(), N); //get |Iy|;
  44 + //stim::cpu2image(Iy.data(), "data_output/Pb_Iy_0924.bmp", w, h, stim::cmBrewer);
  45 +
  46 + array_divide(Iy.data(), Ix.data(), temp.data(), N); //temp = Iy./Ix
  47 + array_atan(temp.data(), theta.data(), N); //theta = atan(temp)
  48 + array_cos(theta.data(), cos.data(), N); //cos = cos(theta)
  49 + array_sin(theta.data(), sin.data(), N); //sin = sin(theta)
  50 + array_multiply(Ix.data(), cos.data(), Ix.data(), N); //Ix = Ix.*cos
  51 + array_multiply(Iy.data(), sin.data(), Iy.data(), N); //Iy = Iy.*sin
  52 + array_add(Ix.data(), Iy.data(), Pb.data(), N); //Pb = Ix + Iy;
  53 +
  54 + float max = Pb.maxv(); // get the maximum of Pb used for normalization
  55 + array_multiply(Pb.data(), 1/max, N); // normalize the Pb
  56 +
  57 + //stim::cpu2image(Pb.data(), "data_output/Pb_0924.bmp", w, h, stim::cmBrewer); show the Pb(optional)
  58 +
  59 + return Pb;
  60 +
  61 +}
0 62 \ No newline at end of file
... ...
cPb.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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* 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> 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 + 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/cPb_slice"<< i*s + j << ".bmp"; // set the name for test result file (optional)
  41 + std::string sss = ss.str();
  42 +
  43 + // get the gaussian gradient by convolving each image slice with the mask
  44 + temp = Pb(img.channel(i), r[i*s + j], sigma_n);
  45 +
  46 + // output the test result of each slice (optional)
  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(cPb.data(), temp.data(), cPb.data(), N);
  54 +
  55 + ss.str(""); //(optional) clear the space for stream
  56 +
  57 + }
  58 + }
  59 +
  60 + float max = cPb.maxv(); // get the maximum of cPb used for normalization
  61 + array_multiply(cPb.data(), 1/max, N); // normalize the cPb
  62 +
  63 + // output the test result of cPb (optional)
  64 + //stim::cpu2image(cPb.data(), "data_output/cPb_0916.bmp", w, h, stim::cmBrewer);
  65 +
  66 + return cPb;
  67 +}
... ...
cudafunc.cu
  Show/Hide comments   View file @6dcc460
1 1 #include <stim/cuda/arraymath.cuh>
2   -
3   -/*void blur(float* image, float sigma, unsigned int x, unsigned int y){
4   -
5   - stim::cuda::cpu_gaussian_blur_2d<float>(image, sigma, x, y);
6   -}*/
  2 +#include <stim/cuda/templates/conv2.cuh>
  3 +#include <stim/cuda/templates/conv2sep.cuh>
  4 +#include <stim/cuda/templates/chi_gradient.cuh>
7 5  
8 6 void array_multiply(float* lhs, float rhs, unsigned int N){
9 7  
... ... @@ -16,6 +14,12 @@ void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N){
16 14  
17 15 }
18 16  
  17 +void array_multiply(float* ptr1, float* ptr2, float* product, unsigned int N){
  18 +
  19 + stim::cuda::cpu_add(ptr1, ptr2, product, N);
  20 +
  21 +}
  22 +
19 23 void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M){
20 24  
21 25 stim::cuda::cpu_conv2(img, mask, cpu_copy, w, h, M);
... ... @@ -28,3 +32,41 @@ void array_abs(float* img, unsigned int N){
28 32  
29 33 }
30 34  
  35 +void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1){
  36 +
  37 + stim::cuda::cpu_conv2sep(img, x, y, kernel0, k0, kernel1, k1);
  38 +
  39 +}
  40 +
  41 +void array_cos(float* ptr1, float* cpu_out, unsigned int N){
  42 +
  43 + stim::cuda::cpu_cos(ptr1, cpu_out, N);
  44 +
  45 +}
  46 +
  47 +void array_sin(float* ptr1, float* cpu_out, unsigned int N){
  48 +
  49 + stim::cuda::cpu_sin(ptr1, cpu_out, N);
  50 +
  51 +}
  52 +
  53 +void array_atan(float* ptr1, float* cpu_out, unsigned int N){
  54 +
  55 + stim::cuda::cpu_atan(ptr1, cpu_out, N);
  56 +
  57 +}
  58 +
  59 +void array_divide(float* ptr1, float* ptr2,float* cpu_quotient, unsigned int N){
  60 +
  61 + stim::cuda::cpu_divide(ptr1, ptr2, cpu_quotient, N);
  62 +
  63 +}
  64 +
  65 +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){
  66 +
  67 + stim::cuda::cpu_chi_grad(img, cpu_copy, w, h, r, bin_n, bin_size, theta);
  68 +
  69 +}
  70 +
  71 +
  72 +
... ...
gauss_deriv1.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 generates the first-order gaussian derivative filter gx gy,
  7 +/// convolves the image with gx gy,
  8 +/// and returns an image class which channel(0) is Ix and channel(1) is Iy
  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 +
  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> Gd1(stim::image<float> image, int r, unsigned int sigma_n){
  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 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter)
  23 + float sigma = float(r)/float(sigma_n); // calculate the sigma used in gaussian function
  24 +
  25 + stim::image<float> I(w, h, 1, 2); // allocate space for return image class
  26 + stim::image<float> Ix(w, h); // allocate space for Ix
  27 + stim::image<float> Iy(w, h); // allocate space for Iy
  28 + Ix = image; // initialize Ix
  29 + Iy = image; // initialize Iy
  30 +
  31 + float* array_x1;
  32 + array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gx
  33 + float* array_y1;
  34 + array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gx
  35 + float* array_x2;
  36 + array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gy
  37 + float* array_y2;
  38 + array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gy
  39 +
  40 +
  41 + for (int i = 0; i < winsize; i++){
  42 +
  43 + int x = i - r; //range of x
  44 + int y = i - r; //range of y
  45 +
  46 + // create the 1D x-oriented gaussian derivative filter array_x1 for gx
  47 + array_x1[i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  48 + // create the 1D y-oriented gaussian derivative filter array_y1 for gx
  49 + array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  50 + // create the 1D x-oriented gaussian derivative filter array_x2 for gy
  51 + array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  52 + // create the 1D y-oriented gaussian derivative filter array_y2 for gy
  53 + array_y2[i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  54 + }
  55 +
  56 + //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  57 + //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  58 + //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  59 + //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  60 +
  61 + // get Ix by convolving the image with gx
  62 + conv2_sep(Ix.data(), w, h, array_x1, winsize, array_y1, winsize);
  63 +
  64 + //stim::cpu2image(Ix.data(), "data_output/Ix_0915.bmp", w, h, stim::cmBrewer);
  65 + // get Iy by convolving the image with gy
  66 + conv2_sep(Iy.data(), w, h, array_x2, winsize, array_y2, winsize);
  67 +
  68 + //stim::cpu2image(Iy.data(), "data_output/Iy_0915.bmp", w, h, stim::cmBrewer);
  69 +
  70 + delete [] array_x1; //free the memory of array_x1
  71 + delete [] array_y1; //free the memory of array_y1
  72 + delete [] array_x2; //free the memory of array_x2
  73 + delete [] array_y2; //free the memory of array_y2
  74 +
  75 + I.set_channel(0, Ix.data());
  76 + I.set_channel(1, Iy.data());
  77 +
  78 + return I;
  79 +
  80 +}
0 81 \ No newline at end of file
... ...
gauss_deriv2.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 generates the second-order gaussian derivative filter gxx gyy,
  7 +/// convolves the image with gxx gyy,
  8 +/// and returns an image class which channel(0) is Ixx and channel(1) is Iyy
  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 +
  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> Gd2(stim::image<float> image, int r, unsigned int sigma_n){
  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 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter)
  23 + float sigma = float(r)/float(sigma_n); // calculate the sigma used in gaussian function
  24 +
  25 + stim::image<float> I(w, h, 1, 2); // allocate space for return image class
  26 + stim::image<float> Ixx(w, h); // allocate space for Ixx
  27 + stim::image<float> Iyy(w, h); // allocate space for Iyy
  28 + Ixx = image; // initialize Ixx
  29 + Iyy = image; // initialize Iyy
  30 +
  31 + float* array_x1;
  32 + array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gxx
  33 + float* array_y1;
  34 + array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gxx
  35 + float* array_x2;
  36 + array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gyy
  37 + float* array_y2;
  38 + array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gyy
  39 +
  40 +
  41 + for (int i = 0; i < winsize; i++){
  42 +
  43 + int x = i - r; //range of x
  44 + int y = i - r; //range of y
  45 +
  46 + // create the 1D x-oriented gaussian derivative filter array_x1 for gxx
  47 + array_x1[i] = (-1) * (1 - pow(x, 2)) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  48 + // create the 1D y-oriented gaussian derivative filter array_y1 for gxx
  49 + array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  50 + // create the 1D x-oriented gaussian derivative filter array_x2 for gyy
  51 + array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  52 + // create the 1D y-oriented gaussian derivative filter array_y2 for gyy
  53 + array_y2[i] = (-1) * (1 - pow(y, 2)) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  54 + }
  55 +
  56 + //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  57 + //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  58 + //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  59 + //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  60 +
  61 + // get Ixx by convolving the image with gxx
  62 + conv2_sep(Ixx.data(), w, h, array_x1, winsize, array_y1, winsize);
  63 +
  64 + //stim::cpu2image(Ixx.data(), "data_output/Ixx_0915.bmp", w, h, stim::cmBrewer);
  65 + // get Iyy by convolving the image with gyy
  66 + conv2_sep(Iyy.data(), w, h, array_x2, winsize, array_y2, winsize);
  67 +
  68 + //stim::cpu2image(Iyy.data(), "data_output/Iyy_0915.bmp", w, h, stim::cmBrewer);
  69 +
  70 + delete [] array_x1; //free the memory of array_x1
  71 + delete [] array_y1; //free the memory of array_y1
  72 + delete [] array_x2; //free the memory of array_x2
  73 + delete [] array_y2; //free the memory of array_y2
  74 +
  75 + I.set_channel(0, Ixx.data());
  76 + I.set_channel(1, Iyy.data());
  77 +
  78 + return I;
  79 +
  80 +}
0 81 \ No newline at end of file
... ...
gauss_deriv_center.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 +
  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 array_abs(float* img, unsigned int N);
  11 +
  12 +/// This function evaluates the center-surround(Laplacian of Gaussian) gaussian derivative gradient of an one-channel image
  13 +
  14 +/// @param img is the one-channel image
  15 +/// @param r is an array of radii for different scaled discs(filters)
  16 +/// @param sigma_n is the number of standard deviations used to define the sigma
  17 +
  18 +stim::image<float> Gd_center(stim::image<float> image, int r, unsigned int sigma_n){
  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 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter)
  24 +
  25 + stim::image<float> I(w, h, 1, 2); // allocate space for return image of Gd2
  26 + stim::image<float> Ixx(w, h); // allocate space for Ixx
  27 + stim::image<float> Iyy(w, h); // allocate space for Iyy
  28 + stim::image<float> Gd_center(w, h); // allocate space for Pb
  29 +
  30 + I = Gd2(image, r, sigma_n); // calculate the Ixx, Iyy
  31 + Ixx = I.channel(0);
  32 + Iyy = I.channel(1);
  33 +
  34 + array_add(Ixx.data(), Iyy.data(), Gd_center.data(), N); //Gd_center = Ixx + Iyy;
  35 + array_abs(Gd_center.data(), N);
  36 +
  37 + //stim::cpu2image(Gd_center.data(), "data_output/Gd_center_0919.bmp", w, h, stim::cmBrewer);
  38 +
  39 + return Gd_center;
  40 +
  41 +}
... ...
gauss_deriv_even.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 +}
0 53 \ No newline at end of file
... ...
gauss_deriv_odd.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 +
  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 array_abs(float* img, unsigned int N);
  11 +
  12 +/// This function evaluates the theta-dependent odd symmetric gaussian derivative gradient of an one-channel image
  13 +
  14 +/// @param img is the one-channel image
  15 +/// @param r is an array of radii for different scaled discs(filters)
  16 +/// @param sigma_n is the number of standard deviations used to define the sigma
  17 +/// @param theta is angle used for computing the gradient
  18 +
  19 +stim::image<float> Gd_odd(stim::image<float> image, int r, unsigned int sigma_n, float theta){
  20 +
  21 + float theta_r = (theta * PI)/180; //change angle unit from degree to rad
  22 +
  23 + unsigned int w = image.width(); // get the width of picture
  24 + unsigned int h = image.height(); // get the height of picture
  25 + unsigned N = w * h; // get the number of pixels of picture
  26 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter)
  27 +
  28 + stim::image<float> I(w, h, 1, 2); // allocate space for return image of Gd1
  29 + stim::image<float> Ix(w, h); // allocate space for Ix
  30 + stim::image<float> Iy(w, h); // allocate space for Iy
  31 + stim::image<float> Gd_odd_theta(w, h); // allocate space for Pb
  32 +
  33 + I = Gd1(image, r, sigma_n); // calculate the Ix, Iy
  34 + Ix = I.channel(0);
  35 + Iy = I.channel(1);
  36 +
  37 + array_multiply(Ix.data(), cos(theta_r), N); //Ix = Ix*cos(theta_r)
  38 + array_multiply(Iy.data(), sin(theta_r), N); //Iy = Iy*sin(theta_r)
  39 + array_add(Ix.data(), Iy.data(), Gd_odd_theta.data(), N); //Gd_odd_theta = Ix + Iy;
  40 + array_abs(Gd_odd_theta.data(), N);
  41 +
  42 + //stim::cpu2image(I.channel(0).data(), "data_output/Gd_odd_x_0919.bmp", w, h, stim::cmBrewer);
  43 + //stim::cpu2image(I.channel(1).data(), "data_output/Gd_odd_y_0919.bmp", w, h, stim::cmBrewer);
  44 + //stim::cpu2image(Gd_odd_theta.data(), "data_output/Gd_odd_theta_0919.bmp", w, h, stim::cmBrewer);
  45 +
  46 + return Gd_odd_theta;
  47 +
  48 +}
... ...
kmeans.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 +}
0 68 \ No newline at end of file
... ...
main.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 +
  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 + image.load("101085.bmp"); // load the input image
  18 + //image.load("slice00_500_500.bmp"); // load the input image
  19 + image = image.channel(0); // get the first channel of black-and-white image
  20 +
  21 +
  22 + unsigned int w = image.width(); // get the width of picture
  23 + unsigned int h = image.height(); // get the height of picture
  24 + unsigned int N = w * h; // get the number of pixels
  25 + int c = image.channels(); // get the number if channels of picture
  26 + int s = 3; // set the number of scales
  27 + int r1[3] = {3,5,10}; // set an array of radii for different scaled discs(filters) for cPb, the length of r is c*s
  28 + int r2[3] = {5,10,20}; // set an array of radii for different scaled discs(filters) for tPb, the length of r is c*s
  29 + float alpha[3] = {1,1,1}; // set an array of weights for different scaled discs(filters)
  30 + unsigned int theta_n = 8; // set the number of angles used in filter orientation
  31 + unsigned int bin_n = 16; // set the number of bins used in chi-distance
  32 + unsigned int K = 16; // set the number of cludters, K should be multiple of bin_n
  33 +
  34 + stim::image<float> img_cPb(w, h); // allocate the space for cPb
  35 + stim::image<float> img_tPb(w, h); // allocate the space for tPb
  36 + stim::image<float> img_mPb(w, h); // allocate the space for tPb
  37 +
  38 + std::cout<<"imagesize: "<< w <<"*"<< h <<'\n';
  39 +
  40 + //*******************cPb********************
  41 + std::cout<<"begin cPb"<<'\n';
  42 +
  43 +
  44 + std::clock_t start1; // (optional) set the timer to calculate the total time
  45 + start1 = std::clock(); // (optional) set timer start point
  46 +
  47 + img_cPb = cPb(image, r1, alpha, s);
  48 +
  49 + // show the cPb (optional)
  50 + stim::cpu2image(img_cPb.data(), "data_output/img_cPb_0925.bmp", w, h, stim::cmBrewer);
  51 +
  52 + double duration1 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
  53 + std::cout<<"cPb time: "<< duration1 <<"s"<<'\n'; // (optional) show the total time
  54 +
  55 +
  56 + //*******************tPb********************
  57 + std::cout<<"begin tPb"<<'\n';
  58 +
  59 +
  60 + std::clock_t start2; // (optional) set the timer to calculate the total time
  61 + start2 = std::clock(); // (optional) set timer start point
  62 +
  63 + img_tPb = tPb(image, r2, alpha, theta_n, bin_n, s, K);
  64 +
  65 + // show the tPb (optional)
  66 + stim::cpu2image(img_tPb.data(), "data_output/img_tPb_0925.bmp", w, h, stim::cmBrewer);
  67 +
  68 + double duration2 = ( std::clock() - start2 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
  69 + std::cout<<"tPb time: "<< duration2 <<"s"<<'\n'; // (optional) show the total time
  70 +
  71 +
  72 + //*******************************************
  73 +
  74 + double duration3 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
  75 + std::cout<<"total running time: "<< duration3 <<"s"<<'\n'; // (optional) show the total time
  76 +
  77 + //******************mPb**********************
  78 + // set parameters for linear combination
  79 + float a = 1;
  80 + float b = 0.5;
  81 +
  82 + // create mPb by linearly combined cPb and tPb
  83 + array_multiply(img_cPb.data(), a, N);
  84 + array_multiply(img_tPb.data(), b, N);
  85 + array_add(img_cPb.data(), img_tPb.data(), img_mPb.data(), N);
  86 +
  87 + // show the mPb (optional)
  88 + stim::cpu2image(img_mPb.data(), "data_output/img_mPb_0925.bmp", w, h, stim::cmBrewer);
  89 +
  90 + return 0;
  91 +
  92 +}
... ...
fun_mPb_theta.cpp renamed to old version/fun_mPb_theta.cpp
  Show/Hide comments   View file @6dcc460
... ... @@ -4,6 +4,7 @@
4 4 #include <stim/image/image_contour_detection.h>
5 5 #include <sstream>
6 6  
  7 +
7 8 void array_multiply(float* lhs, float rhs, unsigned int N);
8 9 void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
9 10  
... ...
func_mPb.cpp renamed to old version/func_mPb.cpp
  Show/Hide comments   View file @6dcc460
... ... @@ -11,6 +11,7 @@
11 11 /// @param r is an array of radii for different scaled discs(filters)
12 12 /// @param alpha is an array of weights for different scaled discs(filters)
13 13 /// @param s is the number of scales
  14 +
14 15 stim::image<float> func_mPb(stim::image<float> img, unsigned int theta_n, int* r, float* alpha, int s){
15 16  
16 17 std::clock_t start; // (optional) set the timer to calculate the total time
... ... @@ -34,8 +35,8 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; img, unsigned int theta_n, int* r
34 35  
35 36 for (unsigned int n = 0; n < theta_n; n++){
36 37  
37   - ss << "data_output/mPb_theta"<< n << "_0911.bmp"; // (optional) set the name for test result file
38   - std::string sss = ss.str(); // (optional)
  38 + //ss << "data_output/mPb_theta"<< n << "_0911.bmp"; // (optional) set the name for test result file
  39 + //std::string sss = ss.str(); // (optional)
39 40 float theta = 180 * ((float)n/theta_n); // calculate the even-splited angle for each mPb_theta
40 41  
41 42 mPb_theta = func_mPb_theta(img, theta, r, alpha, s); // calculate the mPb_theta
... ... @@ -50,7 +51,7 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; img, unsigned int theta_n, int* r
50 51 unsigned long idx = n * w * h * 1; //index for the nth mPb_theta
51 52  
52 53  
53   - stim::cpu2image(mPb_theta.data(), sss, w, h, stim::cmBrewer); // (optional) output the nth mPb_theta
  54 + //stim::cpu2image(mPb_theta.data(), sss, w, h, stim::cmBrewer); // (optional) output the nth mPb_theta
54 55  
55 56 for(unsigned long i = 0; i < N; i++){
56 57  
... ... @@ -65,11 +66,11 @@ stim::image&lt;float&gt; func_mPb(stim::image&lt;float&gt; img, unsigned int theta_n, int* r
65 66 }
66 67  
67 68  
68   - ss.str(""); //(optional) clear the space for stream
  69 + //ss.str(""); //(optional) clear the space for stream
69 70  
70 71 }
71 72  
72   - stim::cpu2image(mPb.data(), "data_output/mPb_500_0911_neat.bmp", w, h, stim::cmBrewer); // output the mPb
  73 + stim::cpu2image(mPb.data(), "data_output/mPb_500_0914_neat.bmp", w, h, stim::cmBrewer); // output the mPb
73 74  
74 75 double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time
75 76 std::cout<<"total time:"<< duration2 <<"s"<<'\n'; // (optional) show the total time
... ...
gauss_derivative_odd.cpp renamed to old version/gauss_derivative_odd.cpp
  Show/Hide comments   View file @6dcc460
... ... @@ -12,7 +12,7 @@ void array_multiply(float* lhs, float rhs, unsigned int N);
12 12  
13 13 /// @param img is the one-channel image
14 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
  15 +/// @param sigma_n is the number of standard deviations used to define the sigma
16 16 /// @param theta is angle used for computing the gradient
17 17  
18 18 stim::image<float> gaussian_derivative_filter_odd(stim::image<float> image, int r, unsigned int sigma_n, float theta){
... ... @@ -28,8 +28,12 @@ stim::image&lt;float&gt; gaussian_derivative_filter_odd(stim::image&lt;float&gt; image, int
28 28 stim::image<float> mask_theta(winsize, winsize);// allocate space for theta-oriented filter
29 29 stim::image<float> derivative_theta(w, h); // allocate space for theta-oriented gradient
30 30  
  31 + stim::image<float> mask_r;
  32 +
  33 +
31 34 float theta_r = (theta * PI)/180; //change angle unit from degree to rad
32 35  
  36 + //*************inseparable convolution****************
33 37 for (int j = 0; j < winsize; j++){
34 38 for (int i = 0; i< winsize; i++){
35 39  
... ... @@ -46,14 +50,37 @@ stim::image&lt;float&gt; gaussian_derivative_filter_odd(stim::image&lt;float&gt; image, int
46 50 }
47 51 }
48 52  
  53 + mask_r = mask_x.rotate(90, r, r);
  54 + stim::cpu2image(mask_r.data(), "data_output/mask_r90_0915.bmp", winsize, winsize, stim::cmBrewer);
  55 + //stim::cpu2image(mask_theta.data(), "data_output/mask_0911_2.bmp", winsize, winsize, stim::cmBrewer); // (optional) show the mask result
  56 +
  57 + // do the 2D convolution with image and mask
  58 + conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize);
  59 + //*********************************************************/
  60 +
  61 + /*************separable convolution****************
  62 + for (int i = 0; i < winsize; i++){
  63 +
  64 + int x = i - r; //range of x
  65 + int y = i - r; //range of y
  66 +
  67 + // create the 1D x-oriented gaussian derivative filter array_x
  68 + array_x[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  69 + // create the 1D y-oriented gaussian derivative filter array_y
  70 + array_y[i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  71 +
  72 +
  73 + }
  74 +
49 75 //stim::cpu2image(mask_theta.data(), "data_output/mask_0911_2.bmp", winsize, winsize, stim::cmBrewer); // (optional) show the mask result
50 76  
51 77 // do the 2D convolution with image and mask
52 78 conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize);
  79 + //*********************************************************/
53 80  
54 81 array_abs(derivative_theta.data(), N); // get the absolute value for each pixel (why slower than the "for loop" method sometimes?)
55 82  
56   - float max = derivative_theta.max(); // get the maximum of gradient used for normalization
  83 + float max = derivative_theta.maxv(); // get the maximum of gradient used for normalization
57 84 array_multiply(derivative_theta.data(), 1/max, N); // normalize the gradient
58 85  
59 86  
... ...
test_main.cpp renamed to old version/test_main - Copy.cpp
  Show/Hide comments   View file @6dcc460
... ... @@ -3,19 +3,27 @@
3 3 #include <stim/visualization/colormap.h>
4 4 #include <stim/image/image_contour_detection.h>
5 5 #include <iostream>
  6 +#include <stim/visualization/colormap.h>
  7 +
6 8 /// calculate the mPb given a multi-channel image
7 9  
8 10 int main()
9 11 {
10 12 stim::image<float> img; // generate an image object
  13 + //stim::image<float> img_r;
11 14  
12   - img.load("slice00_500_500.bmp"); // load the input image
  15 + //img.load("slice00_500_500.bmp"); // load the input image
  16 + img.load("101087_square.bmp"); // load the input image
13 17 img = img.channel(0); // get the first channel of black-and-white image
14 18  
15 19 unsigned int w = img.width(); // get the width of picture
16 20 unsigned int h = img.height(); // get the height of picture
17 21 int c = img.channels(); // get the number if channels of picture
18 22 int s = 3; // set the number of scales
  23 +
  24 + //img_r = img.rotate(45, w/2, h/2);
  25 + //stim::cpu2image(img_r.data(), "data_output/rotate_0915.bmp", w, h, stim::cmBrewer);
  26 +
19 27  
20 28 int r[3] = {3,5,10}; // set an array of radii for different scaled discs(filters)
21 29 float alpha[3] = {1,1,1}; // set an array of weights for different scaled discs(filters)
... ...
tPb.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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_0925.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 +}
... ...
textons.cpp 0 → 100644
  Show/Hide comments   View file @6dcc460
  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 + unsigned int r_odd = 3; // set disc radii for odd-symmetric filter
  23 + unsigned int sigma_n_odd = 3; // set sigma_n for odd-symmetric filter
  24 + unsigned int r_even = 3; // set disc radii for even-symmetric filter
  25 + unsigned int sigma_n_even = 3; // set sigma_n for even-symmetric filter
  26 + unsigned int r_center = 3; // set disc radii for center-surround filter
  27 + unsigned int sigma_n_center = 3; // set sigma_n for center-surround filter
  28 +
  29 + //std::ostringstream ss1, ss2; // (optional) set the stream to designate the test result file
  30 +
  31 + for (unsigned int i = 0; i < theta_n; i++){
  32 +
  33 + //ss1 << "data_output/textons_channel_"<< i << ".bmp"; // set the name for test result file (optional)
  34 + //std::string sss1 = ss1.str();
  35 + //ss2 << "data_output/textons_channel_"<< i+theta_n << ".bmp"; // set the name for test result file (optional)
  36 + //std::string sss2 = ss2.str();
  37 +
  38 + float theta = 180 * ((float)i/theta_n); // calculate the even-splited angle for each oriented filter
  39 +
  40 + temp = Gd_odd(image, r_odd, sigma_n_odd, theta); // return Gd_odd to temp
  41 + //stim::cpu2image(temp.data(), sss1, w, h, stim::cmBrewer);
  42 + textons.set_channel(i, temp.data()); // copy temp to ith channel of textons
  43 +
  44 + temp = Gd_even(image, r_even, sigma_n_even, theta); // return Gd_even to temp
  45 + //stim::cpu2image(temp.data(), sss2, w, h, stim::cmBrewer);
  46 + textons.set_channel(i + theta_n, temp.data()); // copy temp to (i+theta_n)th channel of textons
  47 +
  48 + //ss1.str(""); //(optional) clear the space for stream
  49 + //ss2.str(""); //(optional) clear the space for stream
  50 +
  51 + }
  52 +
  53 + temp = Gd_center(image, r_center, sigma_n_center); // return Gd_center to temp
  54 + //stim::cpu2image(temp.data(), "data_output/textons_channel_16.bmp", w, h, stim::cmBrewer);
  55 + textons.set_channel(theta_n*2, temp.data()); // copy temp to (theta_n*2)th channel of textons
  56 +
  57 + return textons;
  58 +
  59 +}
  60 +
  61 +
0 62 \ No newline at end of file
... ...