diff --git a/CMakeLists.txt b/CMakeLists.txt index f4ef136..3350281 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,9 @@ project(bsds500) #set the module directory set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}") +#find OpenCV +find_package(OpenCV REQUIRED ) + #set up CUDA find_package(CUDA REQUIRED) @@ -20,6 +23,7 @@ find_package(Threads) find_package(X11) include_directories( + ${OpenCV_INCLUDE_DIRS} ${STIM_INCLUDE_DIRS} ) @@ -41,6 +45,7 @@ cuda_add_executable(bsds500 target_link_libraries(bsds500 #${CUDA_cufft_LIBRARY} #${CUDA_cublas_LIBRARY} + ${OpenCV_LIBS} ${CMAKE_THREAD_LIBS_INIT} ${X11_LIBRARIES} ) @@ -48,5 +53,6 @@ target_link_libraries(bsds500 #copy an image test case configure_file(data/101085.bmp 101085.bmp COPYONLY) configure_file(data/101087.bmp 101087.bmp COPYONLY) +configure_file(data/101087_square.bmp 101087_square.bmp COPYONLY) configure_file(data/slice00.bmp slice00.bmp COPYONLY) configure_file(data/slice00_500_500.bmp slice00_500_500.bmp COPYONLY) diff --git a/Pb.cpp b/Pb.cpp new file mode 100644 index 0000000..96683ba --- /dev/null +++ b/Pb.cpp @@ -0,0 +1,61 @@ +#include +//#include +#include +#include + +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 r, unsigned int sigma_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 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter) + + stim::image I(w, h, 1, 2); // allocate space for return image of Gd1 + 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 = Gd1(image, r, sigma_n); // 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/cPb.cpp b/cPb.cpp new file mode 100644 index 0000000..3d9f23c --- /dev/null +++ b/cPb.cpp @@ -0,0 +1,67 @@ +#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* r, 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 + int sigma_n = 3; // set the number of standard deviations used to define the sigma + + 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), r[i*s + j], sigma_n); + + // 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/cudafunc.cu b/cudafunc.cu index afb7cdb..58921dd 100644 --- a/cudafunc.cu +++ b/cudafunc.cu @@ -1,9 +1,7 @@ #include - -/*void blur(float* image, float sigma, unsigned int x, unsigned int y){ - - stim::cuda::cpu_gaussian_blur_2d(image, sigma, x, y); -}*/ +#include +#include +#include void array_multiply(float* lhs, float rhs, unsigned int N){ @@ -16,6 +14,12 @@ void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N){ } +void array_multiply(float* ptr1, float* ptr2, float* product, unsigned int N){ + + stim::cuda::cpu_add(ptr1, ptr2, product, N); + +} + void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M){ stim::cuda::cpu_conv2(img, mask, cpu_copy, w, h, M); @@ -28,3 +32,41 @@ void array_abs(float* img, unsigned int N){ } +void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1){ + + stim::cuda::cpu_conv2sep(img, x, y, kernel0, k0, kernel1, k1); + +} + +void array_cos(float* ptr1, float* cpu_out, unsigned int N){ + + stim::cuda::cpu_cos(ptr1, cpu_out, N); + +} + +void array_sin(float* ptr1, float* cpu_out, unsigned int N){ + + stim::cuda::cpu_sin(ptr1, cpu_out, N); + +} + +void array_atan(float* ptr1, float* cpu_out, unsigned int N){ + + stim::cuda::cpu_atan(ptr1, cpu_out, N); + +} + +void array_divide(float* ptr1, float* ptr2,float* cpu_quotient, unsigned int N){ + + stim::cuda::cpu_divide(ptr1, ptr2, cpu_quotient, 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){ + + stim::cuda::cpu_chi_grad(img, cpu_copy, w, h, r, bin_n, bin_size, theta); + +} + + + diff --git a/fun_mPb_theta.cpp b/fun_mPb_theta.cpp deleted file mode 100644 index 3bff266..0000000 --- a/fun_mPb_theta.cpp +++ /dev/null @@ -1,62 +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 theta-dependent multicale Pb given an multi-channel image - -/// @param img is the multi-channel image -/// @param theta is the angle used for computing the 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 - -stim::image func_mPb_theta(stim::image img, float theta, int* r, 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 mPb_theta(w, h, 1); // allocate space for theta-dependent multicale Pb - unsigned size = mPb_theta.size(); // get the size of theta-dependent multicale Pb - memset ( mPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of theta-dependent multicale Pb to 0 - - - unsigned int N = w * h; // get the number of pixels - int sigma_n = 3; // set the number of standard deviations used to define the sigma - - 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/mPb_theta_slice"<< i*s + j << ".bmp"; // set the name for test result file - //std::string sss = ss.str(); - - // get the gaussian gradient by convolving each image slice with the mask - temp = gaussian_derivative_filter_odd(img.channel(i), r[i*s + j], sigma_n, theta); - - // (optional) output the test result of each slice - //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(mPb_theta.data(), temp.data(), mPb_theta.data(), N); - - //ss.str(""); //(optional) clear the space for stream - - } - } - - - return mPb_theta; -} diff --git a/func_mPb.cpp b/func_mPb.cpp deleted file mode 100644 index c8011b7..0000000 --- a/func_mPb.cpp +++ /dev/null @@ -1,79 +0,0 @@ -#include -#include -#include -#include -#include - -/// This function evaluates the multicale Pb given an multi-channel image - -/// @param img is the multi-channel image -/// @param theta_n is the number of angles used for computing the gradient -/// @param r is an array of radii for different scaled discs(filters) -/// @param alpha is an array of weights for different scaled discs(filters) -/// @param s is the number of scales -stim::image func_mPb(stim::image img, unsigned int theta_n, int* r, float* alpha, int s){ - - std::clock_t start; // (optional) set the timer to calculate the total time - start = std::clock(); // (optional) set timer start point - - - std::ostringstream ss; // (optional) set the stream to designate the test result file - - 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 mPb_theta(w,h); // allocate space for theta-dependent multicale Pb (mPb_theta) - stim::image mPb(w,h); // allocate space for multicale Pb (mPb) - - unsigned size = mPb.size(); // get the size of mPb - memset ( mPb.data(), 0, size * sizeof(float)); // initialize all the pixels of mPb to 0 - - float* ptr; // set a pointer - ptr = (float*) malloc(size * sizeof(float) * theta_n); // this pointer points to a continuous space allocated to store all the mPb_theta - - for (unsigned int n = 0; n < theta_n; n++){ - - ss << "data_output/mPb_theta"<< n << "_0911.bmp"; // (optional) set the name for test result file - std::string sss = ss.str(); // (optional) - float theta = 180 * ((float)n/theta_n); // calculate the even-splited angle for each mPb_theta - - mPb_theta = func_mPb_theta(img, theta, r, alpha, s); // calculate the mPb_theta - - float* ptr_n = &ptr[ n * w * h * 1 ]; // set a pointer which points to the space for each mPb_theta - mPb_theta.data_noninterleaved(ptr_n); // set this pointer to point to the each mPb_theta - - double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the time for generating each mPb_theta - std::cout<<"mPb_theta_"<< theta <<" complished time:"<< duration1 <<"s"<<'\n'; // (optional) show this time - - - unsigned long idx = n * w * h * 1; //index for the nth mPb_theta - - - stim::cpu2image(mPb_theta.data(), sss, w, h, stim::cmBrewer); // (optional) output the nth mPb_theta - - for(unsigned long i = 0; i < N; i++){ - - float pixel = ptr[i+idx]; //get the ith pixel in nth mPb_theta - - if(pixel > mPb.data()[i]){ //get the maximum value among all mPb_theta for ith pixel - mPb.data()[i] = pixel; - } - - else{ - } - } - - - ss.str(""); //(optional) clear the space for stream - - } - - stim::cpu2image(mPb.data(), "data_output/mPb_500_0911_neat.bmp", w, h, stim::cmBrewer); // output the mPb - - double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time - std::cout<<"total time:"<< duration2 <<"s"<<'\n'; // (optional) show the total time - - return mPb; - -} diff --git a/gauss_deriv1.cpp b/gauss_deriv1.cpp new file mode 100644 index 0000000..ed50990 --- /dev/null +++ b/gauss_deriv1.cpp @@ -0,0 +1,80 @@ +#include +//#include +#include +#include + +/// 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 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 + +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 Gd1(stim::image image, int r, unsigned int sigma_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 + 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 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/gauss_deriv2.cpp b/gauss_deriv2.cpp new file mode 100644 index 0000000..85414b1 --- /dev/null +++ b/gauss_deriv2.cpp @@ -0,0 +1,80 @@ +#include +//#include +#include +#include + +/// 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 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 + +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 Gd2(stim::image image, int r, unsigned int sigma_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 + 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 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/gauss_deriv_center.cpp b/gauss_deriv_center.cpp new file mode 100644 index 0000000..4b41d0d --- /dev/null +++ b/gauss_deriv_center.cpp @@ -0,0 +1,41 @@ +#include +#include +#include +#include + +#define PI 3.1415926 + +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 Gd_center(stim::image image, int r, unsigned int sigma_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 + int winsize = 2 * r + 1; // set the winsdow size of disc(filter) + + stim::image I(w, h, 1, 2); // allocate space for return image of Gd2 + stim::image Ixx(w, h); // allocate space for Ixx + stim::image Iyy(w, h); // allocate space for Iyy + stim::image Gd_center(w, h); // allocate space for Pb + + I = Gd2(image, r, sigma_n); // calculate the Ixx, Iyy + Ixx = I.channel(0); + Iyy = I.channel(1); + + array_add(Ixx.data(), Iyy.data(), Gd_center.data(), N); //Gd_center = Ixx + Iyy; + array_abs(Gd_center.data(), N); + + //stim::cpu2image(Gd_center.data(), "data_output/Gd_center_0919.bmp", w, h, stim::cmBrewer); + + return Gd_center; + +} diff --git a/gauss_deriv_even.cpp b/gauss_deriv_even.cpp new file mode 100644 index 0000000..a268d95 --- /dev/null +++ b/gauss_deriv_even.cpp @@ -0,0 +1,52 @@ +#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/gauss_deriv_odd.cpp b/gauss_deriv_odd.cpp new file mode 100644 index 0000000..b1ea9aa --- /dev/null +++ b/gauss_deriv_odd.cpp @@ -0,0 +1,48 @@ +#include +#include +#include +#include + +#define PI 3.1415926 + +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 Gd_odd(stim::image image, int r, unsigned int sigma_n, 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 + unsigned N = w * h; // get the number of pixels of picture + int winsize = 2 * r + 1; // set the winsdow size of disc(filter) + + stim::image I(w, h, 1, 2); // allocate space for return image of Gd1 + stim::image Ix(w, h); // allocate space for Ix + stim::image Iy(w, h); // allocate space for Iy + stim::image Gd_odd_theta(w, h); // allocate space for Pb + + I = Gd1(image, r, sigma_n); // 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(), Gd_odd_theta.data(), N); //Gd_odd_theta = Ix + Iy; + array_abs(Gd_odd_theta.data(), N); + + //stim::cpu2image(I.channel(0).data(), "data_output/Gd_odd_x_0919.bmp", w, h, stim::cmBrewer); + //stim::cpu2image(I.channel(1).data(), "data_output/Gd_odd_y_0919.bmp", w, h, stim::cmBrewer); + //stim::cpu2image(Gd_odd_theta.data(), "data_output/Gd_odd_theta_0919.bmp", w, h, stim::cmBrewer); + + return Gd_odd_theta; + +} diff --git a/gauss_derivative_odd.cpp b/gauss_derivative_odd.cpp deleted file mode 100644 index ca98c41..0000000 --- a/gauss_derivative_odd.cpp +++ /dev/null @@ -1,64 +0,0 @@ -#include -#include -#include - -#define PI 3.1415926 - -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); -void array_multiply(float* lhs, float rhs, unsigned int N); - -/// This function evaluates the gaussian derivative gradient given 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 gaussian_derivative_filter_odd(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 mask_x(winsize, winsize); // allocate space for x-axis-oriented filter - stim::image mask_y(winsize, winsize); // allocate space for y-axis-oriented filter - stim::image mask_theta(winsize, winsize);// allocate space for theta-oriented filter - stim::image derivative_theta(w, h); // allocate space for theta-oriented gradient - - float theta_r = (theta * PI)/180; //change angle unit from degree to rad - - 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) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); - // create the y-oriented gaussian derivative filter mask_y - mask_y.data()[j*winsize + i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); - // create the mask_theta - mask_theta.data()[j*winsize + i] = cos(theta_r) * mask_x.data()[j*winsize + i] + sin(theta_r) * mask_y.data()[j*winsize + i] ; - - } - } - - //stim::cpu2image(mask_theta.data(), "data_output/mask_0911_2.bmp", winsize, winsize, stim::cmBrewer); // (optional) show the mask result - - // do the 2D convolution with image and mask - conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize); - - array_abs(derivative_theta.data(), N); // get the absolute value for each pixel (why slower than the "for loop" method sometimes?) - - float max = derivative_theta.max(); // get the maximum of gradient used for normalization - array_multiply(derivative_theta.data(), 1/max, N); // normalize the gradient - - - //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta_0911.bmp", w, h, stim::cmBrewer); // (optional) show the gradient result - - return derivative_theta; - -} diff --git a/kmeans.cpp b/kmeans.cpp new file mode 100644 index 0000000..fee17b7 --- /dev/null +++ b/kmeans.cpp @@ -0,0 +1,67 @@ +#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/main.cpp b/main.cpp new file mode 100644 index 0000000..60bb004 --- /dev/null +++ b/main.cpp @@ -0,0 +1,92 @@ +#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 + + 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 r1[3] = {3,5,10}; // set an array of radii for different scaled discs(filters) for cPb, the length of r is c*s + 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 + 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::cout<<"begin cPb"<<'\n'; + + + std::clock_t start1; // (optional) set the timer to calculate the total time + start1 = std::clock(); // (optional) set timer start point + + img_cPb = cPb(image, r1, alpha, s); + + // show the cPb (optional) + stim::cpu2image(img_cPb.data(), "data_output/img_cPb_0925.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 + + + //*******************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, r2, alpha, theta_n, bin_n, s, K); + + // show the tPb (optional) + stim::cpu2image(img_tPb.data(), "data_output/img_tPb_0925.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); + + // show the mPb (optional) + stim::cpu2image(img_mPb.data(), "data_output/img_mPb_0925.bmp", w, h, stim::cmBrewer); + + return 0; + +} diff --git a/old version/fun_mPb_theta.cpp b/old version/fun_mPb_theta.cpp new file mode 100644 index 0000000..da603b9 --- /dev/null +++ b/old version/fun_mPb_theta.cpp @@ -0,0 +1,63 @@ +#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 theta-dependent multicale Pb given an multi-channel image + +/// @param img is the multi-channel image +/// @param theta is the angle used for computing the 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 + +stim::image func_mPb_theta(stim::image img, float theta, int* r, 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 mPb_theta(w, h, 1); // allocate space for theta-dependent multicale Pb + unsigned size = mPb_theta.size(); // get the size of theta-dependent multicale Pb + memset ( mPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of theta-dependent multicale Pb to 0 + + + unsigned int N = w * h; // get the number of pixels + int sigma_n = 3; // set the number of standard deviations used to define the sigma + + 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/mPb_theta_slice"<< i*s + j << ".bmp"; // set the name for test result file + //std::string sss = ss.str(); + + // get the gaussian gradient by convolving each image slice with the mask + temp = gaussian_derivative_filter_odd(img.channel(i), r[i*s + j], sigma_n, theta); + + // (optional) output the test result of each slice + //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(mPb_theta.data(), temp.data(), mPb_theta.data(), N); + + //ss.str(""); //(optional) clear the space for stream + + } + } + + + return mPb_theta; +} diff --git a/old version/func_mPb.cpp b/old version/func_mPb.cpp new file mode 100644 index 0000000..f104934 --- /dev/null +++ b/old version/func_mPb.cpp @@ -0,0 +1,80 @@ +#include +#include +#include +#include +#include + +/// This function evaluates the multicale Pb given an multi-channel image + +/// @param img is the multi-channel image +/// @param theta_n is the number of angles used for computing the gradient +/// @param r is an array of radii for different scaled discs(filters) +/// @param alpha is an array of weights for different scaled discs(filters) +/// @param s is the number of scales + +stim::image func_mPb(stim::image img, unsigned int theta_n, int* r, float* alpha, int s){ + + std::clock_t start; // (optional) set the timer to calculate the total time + start = std::clock(); // (optional) set timer start point + + + std::ostringstream ss; // (optional) set the stream to designate the test result file + + 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 mPb_theta(w,h); // allocate space for theta-dependent multicale Pb (mPb_theta) + stim::image mPb(w,h); // allocate space for multicale Pb (mPb) + + unsigned size = mPb.size(); // get the size of mPb + memset ( mPb.data(), 0, size * sizeof(float)); // initialize all the pixels of mPb to 0 + + float* ptr; // set a pointer + ptr = (float*) malloc(size * sizeof(float) * theta_n); // this pointer points to a continuous space allocated to store all the mPb_theta + + for (unsigned int n = 0; n < theta_n; n++){ + + //ss << "data_output/mPb_theta"<< n << "_0911.bmp"; // (optional) set the name for test result file + //std::string sss = ss.str(); // (optional) + float theta = 180 * ((float)n/theta_n); // calculate the even-splited angle for each mPb_theta + + mPb_theta = func_mPb_theta(img, theta, r, alpha, s); // calculate the mPb_theta + + float* ptr_n = &ptr[ n * w * h * 1 ]; // set a pointer which points to the space for each mPb_theta + mPb_theta.data_noninterleaved(ptr_n); // set this pointer to point to the each mPb_theta + + double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the time for generating each mPb_theta + std::cout<<"mPb_theta_"<< theta <<" complished time:"<< duration1 <<"s"<<'\n'; // (optional) show this time + + + unsigned long idx = n * w * h * 1; //index for the nth mPb_theta + + + //stim::cpu2image(mPb_theta.data(), sss, w, h, stim::cmBrewer); // (optional) output the nth mPb_theta + + for(unsigned long i = 0; i < N; i++){ + + float pixel = ptr[i+idx]; //get the ith pixel in nth mPb_theta + + if(pixel > mPb.data()[i]){ //get the maximum value among all mPb_theta for ith pixel + mPb.data()[i] = pixel; + } + + else{ + } + } + + + //ss.str(""); //(optional) clear the space for stream + + } + + stim::cpu2image(mPb.data(), "data_output/mPb_500_0914_neat.bmp", w, h, stim::cmBrewer); // output the mPb + + double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time + std::cout<<"total time:"<< duration2 <<"s"<<'\n'; // (optional) show the total time + + return mPb; + +} diff --git a/old version/gauss_derivative_odd.cpp b/old version/gauss_derivative_odd.cpp new file mode 100644 index 0000000..dbb37dc --- /dev/null +++ b/old version/gauss_derivative_odd.cpp @@ -0,0 +1,91 @@ +#include +#include +#include + +#define PI 3.1415926 + +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); +void array_multiply(float* lhs, float rhs, unsigned int N); + +/// This function evaluates the gaussian derivative gradient given 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 gaussian_derivative_filter_odd(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 mask_x(winsize, winsize); // allocate space for x-axis-oriented filter + stim::image mask_y(winsize, winsize); // allocate space for y-axis-oriented filter + stim::image mask_theta(winsize, winsize);// allocate space for theta-oriented filter + stim::image derivative_theta(w, h); // allocate space for theta-oriented gradient + + stim::image mask_r; + + + float theta_r = (theta * PI)/180; //change angle unit from degree to rad + + //*************inseparable convolution**************** + 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) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); + // create the y-oriented gaussian derivative filter mask_y + mask_y.data()[j*winsize + i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); + // create the mask_theta + mask_theta.data()[j*winsize + i] = cos(theta_r) * mask_x.data()[j*winsize + i] + sin(theta_r) * mask_y.data()[j*winsize + i] ; + + } + } + + mask_r = mask_x.rotate(90, r, r); + stim::cpu2image(mask_r.data(), "data_output/mask_r90_0915.bmp", winsize, winsize, stim::cmBrewer); + //stim::cpu2image(mask_theta.data(), "data_output/mask_0911_2.bmp", winsize, winsize, stim::cmBrewer); // (optional) show the mask result + + // do the 2D convolution with image and mask + conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize); + //*********************************************************/ + + /*************separable convolution**************** + 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_x + array_x[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); + // create the 1D y-oriented gaussian derivative filter array_y + array_y[i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); + + + } + + //stim::cpu2image(mask_theta.data(), "data_output/mask_0911_2.bmp", winsize, winsize, stim::cmBrewer); // (optional) show the mask result + + // do the 2D convolution with image and mask + conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize); + //*********************************************************/ + + array_abs(derivative_theta.data(), N); // get the absolute value for each pixel (why slower than the "for loop" method sometimes?) + + float max = derivative_theta.maxv(); // get the maximum of gradient used for normalization + array_multiply(derivative_theta.data(), 1/max, N); // normalize the gradient + + + //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta_0911.bmp", w, h, stim::cmBrewer); // (optional) show the gradient result + + return derivative_theta; + +} diff --git a/old version/test_main - Copy.cpp b/old version/test_main - Copy.cpp new file mode 100644 index 0000000..8980a02 --- /dev/null +++ b/old version/test_main - Copy.cpp @@ -0,0 +1,37 @@ +#include +#include +#include +#include +#include +#include + +/// calculate the mPb given a multi-channel image + +int main() +{ + stim::image img; // generate an image object + //stim::image img_r; + + //img.load("slice00_500_500.bmp"); // load the input image + img.load("101087_square.bmp"); // load the input image + img = img.channel(0); // get the first channel of black-and-white image + + unsigned int w = img.width(); // get the width of picture + unsigned int h = img.height(); // get the height of picture + int c = img.channels(); // get the number if channels of picture + int s = 3; // set the number of scales + + //img_r = img.rotate(45, w/2, h/2); + //stim::cpu2image(img_r.data(), "data_output/rotate_0915.bmp", w, h, stim::cmBrewer); + + + int r[3] = {3,5,10}; // set an array of radii for different scaled discs(filters) + 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 for computing the gradient + + stim::image mPb; // allocate the space for mPb + mPb = func_mPb(img, theta_n, r, alpha, s); // calculate the mPb + + return 0; + +} diff --git a/tPb.cpp b/tPb.cpp new file mode 100644 index 0000000..e066dae --- /dev/null +++ b/tPb.cpp @@ -0,0 +1,97 @@ +#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_0925.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/test_main.cpp b/test_main.cpp deleted file mode 100644 index 49cb37c..0000000 --- a/test_main.cpp +++ /dev/null @@ -1,29 +0,0 @@ -#include -#include -#include -#include -#include -/// calculate the mPb given a multi-channel image - -int main() -{ - stim::image img; // generate an image object - - img.load("slice00_500_500.bmp"); // load the input image - img = img.channel(0); // get the first channel of black-and-white image - - unsigned int w = img.width(); // get the width of picture - unsigned int h = img.height(); // get the height of picture - int c = img.channels(); // get the number if channels of picture - int s = 3; // set the number of scales - - int r[3] = {3,5,10}; // set an array of radii for different scaled discs(filters) - 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 for computing the gradient - - stim::image mPb; // allocate the space for mPb - mPb = func_mPb(img, theta_n, r, alpha, s); // calculate the mPb - - return 0; - -} diff --git a/textons.cpp b/textons.cpp new file mode 100644 index 0000000..dbe66f1 --- /dev/null +++ b/textons.cpp @@ -0,0 +1,61 @@ +#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 + + unsigned int r_odd = 3; // set disc radii for odd-symmetric filter + unsigned int sigma_n_odd = 3; // set sigma_n for odd-symmetric filter + unsigned int r_even = 3; // set disc radii for even-symmetric filter + unsigned int sigma_n_even = 3; // set sigma_n for even-symmetric filter + unsigned int r_center = 3; // set disc radii for center-surround filter + unsigned int sigma_n_center = 3; // set sigma_n for center-surround 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 = Gd_odd(image, r_odd, sigma_n_odd, theta); // return Gd_odd 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 = Gd_even(image, r_even, sigma_n_even, theta); // return Gd_even 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 = Gd_center(image, r_center, sigma_n_center); // return Gd_center 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 -- libgit2 0.21.4