diff --git a/stim/cuda/arraymath/array_multiply.cuh b/stim/cuda/arraymath/array_multiply.cuh index 4385e77..39a94d2 100644 --- a/stim/cuda/arraymath/array_multiply.cuh +++ b/stim/cuda/arraymath/array_multiply.cuh @@ -54,6 +54,63 @@ namespace stim{ //free allocated memory cudaFree(gpuLHS); } + + +// array .* array multiplication + + template + __global__ void cuda_multiply(T* ptr1, T* ptr2, T* product, unsigned int N){ + + //calculate the 1D index for this thread + int idx = blockIdx.x * blockDim.x + threadIdx.x; + + if(idx < N){ + product[idx] = ptr1[idx] * ptr2[idx]; + } + + } + + template + void gpu_multiply(T* ptr1, T* ptr2, T* product, unsigned int N){ + + //get the maximum number of threads per block for the CUDA device + int threads = stim::maxThreadsPerBlock(); + + //calculate the number of blocks + int blocks = N / threads + 1; + + //call the kernel to do the multiplication + cuda_multiply <<< blocks, threads >>>(ptr1, ptr2, product, N); + + } + + template + void cpu_multiply(T* ptr1, T* ptr2, T* cpu_product, unsigned int N){ + + //allocate memory on the GPU for the array + T* gpu_ptr1; + T* gpu_ptr2; + T* gpu_product; + HANDLE_ERROR( cudaMalloc( &gpu_ptr1, N * sizeof(T) ) ); + HANDLE_ERROR( cudaMalloc( &gpu_ptr2, N * sizeof(T) ) ); + HANDLE_ERROR( cudaMalloc( &gpu_product, N * sizeof(T) ) ); + + //copy the array to the GPU + HANDLE_ERROR( cudaMemcpy( gpu_ptr1, ptr1, N * sizeof(T), cudaMemcpyHostToDevice) ); + HANDLE_ERROR( cudaMemcpy( gpu_ptr2, ptr2, N * sizeof(T), cudaMemcpyHostToDevice) ); + + //call the GPU version of this function + gpu_multiply(gpu_ptr1, gpu_ptr2 ,gpu_product, N); + + //copy the array back to the CPU + HANDLE_ERROR( cudaMemcpy( cpu_product, gpu_product, N * sizeof(T), cudaMemcpyDeviceToHost) ); + + //free allocated memory + cudaFree(gpu_ptr1); + cudaFree(gpu_ptr2); + cudaFree(gpu_product); + + } } } diff --git a/stim/cuda/bsds500/Pb.cpp b/stim/cuda/bsds500/Pb.cpp deleted file mode 100644 index c4ec506..0000000 --- a/stim/cuda/bsds500/Pb.cpp +++ /dev/null @@ -1,65 +0,0 @@ -#include -//#include -#include -#include - -#define SIGMA_N 3 - -void array_abs(float* img, unsigned int N); -void array_multiply(float* lhs, float rhs, unsigned int N); -void array_cos(float* ptr1, float* cpu_out, unsigned int N); -void array_sin(float* ptr1, float* cpu_out, unsigned int N); -void array_atan(float* ptr1, float* cpu_out, unsigned int N); -void array_divide(float* ptr1, float* ptr2,float* cpu_quotient, unsigned int N); -void array_multiply(float* ptr1, float* ptr2, float* product, unsigned int N); -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N); - -/// This function uses odd-symmetric gaussian derivative filter to evaluate -/// the max probability of a contour on one scale, given an one-channel image - -/// @param img is an one-channel image -/// @param r is an array of radii for different scaled discs(filters) -/// @param sigma_n is the number of standard deviations used to define the sigma - -stim::image Pb(stim::image image, int sigma){ - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned N = w * h; // get the number of pixels of picture - - int r = SIGMA_N * sigma; // set the radius of filter - int winsize = 2 * r + 1; // set the winsdow size of filter - - stim::image I(w, h, 1, 2); // allocate space for return image of dG1_conv2 - stim::image theta(w, h); // allocate space for theta matrix - stim::image cos(w, h); // allocate space for cos(theta) - stim::image sin(w, h); // allocate space for sin(theta) - stim::image temp(w, h); // allocate space for temp - stim::image Ix(w, h); // allocate space for Ix - stim::image Iy(w, h); // allocate space for Iy - stim::image Pb(w, h); // allocate space for Pb - - I = dG1_conv2(image, sigma); // calculate the Ix, Iy - Ix = I.channel(0); - array_abs(Ix.data(), N); //get |Ix|; - //stim::cpu2image(Ix.data(), "data_output/Pb_Ix_0924.bmp", w, h, stim::cmBrewer); - Iy = I.channel(1); - array_abs(Iy.data(), N); //get |Iy|; - //stim::cpu2image(Iy.data(), "data_output/Pb_Iy_0924.bmp", w, h, stim::cmBrewer); - - array_divide(Iy.data(), Ix.data(), temp.data(), N); //temp = Iy./Ix - array_atan(temp.data(), theta.data(), N); //theta = atan(temp) - array_cos(theta.data(), cos.data(), N); //cos = cos(theta) - array_sin(theta.data(), sin.data(), N); //sin = sin(theta) - array_multiply(Ix.data(), cos.data(), Ix.data(), N); //Ix = Ix.*cos - array_multiply(Iy.data(), sin.data(), Iy.data(), N); //Iy = Iy.*sin - array_add(Ix.data(), Iy.data(), Pb.data(), N); //Pb = Ix + Iy; - - float max = Pb.maxv(); // get the maximum of Pb used for normalization - array_multiply(Pb.data(), 1/max, N); // normalize the Pb - - //stim::cpu2image(Pb.data(), "data_output/Pb_0924.bmp", w, h, stim::cmBrewer); show the Pb(optional) - - return Pb; - -} \ No newline at end of file diff --git a/stim/cuda/bsds500/cPb.cpp b/stim/cuda/bsds500/cPb.cpp deleted file mode 100644 index f8f1388..0000000 --- a/stim/cuda/bsds500/cPb.cpp +++ /dev/null @@ -1,66 +0,0 @@ -#include -#include -#include -#include -#include - - -void array_multiply(float* lhs, float rhs, unsigned int N); -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N); - -/// This function evaluates the cPb given an multi-channel image - -/// @param img is the multi-channel image -/// @param r is an array of radii for different scaled discs(filters) -/// @param alpha is is an array of weights for different scaled discs(filters) -/// @param s is the number of scales - -stim::image cPb(stim::image img, int* sigma, float* alpha, int s){ - - unsigned int w = img.width(); // get the width of picture - unsigned int h = img.height(); // get the height of picture - unsigned int c = img.channels(); // get the channels of picture - - - stim::image cPb(w, h, 1); // allocate space for cPb - unsigned size = cPb.size(); // get the size of cPb - memset ( cPb.data(), 0, size * sizeof(float)); // initialize all the pixels of cPb to 0 - - - unsigned int N = w * h; // get the number of pixels - - //std::ostringstream ss; // (optional) set the stream to designate the test result file - - stim::image temp; // set the temporary image to store the addtion result - - for (int i = 0; i < c; i++){ - for (int j = 0; j < s; j++){ - - //ss << "data_output/cPb_slice"<< i*s + j << ".bmp"; // set the name for test result file (optional) - //std::string sss = ss.str(); - - // get the gaussian gradient by convolving each image slice with the mask - temp = Pb(img.channel(i), sigma[i*s + j]); - - // output the test result of each slice (optional) - //stim::cpu2image(temp.data(), sss, w, h, stim::cmBrewer); - - // multiply each gaussian gradient with its weight - array_multiply(temp.data(), alpha[i*s + j], N); - - // add up all the weighted gaussian gradients - array_add(cPb.data(), temp.data(), cPb.data(), N); - - //ss.str(""); //(optional) clear the space for stream - - } - } - - float max = cPb.maxv(); // get the maximum of cPb used for normalization - array_multiply(cPb.data(), 1/max, N); // normalize the cPb - - // output the test result of cPb (optional) - //stim::cpu2image(cPb.data(), "data_output/cPb_0916.bmp", w, h, stim::cmBrewer); - - return cPb; -} diff --git a/stim/cuda/bsds500/dG1_conv2.cpp b/stim/cuda/bsds500/dG1_conv2.cpp deleted file mode 100644 index 4e98709..0000000 --- a/stim/cuda/bsds500/dG1_conv2.cpp +++ /dev/null @@ -1,81 +0,0 @@ -#include -//#include -#include -#include -#define SIGMA_N 3 - -/// This function generates the first-order gaussian derivative filter gx gy, -/// convolves the image with gx gy, -/// and returns an image class which channel(0) is Ix and channel(1) is Iy - -/// @param img is the one-channel image -/// @param sigma is the parameter for gaussian function - -void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1); -//void array_abs(float* img, unsigned int N); - -stim::image dG1_conv2(stim::image image, int sigma){ - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned N = w * h; // get the number of pixels of picture - - int r = SIGMA_N * sigma; - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter - - stim::image I(w, h, 1, 2); // allocate space for return image class - stim::image Ix(w, h); // allocate space for Ix - stim::image Iy(w, h); // allocate space for Iy - Ix = image; // initialize Ix - Iy = image; // initialize Iy - - float* array_x1; - array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gx - float* array_y1; - array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gx - float* array_x2; - array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gy - float* array_y2; - array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gy - - - for (int i = 0; i < winsize; i++){ - - int x = i - r; //range of x - int y = i - r; //range of y - - // create the 1D x-oriented gaussian derivative filter array_x1 for gx - array_x1[i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); - // create the 1D y-oriented gaussian derivative filter array_y1 for gx - array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); - // create the 1D x-oriented gaussian derivative filter array_x2 for gy - array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); - // create the 1D y-oriented gaussian derivative filter array_y2 for gy - array_y2[i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); - } - - //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - - // get Ix by convolving the image with gx - conv2_sep(Ix.data(), w, h, array_x1, winsize, array_y1, winsize); - - //stim::cpu2image(Ix.data(), "data_output/Ix_0915.bmp", w, h, stim::cmBrewer); - // get Iy by convolving the image with gy - conv2_sep(Iy.data(), w, h, array_x2, winsize, array_y2, winsize); - - //stim::cpu2image(Iy.data(), "data_output/Iy_0915.bmp", w, h, stim::cmBrewer); - - delete [] array_x1; //free the memory of array_x1 - delete [] array_y1; //free the memory of array_y1 - delete [] array_x2; //free the memory of array_x2 - delete [] array_y2; //free the memory of array_y2 - - I.set_channel(0, Ix.data()); - I.set_channel(1, Iy.data()); - - return I; - -} \ No newline at end of file diff --git a/stim/cuda/bsds500/dG1_theta_conv2.cpp b/stim/cuda/bsds500/dG1_theta_conv2.cpp deleted file mode 100644 index 99c29fb..0000000 --- a/stim/cuda/bsds500/dG1_theta_conv2.cpp +++ /dev/null @@ -1,52 +0,0 @@ -#include -#include -#include -#include - -#define PI 3.1415926 -#define SIGMA_N 3 - -void array_multiply(float* lhs, float rhs, unsigned int N); -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N); -void array_abs(float* img, unsigned int N); - -/// This function evaluates the theta-dependent odd symmetric gaussian derivative gradient of an one-channel image - -/// @param img is the one-channel image -/// @param r is an array of radii for different scaled discs(filters) -/// @param sigma_n is the number of standard deviations used to define the sigma -/// @param theta is angle used for computing the gradient - -stim::image dG1_theta_conv2(stim::image image, int sigma, float theta){ - - float theta_r = (theta * PI)/180; //change angle unit from degree to rad - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - - int r = SIGMA_N * sigma; - unsigned N = w * h; // get the number of pixels of picture - - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter - - stim::image I(w, h, 1, 2); // allocate space for return image of dG1_conv2 - stim::image Ix(w, h); // allocate space for Ix - stim::image Iy(w, h); // allocate space for Iy - stim::image dG1_theta(w, h); // allocate space for Pb - - I = dG1_conv2(image, sigma); // calculate the Ix, Iy - Ix = I.channel(0); - Iy = I.channel(1); - - array_multiply(Ix.data(), cos(theta_r), N); //Ix = Ix*cos(theta_r) - array_multiply(Iy.data(), sin(theta_r), N); //Iy = Iy*sin(theta_r) - array_add(Ix.data(), Iy.data(), dG1_theta.data(), N); //dG1_theta = Ix + Iy; - array_abs(dG1_theta.data(), N); - - //stim::cpu2image(I.channel(0).data(), "data_output/dG1_theta_x_0919.bmp", w, h, stim::cmBrewer); - //stim::cpu2image(I.channel(1).data(), "data_output/dG1_theta_y_0919.bmp", w, h, stim::cmBrewer); - //stim::cpu2image(dG1_theta.data(), "data_output/dG1_theta_0919.bmp", w, h, stim::cmBrewer); - - return dG1_theta; - -} diff --git a/stim/cuda/bsds500/dG2_conv2.cpp b/stim/cuda/bsds500/dG2_conv2.cpp deleted file mode 100644 index cb92322..0000000 --- a/stim/cuda/bsds500/dG2_conv2.cpp +++ /dev/null @@ -1,81 +0,0 @@ -#include -//#include -#include -#include -#define SIGMA_N 3 - -/// This function generates the second-order gaussian derivative filter gxx gyy, -/// convolves the image with gxx gyy, -/// and returns an image class which channel(0) is Ixx and channel(1) is Iyy - -/// @param img is the one-channel image -/// @param sigma is the parameter for gaussian function - -void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1); -//void array_abs(float* img, unsigned int N); - -stim::image dG2_conv2(stim::image image, int sigma){ - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned N = w * h; // get the number of pixels of picture - - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter - int r = SIGMA_N * sigma; - - stim::image I(w, h, 1, 2); // allocate space for return image class - stim::image Ixx(w, h); // allocate space for Ixx - stim::image Iyy(w, h); // allocate space for Iyy - Ixx = image; // initialize Ixx - Iyy = image; // initialize Iyy - - float* array_x1; - array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gxx - float* array_y1; - array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gxx - float* array_x2; - array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gyy - float* array_y2; - array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gyy - - - for (int i = 0; i < winsize; i++){ - - int x = i - r; //range of x - int y = i - r; //range of y - - // create the 1D x-oriented gaussian derivative filter array_x1 for gxx - array_x1[i] = (-1) * (1 - pow(x, 2)) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); - // create the 1D y-oriented gaussian derivative filter array_y1 for gxx - array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); - // create the 1D x-oriented gaussian derivative filter array_x2 for gyy - array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); - // create the 1D y-oriented gaussian derivative filter array_y2 for gyy - array_y2[i] = (-1) * (1 - pow(y, 2)) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); - } - - //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result - - // get Ixx by convolving the image with gxx - conv2_sep(Ixx.data(), w, h, array_x1, winsize, array_y1, winsize); - - //stim::cpu2image(Ixx.data(), "data_output/Ixx_0915.bmp", w, h, stim::cmBrewer); - // get Iyy by convolving the image with gyy - conv2_sep(Iyy.data(), w, h, array_x2, winsize, array_y2, winsize); - - //stim::cpu2image(Iyy.data(), "data_output/Iyy_0915.bmp", w, h, stim::cmBrewer); - - delete [] array_x1; //free the memory of array_x1 - delete [] array_y1; //free the memory of array_y1 - delete [] array_x2; //free the memory of array_x2 - delete [] array_y2; //free the memory of array_y2 - - I.set_channel(0, Ixx.data()); - I.set_channel(1, Iyy.data()); - - return I; - -} \ No newline at end of file diff --git a/stim/cuda/bsds500/dG2_d2x_theta_conv2.cpp b/stim/cuda/bsds500/dG2_d2x_theta_conv2.cpp deleted file mode 100644 index 13ce010..0000000 --- a/stim/cuda/bsds500/dG2_d2x_theta_conv2.cpp +++ /dev/null @@ -1,55 +0,0 @@ -#include -#include -#include -#include - -#define SIGMA_N 3 - -/// This function evaluates the theta-dependent even-symmetric gaussian derivative gradient of an one-channel image - -/// @param img is the one-channel image -/// @param r is an array of radii for different scaled discs(filters) -/// @param sigma_n is the number of standard deviations used to define the sigma -/// @param theta is angle used for computing the gradient - -void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M); -void array_abs(float* img, unsigned int N); - -stim::image dG2_d2x_theta_conv2(stim::image image, int sigma, float theta){ - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned N = w * h; // get the number of pixels of picture - - int r = SIGMA_N * sigma; // set the radius of filter - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter - - stim::image I(w, h, 1, 2); // allocate space for return image class - stim::image dG2_d2x_theta(w, h); // allocate space for dG2_d2x_theta - stim::image mask_x(winsize, winsize); // allocate space for x-axis-oriented filter - stim::image mask_r(winsize, winsize); // allocate space for theta-oriented filter - - for (int j = 0; j < winsize; j++){ - for (int i = 0; i< winsize; i++){ - - int x = i - r; //range of x - int y = j - r; //range of y - - // create the x-oriented gaussian derivative filter mask_x - 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))); - - } - } - - mask_r = mask_x.rotate(theta, r, r); - //mask_r = mask_x.rotate(45, r, r); - //stim::cpu2image(mask_r.data(), "data_output/mask_r_0919.bmp", winsize, winsize, stim::cmBrewer); - - // do the 2D convolution with image and mask - conv2(image.data(), mask_r.data(), dG2_d2x_theta.data(), w, h, winsize); - array_abs(dG2_d2x_theta.data(), N); - - //stim::cpu2image(dG2_d2x_theta.data(), "data_output/dG2_d2x_theta_0919.bmp", w, h, stim::cmGrayscale); - - return dG2_d2x_theta; -} \ No newline at end of file diff --git a/stim/cuda/bsds500/dG_d2x_theta_conv2.cpp b/stim/cuda/bsds500/dG_d2x_theta_conv2.cpp deleted file mode 100644 index a268d95..0000000 --- a/stim/cuda/bsds500/dG_d2x_theta_conv2.cpp +++ /dev/null @@ -1,52 +0,0 @@ -#include -#include -#include -#include - -/// This function evaluates the theta-dependent even-symmetric gaussian derivative gradient of an one-channel image - -/// @param img is the one-channel image -/// @param r is an array of radii for different scaled discs(filters) -/// @param sigma_n is the number of standard deviations used to define the sigma -/// @param theta is angle used for computing the gradient - -void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M); -void array_abs(float* img, unsigned int N); - -stim::image Gd_even(stim::image image, int r, unsigned int sigma_n, float theta){ - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned N = w * h; // get the number of pixels of picture - int winsize = 2 * r + 1; // set the winsdow size of disc(filter) - float sigma = float(r)/float(sigma_n); // calculate the sigma used in gaussian function - - stim::image I(w, h, 1, 2); // allocate space for return image class - stim::image Gd_even_theta(w, h); // allocate space for Gd_even_theta - stim::image mask_x(winsize, winsize); // allocate space for x-axis-oriented filter - stim::image mask_r(winsize, winsize); // allocate space for theta-oriented filter - - for (int j = 0; j < winsize; j++){ - for (int i = 0; i< winsize; i++){ - - int x = i - r; //range of x - int y = j - r; //range of y - - // create the x-oriented gaussian derivative filter mask_x - 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))); - - } - } - - mask_r = mask_x.rotate(theta, r, r); - //mask_r = mask_x.rotate(45, r, r); - //stim::cpu2image(mask_r.data(), "data_output/mask_r_0919.bmp", winsize, winsize, stim::cmBrewer); - - // do the 2D convolution with image and mask - conv2(image.data(), mask_r.data(), Gd_even_theta.data(), w, h, winsize); - array_abs(Gd_even_theta.data(), N); - - //stim::cpu2image(Gd_even_theta.data(), "data_output/Gd_even_theta_0919.bmp", w, h, stim::cmGrayscale); - - return Gd_even_theta; -} \ No newline at end of file diff --git a/stim/cuda/bsds500/kmeans.cpp b/stim/cuda/bsds500/kmeans.cpp deleted file mode 100644 index 083cae8..0000000 --- a/stim/cuda/bsds500/kmeans.cpp +++ /dev/null @@ -1,67 +0,0 @@ -#include -//#include -#include -#include -#include -#include - -/// This function use cvkmeans to cluster given textons - -/// @param testons is a multi-channel image -/// @param K is the number of clusters - -stim::image kmeans(stim::image textons, unsigned int K){ - - unsigned int w = textons.width(); // get the width of picture - unsigned int h = textons.height(); // get the height of picture - unsigned int feature_n = textons.channels(); // get the spectrum of picture - unsigned int N = w * h; // get the number of pixels - - float* sample1 = (float*) malloc(sizeof(float) * N * feature_n); //allocate the space for textons - - //reallocate a multi-channel texton image to a single-channel image - for(unsigned int c = 0; c < feature_n; c++){ - - stim::image temp; - temp = textons.channel(c); - - for(unsigned int j = 0; j < N; j++){ - - sample1[c + j * feature_n] = temp.data()[j]; - } - } - - - cv::Mat sample2(N, feature_n, CV_32F, sample1); //copy image to cv::mat - - //(optional) show the test result - //imshow("sample2", sample2); - - - cv::TermCriteria criteria(CV_TERMCRIT_EPS+CV_TERMCRIT_ITER, 10, 0.1); // set stop-criteria for kmeans iteration - cv::Mat labels(N, 1, CV_8U, cvScalarAll(0)); // allocate space for kmeans output - cv::Mat centers; // allocate space for kmeans output - - unsigned int test_times = 2; // set the number of times of trying kmeans, it will return the best result - - cv::kmeans(sample2, K, labels, criteria, test_times, cv::KMEANS_PP_CENTERS, centers); // kmeans clustering - - //(optional) show the test result - //imwrite( "data_output/labels_1D.bmp", labels); - - stim::image texture(w, h, 1, 1); // allocate space for texture - - for(unsigned int i = 0; i < N; i++){ // reshape the labels from iD array to image - - texture.data()[i] = labels.at(i); - - } - - //texture.save("data_output/kmeans_test0924_2.bmp"); - - //(optional) show the test result - //stim::cpu2image(texture.data(), "data_output/kmeans_test.bmp", w, h, stim::cmBrewer); - - return texture; - -} \ No newline at end of file diff --git a/stim/cuda/bsds500/laplacian_conv2.cpp b/stim/cuda/bsds500/laplacian_conv2.cpp deleted file mode 100644 index e773e94..0000000 --- a/stim/cuda/bsds500/laplacian_conv2.cpp +++ /dev/null @@ -1,43 +0,0 @@ -#include -#include -#include -#include - -#define PI 3.1415926 -#define SIGMA_N 3 - -void array_multiply(float* lhs, float rhs, unsigned int N); -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N); -void array_abs(float* img, unsigned int N); - -/// This function evaluates the center-surround(Laplacian of Gaussian) gaussian derivative gradient of an one-channel image - -/// @param img is the one-channel image -/// @param r is an array of radii for different scaled discs(filters) -/// @param sigma_n is the number of standard deviations used to define the sigma - -stim::image laplacian_conv2(stim::image image, int sigma){ - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned N = w * h; // get the number of pixels of picture - - int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter - - stim::image I(w, h, 1, 2); // allocate space for return image of dG2_conv2 - stim::image Ixx(w, h); // allocate space for Ixx - stim::image Iyy(w, h); // allocate space for Iyy - stim::image laplacian(w, h); // allocate space for Pb - - I = dG2_conv2(image, sigma); // calculate the Ixx, Iyy - Ixx = I.channel(0); - Iyy = I.channel(1); - - array_add(Ixx.data(), Iyy.data(), laplacian.data(), N); //laplacian = Ixx + Iyy; - array_abs(laplacian.data(), N); - - //stim::cpu2image(laplacian.data(), "data_output/laplacian_0919.bmp", w, h, stim::cmBrewer); - - return laplacian; - -} diff --git a/stim/cuda/bsds500/main.cpp b/stim/cuda/bsds500/main.cpp deleted file mode 100644 index ddea1b6..0000000 --- a/stim/cuda/bsds500/main.cpp +++ /dev/null @@ -1,118 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -/// calculate the cPb, tPb and mPb given a multi-channel image - -void array_multiply(float* lhs, float rhs, unsigned int N); -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N); - -int main() -{ - stim::image image; // generate an image object - - //std::ostringstream ss1,ss2; // (optional) set the stream to designate the test result file - - //for(unsigned int i = 27; i < 31; i++){ - - //ss1 << "3d_sample/b0"<< i << ".bmp"; // (optional) set the name for test result file - //std::string sss1 = ss1.str(); // (optional) - - //image.load("bone1001_3.bmp"); - //image.load(sss1); - //image.load("101085.bmp"); // load the input image - image.load("slice00_500_500.bmp"); // load the input image - image = image.channel(0); // get the first channel of black-and-white image - - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned int N = w * h; // get the number of pixels - int c = image.channels(); // get the number if channels of picture - int s = 3; // set the number of scales - 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 - 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 - float alpha[3] = {1,1,1}; // set an array of weights for different scaled discs(filters) - unsigned int theta_n = 8; // set the number of angles used in filter orientation - unsigned int bin_n = 16; // set the number of bins used in chi-distance - unsigned int K = 16; // set the number of cludters, K should be multiple of bin_n - - stim::image img_cPb(w, h); // allocate the space for cPb - stim::image img_tPb(w, h); // allocate the space for tPb - stim::image img_mPb(w, h); // allocate the space for tPb - - std::cout<<"imagesize: "<< w <<"*"<< h <<'\n'; - - //*******************cPb******************** - std::clock_t start1; // (optional) set the timer to calculate the total time - start1 = std::clock(); // (optional) set timer start point - - std::cout<<"begin cPb"<<'\n'; - - img_cPb = cPb(image, sigma, alpha, s); - - // show the cPb (optional) - stim::cpu2image(img_cPb.data(), "data_output/img_cPb_1008.bmp", w, h, stim::cmBrewer); - //ss2 << "3d_sample/cPb0"<< i << ".bmp"; // (optional) set the name for test result file - //std::string sss2 = ss2.str(); - - //stim::cpu2image(img_cPb.data(), "data_output/bone_cPb3.1.bmp", w, h, stim::cmBrewer); - - double duration1 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time - std::cout<<"cPb time: "<< duration1 <<"s"<<'\n'; // (optional) show the total time - - //ss1.str(""); //(optional) clear the space for stream - //ss2.str(""); //(optional) clear the space for stream - - - //*******************tPb******************** - - std::cout<<"begin tPb"<<'\n'; - - - std::clock_t start2; // (optional) set the timer to calculate the total time - start2 = std::clock(); // (optional) set timer start point - - img_tPb = tPb(image, r, alpha, theta_n, bin_n, s, K); - - // show the tPb (optional) - // stim::cpu2image(img_tPb.data(), "data_output/bone_tPb_3.1.bmp", w, h, stim::cmBrewer); - stim::cpu2image(img_tPb.data(), "data_output/img_tPb_1008.bmp", w, h, stim::cmBrewer); - - double duration2 = ( std::clock() - start2 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time - std::cout<<"tPb time: "<< duration2 <<"s"<<'\n'; // (optional) show the total time - - - //******************************************* - - double duration3 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time - std::cout<<"total running time: "<< duration3 <<"s"<<'\n'; // (optional) show the total time - //} - //******************mPb********************** - // set parameters for linear combination - - float a = 1; - float b = 0.5; - - // create mPb by linearly combined cPb and tPb - array_multiply(img_cPb.data(), a, N); - array_multiply(img_tPb.data(), b, N); - array_add(img_cPb.data(), img_tPb.data(), img_mPb.data(), N); - - //ss2 << "3d_sample/mPb0"<< i << ".bmp"; // (optional) set the name for test result file - //std::string sss2 = ss2.str(); - // show the mPb (optional) - stim::cpu2image(img_mPb.data(), "data_output/img_mPb_1008.bmp", w, h, stim::cmBrewer); - - //stim::cpu2image(img_mPb.data(), sss2, w, h, stim::cmBrewer); - - //ss1.str(""); //(optional) clear the space for stream - //ss2.str(""); //(optional) clear the space for stream - - //} - return 0; - -} diff --git a/stim/cuda/bsds500/tPb.cpp b/stim/cuda/bsds500/tPb.cpp deleted file mode 100644 index bc1a44a..0000000 --- a/stim/cuda/bsds500/tPb.cpp +++ /dev/null @@ -1,97 +0,0 @@ -#include -#include -#include -#include -#include - - -void array_multiply(float* lhs, float rhs, unsigned int N); -void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N); -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); - -/// This function evaluates the tPb given a grayscale image - -/// @param img is the multi-channel image -/// @param theta_n is the number of angles used for computing oriented chi-gradient -/// @param r is an array of radii for different scaled discs(filters) -/// @param alpha is is an array of weights for different scaled discs(filters) -/// @param s is the number of scales -/// @param K is the number of clusters - -stim::image tPb(stim::image img, int* r, float* alpha, unsigned int theta_n, unsigned int bin_n, int s, unsigned K){ - - unsigned int w = img.width(); // get the width of picture - unsigned int h = img.height(); // get the height of picture - unsigned int N = w * h; // get the number of pixels - - stim::image img_textons(w, h, 1, theta_n*2+1); // allocate space for img_textons - stim::image img_texture(w, h, 1, 1); // allocate space for img_texture - stim::image tPb_theta(w, h, 1, 1); // allocate space for tPb_theta - stim::image tPb(w, h, 1, 1); // allocate space for tPb - unsigned size = tPb_theta.size(); // get the size of tPb_theta - memset (tPb.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb to 0 - stim::image temp(w, h, 1, 1); // set the temporary image to store the addtion result - - std::ostringstream ss; // (optional) set the stream to designate the test result file - - - img_textons = textons(img, theta_n); - - img_texture = kmeans(img_textons, K); // changing kmeans result into float type is required - - stim::cpu2image(img_texture.data(), "data_output/texture.bmp", w, h, stim::cmBrewer); - - - unsigned int max1 = img_texture.maxv(); // get the maximum of Pb used for normalization - unsigned int bin_size = (max1 + 1)/bin_n; // (whether"+1" or not depends on kmeans result) - - for (int i = 0; i < theta_n; i++){ - - float theta = 180 * ((float)i/theta_n); // calculate the even-splited angle for each tPb_theta - - memset (tPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb_theta to 0 - - //ss << "data_output/0922tPb_theta"<< theta << ".bmp"; // set the name for test result file (optional) - //std::string sss = ss.str(); - - for (int j = 0; j < s; j++){ - - // get the chi-gradient by convolving each image slice with the mask - chi_grad(img_texture.data(), temp.data(), w, h, r[j], bin_n, bin_size, theta); - - float max2 = temp.maxv(); // get the maximum of tPb_theta used for normalization - array_multiply(temp.data(), 1/max2, N); // normalize the tPb_theta - - //output the test result of each slice (optional) - //stim::cpu2image(temp.data(), "data_output/tPb_slice0924_2.bmp", w, h, stim::cmBrewer); - - // multiply each chi-gradient with its weight - array_multiply(temp.data(), alpha[j], N); - - // add up all the weighted chi-gradients - array_add(tPb_theta.data(), temp.data(), tPb_theta.data(), N); - - - } - - //ss.str(""); //(optional) clear the space for stream - - for(unsigned long ti = 0; ti < N; ti++){ - - if(tPb_theta.data()[ti] > tPb.data()[ti]){ //get the maximum value among all tPb_theta for ith pixel - tPb.data()[ti] = tPb_theta.data()[ti]; - } - - else{ - } - } - } - - float max3 = tPb.maxv(); // get the maximum of tPb used for normalization - array_multiply(tPb.data(), 1/max3, N); // normalize the tPb - - //output the test result of tPb (optional) - //stim::cpu2image(tPb.data(), "data_output/tPb_0922.bmp", w, h, stim::cmBrewer); - - return tPb; -} diff --git a/stim/cuda/bsds500/textons.cpp b/stim/cuda/bsds500/textons.cpp deleted file mode 100644 index 484b27d..0000000 --- a/stim/cuda/bsds500/textons.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include -//#include -#include -#include -#include - -/// This function convolve the grayscale image with a set of oriented Gaussian -/// derivative filters, and return a texton image with (theta_n*2+1) channels - -/// @param image is an one-channel grayscale image -/// @param theta_n is the number of angles used for computing the gradient - -stim::image textons(stim::image image, unsigned int theta_n){ - - unsigned int w = image.width(); // get the width of picture - unsigned int h = image.height(); // get the height of picture - unsigned N = w * h; // get the number of pixels of picture - - stim::image textons(w, h, 1, theta_n*2+1); // allocate space for textons - stim::image temp(w, h); // allocate space for temp - - int sigma = 1; // set sigma for odd-symmetric, even-symmetric and center-surround filter filter - - //std::ostringstream ss1, ss2; // (optional) set the stream to designate the test result file - - for (unsigned int i = 0; i < theta_n; i++){ - - //ss1 << "data_output/textons_channel_"<< i << ".bmp"; // set the name for test result file (optional) - //std::string sss1 = ss1.str(); - //ss2 << "data_output/textons_channel_"<< i+theta_n << ".bmp"; // set the name for test result file (optional) - //std::string sss2 = ss2.str(); - - float theta = 180 * ((float)i/theta_n); // calculate the even-splited angle for each oriented filter - - temp = dG1_theta_conv2(image, sigma, theta); // return dG1_theta_conv2 to temp - //stim::cpu2image(temp.data(), sss1, w, h, stim::cmBrewer); - textons.set_channel(i, temp.data()); // copy temp to ith channel of textons - - temp = dG2_d2x_theta_conv2(image, sigma, theta); // return dG2_d2x_theta_conv2 to temp - //stim::cpu2image(temp.data(), sss2, w, h, stim::cmBrewer); - textons.set_channel(i + theta_n, temp.data()); // copy temp to (i+theta_n)th channel of textons - - //ss1.str(""); //(optional) clear the space for stream - //ss2.str(""); //(optional) clear the space for stream - - } - - temp = laplacian_conv2(image, sigma); // return laplacian_conv2 to temp - //stim::cpu2image(temp.data(), "data_output/textons_channel_16.bmp", w, h, stim::cmBrewer); - textons.set_channel(theta_n*2, temp.data()); // copy temp to (theta_n*2)th channel of textons - - return textons; - -} - - \ No newline at end of file diff --git a/stim/image/image_contour_detection.h b/stim/image/image_contour_detection.h deleted file mode 100644 index 60e6d03..0000000 --- a/stim/image/image_contour_detection.h +++ /dev/null @@ -1,16 +0,0 @@ - -//stim::image gaussian_derivative_filter_odd(stim::image image, int r, unsigned int sigma_n, float theta); -//stim::image func_mPb_theta(stim::image img, float theta, int* r, float* alpha, int s); -//stim::image func_mPb(stim::image img, unsigned int theta_n, int* r, float* alpha, int s); - -stim::image dG1_conv2(stim::image image, int sigma); -stim::image dG2_conv2(stim::image image, int sigma); -stim::image dG1_theta_conv2(stim::image image, int sigma, float theta); -stim::image dG2_d2x_theta_conv2(stim::image image, int sigma, float theta); -stim::image laplacian_conv2(stim::image image, int sigma); - -stim::image textons(stim::image image, unsigned int theta_n); -stim::image kmeans(stim::image textons, unsigned int K); -stim::image Pb(stim::image image, int sigma); -stim::image cPb(stim::image img, int* sigma, float* alpha, int s); -stim::image tPb(stim::image img, int* r, float* alpha, unsigned int theta_n, unsigned int bin_n, int s, unsigned int K); \ No newline at end of file -- libgit2 0.21.4