codebase / stimlib

Commit 018b00d675350835a968dba69361027a39b9f3b7

Authored by David Mayerich
1 parent 85025dd1

adapted bsds500 files to the stim library, tested cPb

Showing 12 changed files with 805 additions and 0 deletions   Show diff stats
Inline Side-by-side
stim/cuda/arraymath/array_atan2.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_ARRAY_ATAN2_H
  2 +#define STIM_CUDA_ARRAY_ATAN2_H
  3 +
  4 +#include <iostream>
  5 +#include <cuda.h>
  6 +#include <cmath>
  7 +#include <stim/cuda/cudatools.h>
  8 +
  9 +namespace stim{
  10 + namespace cuda{
  11 +
  12 + template<typename T>
  13 + __global__ void cuda_atan2(T* y, T* x, T* r, unsigned int N){
  14 +
  15 + //calculate the 1D index for this thread
  16 + int idx = blockIdx.x * blockDim.x + threadIdx.x;
  17 +
  18 + if(idx < N){
  19 + r[idx] = atan2(y[idx], x[idx]);
  20 + }
  21 +
  22 + }
  23 +
  24 + template<typename T>
  25 + void gpu_atan2(T* y, T* x, T* r, unsigned int N){
  26 +
  27 + //get the maximum number of threads per block for the CUDA device
  28 + int threads = stim::maxThreadsPerBlock();
  29 +
  30 + //calculate the number of blocks
  31 + int blocks = N / threads + 1;
  32 +
  33 + //call the kernel to do the multiplication
  34 + cuda_atan2 <<< blocks, threads >>>(y, x, r, N);
  35 +
  36 + }
  37 +
  38 + template<typename T>
  39 + void cpu_atan2(T* y, T* x, T* cpu_r, unsigned int N){
  40 +
  41 + //allocate memory on the GPU for the array
  42 + T* gpu_x;
  43 + T* gpu_y;
  44 + T* gpu_r;
  45 + HANDLE_ERROR( cudaMalloc( &gpu_x, N * sizeof(T) ) );
  46 + HANDLE_ERROR( cudaMalloc( &gpu_y, N * sizeof(T) ) );
  47 + HANDLE_ERROR( cudaMalloc( &gpu_r, N * sizeof(T) ) );
  48 +
  49 + //copy the array to the GPU
  50 + HANDLE_ERROR( cudaMemcpy( gpu_x, x, N * sizeof(T), cudaMemcpyHostToDevice) );
  51 + HANDLE_ERROR( cudaMemcpy( gpu_y, y, N * sizeof(T), cudaMemcpyHostToDevice) );
  52 +
  53 + //call the GPU version of this function
  54 + gpu_atan2<T>(gpu_y, gpu_x ,gpu_r, N);
  55 +
  56 + //copy the array back to the CPU
  57 + HANDLE_ERROR( cudaMemcpy( cpu_r, gpu_r, N * sizeof(T), cudaMemcpyDeviceToHost) );
  58 +
  59 + //free allocated memory
  60 + cudaFree(gpu_x);
  61 + cudaFree(gpu_y);
  62 + cudaFree(gpu_r);
  63 +
  64 + }
  65 +
  66 + }
  67 +}
  68 +
  69 +
  70 +
  71 +#endif
0 \ No newline at end of file 72 \ No newline at end of file
stim/cuda/bsds500.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
stim/cuda/bsds500/Pb.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_PB_CUH
  2 +#define STIM_CUDA_PB_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +//#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +
  9 +#include "dG1_conv2.cuh"
  10 +#include <stim/cuda/arraymath/array_abs.cuh>
  11 +#include <stim/cuda/arraymath/array_divide.cuh>
  12 +#include <stim/cuda/arraymath/array_atan.cuh>
  13 +#include <stim/cuda/arraymath/array_atan2.cuh>
  14 +#include <stim/cuda/arraymath/array_cos.cuh>
  15 +#include <stim/cuda/arraymath/array_sin.cuh>
  16 +#include <stim/cuda/arraymath/array_multiply.cuh>
  17 +#include <stim/cuda/arraymath/array_add.cuh>
  18 +
  19 +
  20 +
  21 +
  22 +
  23 +
  24 +#define SIGMA_N 3
  25 +
  26 +//void array_abs(float* img, unsigned int N);
  27 +//void array_multiply(float* lhs, float rhs, unsigned int N);
  28 +//void array_cos(float* ptr1, float* cpu_out, unsigned int N);
  29 +//void array_sin(float* ptr1, float* cpu_out, unsigned int N);
  30 +//void array_atan(float* ptr1, float* cpu_out, unsigned int N);
  31 +//void array_divide(float* ptr1, float* ptr2,float* cpu_quotient, unsigned int N);
  32 +//void array_multiply(float* ptr1, float* ptr2, float* product, unsigned int N);
  33 +//void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
  34 +
  35 +/// This function uses odd-symmetric gaussian derivative filter to evaluate
  36 +/// the max probability of a contour on one scale, given an one-channel image
  37 +
  38 +/// @param img is an one-channel image
  39 +/// @param r is an array of radii for different scaled discs(filters)
  40 +/// @param sigma_n is the number of standard deviations used to define the sigma
  41 +
  42 +stim::image<float> Pb(stim::image<float> image, int sigma){
  43 +
  44 + unsigned int w = image.width(); // get the width of picture
  45 + unsigned int h = image.height(); // get the height of picture
  46 + unsigned N = w * h; // get the number of pixels of picture
  47 +
  48 + stim::image<float> I(w, h, 1, 2); // allocate space for return image of dG1_conv2
  49 + stim::image<float> theta(w, h); // allocate space for theta matrix
  50 + stim::image<float> cos(w, h); // allocate space for cos(theta)
  51 + stim::image<float> sin(w, h); // allocate space for sin(theta)
  52 + stim::image<float> temp(w, h); // allocate space for temp
  53 + stim::image<float> Ix(w, h); // allocate space for Ix
  54 + stim::image<float> Iy(w, h); // allocate space for Iy
  55 + stim::image<float> Pb(w, h); // allocate space for Pb
  56 +
  57 + I = dG1_conv2(image, sigma); // calculate the Ix, Iy
  58 + Ix = I.channel(0);
  59 + stim::cuda::cpu_abs(Ix.data(), N); //get |Ix|;
  60 + //stim::cpu2image(Ix.data(), "data_output/Pb_Ix_0924.bmp", w, h, stim::cmBrewer);
  61 + Iy = I.channel(1);
  62 + stim::cuda::cpu_abs(Iy.data(), N); //get |Iy|;
  63 + //stim::cpu2image(Iy.data(), "data_output/Pb_Iy_0924.bmp", w, h, stim::cmBrewer);
  64 +
  65 + //stim::cuda::cpu_divide(Iy.data(), Ix.data(), theta.data(), N); //temp = Iy./Ix
  66 + stim::cuda::cpu_atan2(Iy.data(), Ix.data(), theta.data(), N); //temp = Iy./Ix
  67 + //stim::cuda::cpu_atan(temp.data(), theta.data(), N); //theta = atan(temp)
  68 + stim::cuda::cpu_cos(theta.data(), cos.data(), N); //cos = cos(theta)
  69 + stim::cuda::cpu_sin(theta.data(), sin.data(), N); //sin = sin(theta)
  70 + stim::cuda::cpu_multiply(Ix.data(), cos.data(), Ix.data(), N); //Ix = Ix.*cos
  71 + stim::cuda::cpu_multiply(Iy.data(), sin.data(), Iy.data(), N); //Iy = Iy.*sin
  72 + stim::cuda::cpu_add(Ix.data(), Iy.data(), Pb.data(), N); //Pb = Ix + Iy;
  73 +
  74 + float max = Pb.maxv(); // get the maximum of Pb used for normalization
  75 + stim::cuda::cpu_multiply(Pb.data(), 1/max, N); // normalize the Pb
  76 +
  77 + //stim::cpu2image(Pb.data(), "data_output/Pb_0924.bmp", w, h, stim::cmBrewer); show the Pb(optional)
  78 +
  79 + return Pb;
  80 +
  81 +}
  82 +
  83 +#endif
stim/cuda/bsds500/cPb.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_CPB_CUH
  2 +#define STIM_CUDA_CPB_CUH
  3 +
  4 +#include <cmath>
  5 +#include <sstream>
  6 +
  7 +#include <stim/image/image.h>
  8 +#include <stim/visualization/colormap.h>
  9 +#include <stim/cuda/arraymath/array_multiply.cuh>
  10 +#include <stim/cuda/arraymath/array_add.cuh>
  11 +
  12 +//BSDS500 files
  13 +#include "Pb.cuh"
  14 +
  15 +//void array_multiply(float* lhs, float rhs, unsigned int N);
  16 +//void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
  17 +
  18 +/// This function evaluates the cPb given an multi-channel image
  19 +
  20 +/// @param img is the multi-channel image
  21 +/// @param r is an array of radii for different scaled discs(filters)
  22 +/// @param alpha is is an array of weights for different scaled discs(filters)
  23 +/// @param s is the number of scales
  24 +
  25 +stim::image<float> cPb(stim::image<float> img, int* sigma, float* alpha, int s){
  26 +
  27 + unsigned int w = img.width(); // get the width of picture
  28 + unsigned int h = img.height(); // get the height of picture
  29 + unsigned int c = img.channels(); // get the channels of picture
  30 +
  31 +
  32 + stim::image<float> cPb(w, h, 1); // allocate space for cPb
  33 + unsigned size = cPb.size(); // get the size of cPb
  34 + memset ( cPb.data(), 0, size * sizeof(float)); // initialize all the pixels of cPb to 0
  35 +
  36 +
  37 + unsigned int N = w * h; // get the number of pixels
  38 +
  39 + //std::ostringstream ss; // (optional) set the stream to designate the test result file
  40 +
  41 + stim::image<float> temp; // set the temporary image to store the addtion result
  42 +
  43 + for (int i = 0; i < c; i++){
  44 + for (int j = 0; j < s; j++){
  45 +
  46 + //ss << "data_output/cPb_slice"<< i*s + j << ".bmp"; // set the name for test result file (optional)
  47 + //std::string sss = ss.str();
  48 +
  49 + // get the gaussian gradient by convolving each image slice with the mask
  50 + temp = Pb(img.channel(i), sigma[i*s + j]);
  51 +
  52 + // output the test result of each slice (optional)
  53 + //stim::cpu2image(temp.data(), sss, w, h, stim::cmBrewer);
  54 +
  55 + // multiply each gaussian gradient with its weight
  56 + stim::cuda::cpu_multiply(temp.data(), alpha[i*s + j], N);
  57 +
  58 + // add up all the weighted gaussian gradients
  59 + stim::cuda::cpu_add(cPb.data(), temp.data(), cPb.data(), N);
  60 +
  61 + //ss.str(""); //(optional) clear the space for stream
  62 +
  63 + }
  64 + }
  65 +
  66 + float max = cPb.maxv(); // get the maximum of cPb used for normalization
  67 + stim::cuda::cpu_multiply(cPb.data(), 1/max, N); // normalize the cPb
  68 +
  69 + // output the test result of cPb (optional)
  70 + //stim::cpu2image(cPb.data(), "data_output/cPb_0916.bmp", w, h, stim::cmBrewer);
  71 +
  72 + return cPb;
  73 +}
  74 +
  75 +#endif
stim/cuda/bsds500/dG1_conv2.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_DG1_CONV2_CUH
  2 +#define STIM_CUDA_DG1_CONV2_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +//#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +
  9 +#include <stim/cuda/templates/conv2sep.cuh>
  10 +
  11 +#define SIGMA_N 3
  12 +
  13 +/// This function generates the first-order gaussian derivative filter gx gy,
  14 +/// convolves the image with gx gy,
  15 +/// and returns an image class which channel(0) is Ix and channel(1) is Iy
  16 +
  17 +/// @param img is the one-channel image
  18 +/// @param sigma is the parameter for gaussian function
  19 +
  20 +//void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1);
  21 +//void array_abs(float* img, unsigned int N);
  22 +
  23 +stim::image<float> dG1_conv2(stim::image<float> image, int sigma){
  24 +
  25 + unsigned int w = image.width(); // get the width of picture
  26 + unsigned int h = image.height(); // get the height of picture
  27 +
  28 + int r = SIGMA_N * sigma;
  29 + int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
  30 +
  31 + stim::image<float> I(w, h, 1, 2); // allocate space for return image class
  32 + stim::image<float> Ix(w, h); // allocate space for Ix
  33 + stim::image<float> Iy(w, h); // allocate space for Iy
  34 + Ix = image; // initialize Ix
  35 + Iy = image; // initialize Iy
  36 +
  37 + float* array_x1;
  38 + array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gx
  39 + float* array_y1;
  40 + array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gx
  41 + float* array_x2;
  42 + array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gy
  43 + float* array_y2;
  44 + array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gy
  45 +
  46 +
  47 + for (int i = 0; i < winsize; i++){
  48 +
  49 + int x = i - r; //range of x
  50 + int y = i - r; //range of y
  51 +
  52 + // create the 1D x-oriented gaussian derivative filter array_x1 for gx
  53 + array_x1[i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  54 + // create the 1D y-oriented gaussian derivative filter array_y1 for gx
  55 + array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  56 + // create the 1D x-oriented gaussian derivative filter array_x2 for gy
  57 + array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  58 + // create the 1D y-oriented gaussian derivative filter array_y2 for gy
  59 + array_y2[i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  60 + }
  61 +
  62 + //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  63 + //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  64 + //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  65 + //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  66 +
  67 + // get Ix by convolving the image with gx
  68 + //conv2_sep(Ix.data(), w, h, array_x1, winsize, array_y1, winsize);
  69 + stim::cuda::cpu_conv2sep(Ix.data(), w, h, array_x1, winsize, array_y1, winsize);
  70 +
  71 + //stim::cpu2image(Ix.data(), "data_output/Ix_0915.bmp", w, h, stim::cmBrewer);
  72 + // get Iy by convolving the image with gy
  73 + //conv2_sep(Iy.data(), w, h, array_x2, winsize, array_y2, winsize);
  74 + stim::cuda::cpu_conv2sep(Iy.data(), w, h, array_x2, winsize, array_y2, winsize);
  75 +
  76 + //stim::cpu2image(Iy.data(), "data_output/Iy_0915.bmp", w, h, stim::cmBrewer);
  77 +
  78 + delete [] array_x1; //free the memory of array_x1
  79 + delete [] array_y1; //free the memory of array_y1
  80 + delete [] array_x2; //free the memory of array_x2
  81 + delete [] array_y2; //free the memory of array_y2
  82 +
  83 + I.set_channel(0, Ix.data());
  84 + I.set_channel(1, Iy.data());
  85 +
  86 + return I;
  87 +
  88 +}
  89 +
  90 +#endif
0 \ No newline at end of file 91 \ No newline at end of file
stim/cuda/bsds500/dG1_theta_conv2.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_DG1_THETA_CONV2_CUH
  2 +#define STIM_CUDA_DG1_THETA_CONV2_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +
  9 +#define PI 3.1415926
  10 +#define SIGMA_N 3
  11 +
  12 +//void array_multiply(float* lhs, float rhs, unsigned int N);
  13 +//void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
  14 +//void array_abs(float* img, unsigned int N);
  15 +
  16 +/// This function evaluates the theta-dependent odd symmetric gaussian derivative gradient of an one-channel image
  17 +
  18 +/// @param img is the 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 +/// @param theta is angle used for computing the gradient
  22 +
  23 +stim::image<float> dG1_theta_conv2(stim::image<float> image, int sigma, float theta){
  24 +
  25 + float theta_r = (theta * PI)/180; //change angle unit from degree to rad
  26 +
  27 + unsigned int w = image.width(); // get the width of picture
  28 + unsigned int h = image.height(); // get the height of picture
  29 +
  30 + int r = SIGMA_N * sigma;
  31 + unsigned N = w * h; // get the number of pixels of picture
  32 +
  33 + int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
  34 +
  35 + stim::image<float> I(w, h, 1, 2); // allocate space for return image of dG1_conv2
  36 + stim::image<float> Ix(w, h); // allocate space for Ix
  37 + stim::image<float> Iy(w, h); // allocate space for Iy
  38 + stim::image<float> dG1_theta(w, h); // allocate space for Pb
  39 +
  40 + I = dG1_conv2(image, sigma); // calculate the Ix, Iy
  41 + Ix = I.channel(0);
  42 + Iy = I.channel(1);
  43 +
  44 + array_multiply(Ix.data(), cos(theta_r), N); //Ix = Ix*cos(theta_r)
  45 + array_multiply(Iy.data(), sin(theta_r), N); //Iy = Iy*sin(theta_r)
  46 + array_add(Ix.data(), Iy.data(), dG1_theta.data(), N); //dG1_theta = Ix + Iy;
  47 + array_abs(dG1_theta.data(), N);
  48 +
  49 + //stim::cpu2image(I.channel(0).data(), "data_output/dG1_theta_x_0919.bmp", w, h, stim::cmBrewer);
  50 + //stim::cpu2image(I.channel(1).data(), "data_output/dG1_theta_y_0919.bmp", w, h, stim::cmBrewer);
  51 + //stim::cpu2image(dG1_theta.data(), "data_output/dG1_theta_0919.bmp", w, h, stim::cmBrewer);
  52 +
  53 + return dG1_theta;
  54 +
  55 +}
  56 +
  57 +#endif
0 \ No newline at end of file 58 \ No newline at end of file
stim/cuda/bsds500/dG2_conv2.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_DG2_CONV2_CUH
  2 +#define STIM_CUDA_DG2_CONV2_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +//#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +#define SIGMA_N 3
  9 +
  10 +/// This function generates the second-order gaussian derivative filter gxx gyy,
  11 +/// convolves the image with gxx gyy,
  12 +/// and returns an image class which channel(0) is Ixx and channel(1) is Iyy
  13 +
  14 +/// @param img is the one-channel image
  15 +/// @param sigma is the parameter for gaussian function
  16 +
  17 +//void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1);
  18 +//void array_abs(float* img, unsigned int N);
  19 +
  20 +stim::image<float> dG2_conv2(stim::image<float> image, int sigma){
  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 N = w * h; // get the number of pixels of picture
  25 +
  26 + int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
  27 + int r = SIGMA_N * sigma;
  28 +
  29 + stim::image<float> I(w, h, 1, 2); // allocate space for return image class
  30 + stim::image<float> Ixx(w, h); // allocate space for Ixx
  31 + stim::image<float> Iyy(w, h); // allocate space for Iyy
  32 + Ixx = image; // initialize Ixx
  33 + Iyy = image; // initialize Iyy
  34 +
  35 + float* array_x1;
  36 + array_x1 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gxx
  37 + float* array_y1;
  38 + array_y1 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gxx
  39 + float* array_x2;
  40 + array_x2 = new float[winsize]; //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gyy
  41 + float* array_y2;
  42 + array_y2 = new float[winsize]; //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gyy
  43 +
  44 +
  45 + for (int i = 0; i < winsize; i++){
  46 +
  47 + int x = i - r; //range of x
  48 + int y = i - r; //range of y
  49 +
  50 + // create the 1D x-oriented gaussian derivative filter array_x1 for gxx
  51 + array_x1[i] = (-1) * (1 - pow(x, 2)) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  52 + // create the 1D y-oriented gaussian derivative filter array_y1 for gxx
  53 + array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  54 + // create the 1D x-oriented gaussian derivative filter array_x2 for gyy
  55 + array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
  56 + // create the 1D y-oriented gaussian derivative filter array_y2 for gyy
  57 + array_y2[i] = (-1) * (1 - pow(y, 2)) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
  58 + }
  59 +
  60 + //stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  61 + //stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  62 + //stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  63 + //stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
  64 +
  65 + // get Ixx by convolving the image with gxx
  66 + conv2_sep(Ixx.data(), w, h, array_x1, winsize, array_y1, winsize);
  67 +
  68 + //stim::cpu2image(Ixx.data(), "data_output/Ixx_0915.bmp", w, h, stim::cmBrewer);
  69 + // get Iyy by convolving the image with gyy
  70 + conv2_sep(Iyy.data(), w, h, array_x2, winsize, array_y2, winsize);
  71 +
  72 + //stim::cpu2image(Iyy.data(), "data_output/Iyy_0915.bmp", w, h, stim::cmBrewer);
  73 +
  74 + delete [] array_x1; //free the memory of array_x1
  75 + delete [] array_y1; //free the memory of array_y1
  76 + delete [] array_x2; //free the memory of array_x2
  77 + delete [] array_y2; //free the memory of array_y2
  78 +
  79 + I.set_channel(0, Ixx.data());
  80 + I.set_channel(1, Iyy.data());
  81 +
  82 + return I;
  83 +
  84 +}
  85 +
  86 +#endif
0 \ No newline at end of file 87 \ No newline at end of file
stim/cuda/bsds500/dG2_d2x_theta_conv2.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_DG2_D2X_THETA_CONV2_CUH
  2 +#define STIM_CUDA_DG2_D2X_THETA_CONV2_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +
  9 +#define SIGMA_N 3
  10 +
  11 +/// This function evaluates the theta-dependent even-symmetric gaussian derivative gradient of an one-channel image
  12 +
  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 +
  18 +//void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M);
  19 +//void array_abs(float* img, unsigned int N);
  20 +
  21 +stim::image<float> dG2_d2x_theta_conv2(stim::image<float> image, int sigma, float theta){
  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 +
  27 + int r = SIGMA_N * sigma; // set the radius of filter
  28 + int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
  29 +
  30 + stim::image<float> I(w, h, 1, 2); // allocate space for return image class
  31 + stim::image<float> dG2_d2x_theta(w, h); // allocate space for dG2_d2x_theta
  32 + stim::image<float> mask_x(winsize, winsize); // allocate space for x-axis-oriented filter
  33 + stim::image<float> mask_r(winsize, winsize); // allocate space for theta-oriented filter
  34 +
  35 + for (int j = 0; j < winsize; j++){
  36 + for (int i = 0; i< winsize; i++){
  37 +
  38 + int x = i - r; //range of x
  39 + int y = j - r; //range of y
  40 +
  41 + // create the x-oriented gaussian derivative filter mask_x
  42 + 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)));
  43 +
  44 + }
  45 + }
  46 +
  47 + mask_r = mask_x.rotate(theta, r, r);
  48 + //mask_r = mask_x.rotate(45, r, r);
  49 + //stim::cpu2image(mask_r.data(), "data_output/mask_r_0919.bmp", winsize, winsize, stim::cmBrewer);
  50 +
  51 + // do the 2D convolution with image and mask
  52 + conv2(image.data(), mask_r.data(), dG2_d2x_theta.data(), w, h, winsize);
  53 + array_abs(dG2_d2x_theta.data(), N);
  54 +
  55 + //stim::cpu2image(dG2_d2x_theta.data(), "data_output/dG2_d2x_theta_0919.bmp", w, h, stim::cmGrayscale);
  56 +
  57 + return dG2_d2x_theta;
  58 +}
  59 +
  60 +#endif
0 \ No newline at end of file 61 \ No newline at end of file
stim/cuda/bsds500/kmeans.h 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_KMEANS_CUH
  2 +#define STIM_CUDA_KMEANS_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +//#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +#include <opencv2/opencv.hpp>
  9 +#include <iostream>
  10 +
  11 +/// This function use cvkmeans to cluster given textons
  12 +
  13 +/// @param testons is a multi-channel image
  14 +/// @param K is the number of clusters
  15 +
  16 +stim::image<float> kmeans(stim::image<float> textons, unsigned int K){
  17 +
  18 + unsigned int w = textons.width(); // get the width of picture
  19 + unsigned int h = textons.height(); // get the height of picture
  20 + unsigned int feature_n = textons.channels(); // get the spectrum of picture
  21 + unsigned int N = w * h; // get the number of pixels
  22 +
  23 + float* sample1 = (float*) malloc(sizeof(float) * N * feature_n); //allocate the space for textons
  24 +
  25 + //reallocate a multi-channel texton image to a single-channel image
  26 + for(unsigned int c = 0; c < feature_n; c++){
  27 +
  28 + stim::image<float> temp;
  29 + temp = textons.channel(c);
  30 +
  31 + for(unsigned int j = 0; j < N; j++){
  32 +
  33 + sample1[c + j * feature_n] = temp.data()[j];
  34 + }
  35 + }
  36 +
  37 +
  38 + cv::Mat sample2(N, feature_n, CV_32F, sample1); //copy image to cv::mat
  39 +
  40 + //(optional) show the test result
  41 + //imshow("sample2", sample2);
  42 +
  43 +
  44 + cv::TermCriteria criteria(CV_TERMCRIT_EPS+CV_TERMCRIT_ITER, 10, 0.1); // set stop-criteria for kmeans iteration
  45 + cv::Mat labels(N, 1, CV_8U, cvScalarAll(0)); // allocate space for kmeans output
  46 + cv::Mat centers; // allocate space for kmeans output
  47 +
  48 + unsigned int test_times = 2; // set the number of times of trying kmeans, it will return the best result
  49 +
  50 + cv::kmeans(sample2, K, labels, criteria, test_times, cv::KMEANS_PP_CENTERS, centers); // kmeans clustering
  51 +
  52 + //(optional) show the test result
  53 + //imwrite( "data_output/labels_1D.bmp", labels);
  54 +
  55 + stim::image<float> texture(w, h, 1, 1); // allocate space for texture
  56 +
  57 + for(unsigned int i = 0; i < N; i++){ // reshape the labels from iD array to image
  58 +
  59 + texture.data()[i] = labels.at<int>(i);
  60 +
  61 + }
  62 +
  63 + //texture.save("data_output/kmeans_test0924_2.bmp");
  64 +
  65 + //(optional) show the test result
  66 + //stim::cpu2image(texture.data(), "data_output/kmeans_test.bmp", w, h, stim::cmBrewer);
  67 +
  68 + return texture;
  69 +
  70 +}
  71 +
  72 +#endif
0 \ No newline at end of file 73 \ No newline at end of file
stim/cuda/bsds500/laplacian_conv2.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_LAPLACIAN_CONV2_CUH
  2 +#define STIM_CUDA_LAPLACIAN_CONV2_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +
  9 +#define PI 3.1415926
  10 +#define SIGMA_N 3
  11 +
  12 +//void array_multiply(float* lhs, float rhs, unsigned int N);
  13 +//void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
  14 +//void array_abs(float* img, unsigned int N);
  15 +
  16 +/// This function evaluates the center-surround(Laplacian of Gaussian) gaussian derivative gradient of an one-channel image
  17 +
  18 +/// @param img is the 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> laplacian_conv2(stim::image<float> image, int sigma){
  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 +
  28 + int winsize = 2 * SIGMA_N * sigma + 1; // set the winsdow size of filter
  29 +
  30 + stim::image<float> I(w, h, 1, 2); // allocate space for return image of dG2_conv2
  31 + stim::image<float> Ixx(w, h); // allocate space for Ixx
  32 + stim::image<float> Iyy(w, h); // allocate space for Iyy
  33 + stim::image<float> laplacian(w, h); // allocate space for Pb
  34 +
  35 + I = dG2_conv2(image, sigma); // calculate the Ixx, Iyy
  36 + Ixx = I.channel(0);
  37 + Iyy = I.channel(1);
  38 +
  39 + array_add(Ixx.data(), Iyy.data(), laplacian.data(), N); //laplacian = Ixx + Iyy;
  40 + array_abs(laplacian.data(), N);
  41 +
  42 + //stim::cpu2image(laplacian.data(), "data_output/laplacian_0919.bmp", w, h, stim::cmBrewer);
  43 +
  44 + return laplacian;
  45 +
  46 +}
  47 +
  48 +#endif
stim/cuda/bsds500/tPb.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_TPB_CUH
  2 +#define STIM_CUDA_TPB_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +#include <sstream>
  9 +
  10 +
  11 +//void array_multiply(float* lhs, float rhs, unsigned int N);
  12 +//void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
  13 +//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);
  14 +
  15 +/// This function evaluates the tPb given a grayscale image
  16 +
  17 +/// @param img is the multi-channel image
  18 +/// @param theta_n is the number of angles used for computing oriented chi-gradient
  19 +/// @param r is an array of radii for different scaled discs(filters)
  20 +/// @param alpha is is an array of weights for different scaled discs(filters)
  21 +/// @param s is the number of scales
  22 +/// @param K is the number of clusters
  23 +
  24 +stim::image<float> tPb(stim::image<float> img, int* r, float* alpha, unsigned int theta_n, unsigned int bin_n, int s, unsigned K){
  25 +
  26 + unsigned int w = img.width(); // get the width of picture
  27 + unsigned int h = img.height(); // get the height of picture
  28 + unsigned int N = w * h; // get the number of pixels
  29 +
  30 + stim::image<float> img_textons(w, h, 1, theta_n*2+1); // allocate space for img_textons
  31 + stim::image<float> img_texture(w, h, 1, 1); // allocate space for img_texture
  32 + stim::image<float> tPb_theta(w, h, 1, 1); // allocate space for tPb_theta
  33 + stim::image<float> tPb(w, h, 1, 1); // allocate space for tPb
  34 + unsigned size = tPb_theta.size(); // get the size of tPb_theta
  35 + memset (tPb.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb to 0
  36 + stim::image<float> temp(w, h, 1, 1); // set the temporary image to store the addtion result
  37 +
  38 + std::ostringstream ss; // (optional) set the stream to designate the test result file
  39 +
  40 +
  41 + img_textons = textons(img, theta_n);
  42 +
  43 + img_texture = kmeans(img_textons, K); // changing kmeans result into float type is required
  44 +
  45 + stim::cpu2image(img_texture.data(), "data_output/texture.bmp", w, h, stim::cmBrewer);
  46 +
  47 +
  48 + unsigned int max1 = img_texture.maxv(); // get the maximum of Pb used for normalization
  49 + unsigned int bin_size = (max1 + 1)/bin_n; // (whether"+1" or not depends on kmeans result)
  50 +
  51 + for (int i = 0; i < theta_n; i++){
  52 +
  53 + float theta = 180 * ((float)i/theta_n); // calculate the even-splited angle for each tPb_theta
  54 +
  55 + memset (tPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb_theta to 0
  56 +
  57 + //ss << "data_output/0922tPb_theta"<< theta << ".bmp"; // set the name for test result file (optional)
  58 + //std::string sss = ss.str();
  59 +
  60 + for (int j = 0; j < s; j++){
  61 +
  62 + // get the chi-gradient by convolving each image slice with the mask
  63 + chi_grad(img_texture.data(), temp.data(), w, h, r[j], bin_n, bin_size, theta);
  64 +
  65 + float max2 = temp.maxv(); // get the maximum of tPb_theta used for normalization
  66 + array_multiply(temp.data(), 1/max2, N); // normalize the tPb_theta
  67 +
  68 + //output the test result of each slice (optional)
  69 + //stim::cpu2image(temp.data(), "data_output/tPb_slice0924_2.bmp", w, h, stim::cmBrewer);
  70 +
  71 + // multiply each chi-gradient with its weight
  72 + array_multiply(temp.data(), alpha[j], N);
  73 +
  74 + // add up all the weighted chi-gradients
  75 + array_add(tPb_theta.data(), temp.data(), tPb_theta.data(), N);
  76 +
  77 +
  78 + }
  79 +
  80 + //ss.str(""); //(optional) clear the space for stream
  81 +
  82 + for(unsigned long ti = 0; ti < N; ti++){
  83 +
  84 + if(tPb_theta.data()[ti] > tPb.data()[ti]){ //get the maximum value among all tPb_theta for ith pixel
  85 + tPb.data()[ti] = tPb_theta.data()[ti];
  86 + }
  87 +
  88 + else{
  89 + }
  90 + }
  91 + }
  92 +
  93 + float max3 = tPb.maxv(); // get the maximum of tPb used for normalization
  94 + array_multiply(tPb.data(), 1/max3, N); // normalize the tPb
  95 +
  96 + //output the test result of tPb (optional)
  97 + //stim::cpu2image(tPb.data(), "data_output/tPb_0922.bmp", w, h, stim::cmBrewer);
  98 +
  99 + return tPb;
  100 +}
  101 +
  102 +#endif
0 \ No newline at end of file 103 \ No newline at end of file
stim/cuda/bsds500/textons.cuh 0 → 100644
  Show/Hide comments   View file @018b00d
  1 +#ifndef STIM_CUDA_TEXONS_CUH
  2 +#define STIM_CUDA_TEXONS_CUH
  3 +
  4 +#include <stim/image/image.h>
  5 +//#include <cmath>
  6 +#include <stim/visualization/colormap.h>
  7 +//#include <stim/image/image_contour_detection.h>
  8 +#include <sstream>
  9 +
  10 +/// This function convolve the grayscale image with a set of oriented Gaussian
  11 +/// derivative filters, and return a texton image with (theta_n*2+1) channels
  12 +
  13 +/// @param image is an one-channel grayscale image
  14 +/// @param theta_n is the number of angles used for computing the gradient
  15 +
  16 +stim::image<float> textons(stim::image<float> image, unsigned int theta_n){
  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 +
  22 + stim::image<float> textons(w, h, 1, theta_n*2+1); // allocate space for textons
  23 + stim::image<float> temp(w, h); // allocate space for temp
  24 +
  25 + int sigma = 1; // set sigma for odd-symmetric, even-symmetric and center-surround filter filter
  26 +
  27 + //std::ostringstream ss1, ss2; // (optional) set the stream to designate the test result file
  28 +
  29 + for (unsigned int i = 0; i < theta_n; i++){
  30 +
  31 + //ss1 << "data_output/textons_channel_"<< i << ".bmp"; // set the name for test result file (optional)
  32 + //std::string sss1 = ss1.str();
  33 + //ss2 << "data_output/textons_channel_"<< i+theta_n << ".bmp"; // set the name for test result file (optional)
  34 + //std::string sss2 = ss2.str();
  35 +
  36 + float theta = 180 * ((float)i/theta_n); // calculate the even-splited angle for each oriented filter
  37 +
  38 + temp = dG1_theta_conv2(image, sigma, theta); // return dG1_theta_conv2 to temp
  39 + //stim::cpu2image(temp.data(), sss1, w, h, stim::cmBrewer);
  40 + textons.set_channel(i, temp.data()); // copy temp to ith channel of textons
  41 +
  42 + temp = dG2_d2x_theta_conv2(image, sigma, theta); // return dG2_d2x_theta_conv2 to temp
  43 + //stim::cpu2image(temp.data(), sss2, w, h, stim::cmBrewer);
  44 + textons.set_channel(i + theta_n, temp.data()); // copy temp to (i+theta_n)th channel of textons
  45 +
  46 + //ss1.str(""); //(optional) clear the space for stream
  47 + //ss2.str(""); //(optional) clear the space for stream
  48 +
  49 + }
  50 +
  51 + temp = laplacian_conv2(image, sigma); // return laplacian_conv2 to temp
  52 + //stim::cpu2image(temp.data(), "data_output/textons_channel_16.bmp", w, h, stim::cmBrewer);
  53 + textons.set_channel(theta_n*2, temp.data()); // copy temp to (theta_n*2)th channel of textons
  54 +
  55 + return textons;
  56 +
  57 +}
  58 +
  59 +#endif
  60 +
  61 +
0 \ No newline at end of file 62 \ No newline at end of file