diff --git a/CMakeLists.txt b/CMakeLists.txt index 6b65367..f4ef136 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,3 +48,5 @@ 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/slice00.bmp slice00.bmp COPYONLY) +configure_file(data/slice00_500_500.bmp slice00_500_500.bmp COPYONLY) diff --git a/fun_mPb_theta.cpp b/fun_mPb_theta.cpp index 532e74c..3bff266 100644 --- a/fun_mPb_theta.cpp +++ b/fun_mPb_theta.cpp @@ -1,102 +1,62 @@ #include #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 gradient image given an input image +/// This function evaluates the theta-dependent multicale Pb given an multi-channel image -/// @param lab is the 3-channel image in the LAB color space +/// @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 + + } + } -stim::image func_mPb_theta(stim::image lab, float theta, unsigned int w, unsigned int h){ - - //allocate space for the gradient image - stim::image mPb_theta(w, h, 1); - - //allocate space for each individual channel - stim::image pic_light, pic_colora, pic_colorb; - pic_light = lab.channel(0); - pic_colora = lab.channel(1); - pic_colorb = lab.channel(2); - - unsigned int N = w * h; //calculate the number of pixels in the image - float sigma = 2; //set the sigma value to \sigma = 2 - unsigned int sigma_n = 3; //set the number of standard deviations used to define the kernel size - unsigned r1 = 3; //disk radii - unsigned r2 = 5; - unsigned r3 = 10; - unsigned r4 = 20; - float alpha[9] = {1,1,1,1,1,1,1,1,1}; //weighting for each channel - - - stim::image l1,l2,l3,a1,a2,a3,b1,b2,b3; - - l1 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r3 * 2, theta, w, h); - stim::cpu2image(l1.data(), "data_output/testl_tex5.bmp", w, h, stim::cmBrewer); - - exit(0); - /*l2 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r2 * 2, theta, w, h); - stim::cpu2image(l2.data(), "data_output/l2_tex.bmp", w, h, stim::cmBrewer); - l3 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r3 * 2, theta, w, h); - stim::cpu2image(l3.data(), "data_output/l3_tex.bmp", w, h, stim::cmBrewer); - a1 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r2 * 2, theta, w, h); - stim::cpu2image(a1.data(), "data_output/a1_tex.bmp", w, h, stim::cmBrewer); - a2 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r3 * 2, theta, w, h); - stim::cpu2image(a2.data(), "data_output/a2_tex.bmp", w, h, stim::cmBrewer); - a3 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r4 * 2, theta, w, h); - stim::cpu2image(a3.data(), "data_output/a3_tex.bmp", w, h, stim::cmBrewer); - b1 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r2 * 2, theta, w, h); - stim::cpu2image(b1.data(), "data_output/b1_tex.bmp", w, h, stim::cmBrewer); - b2 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r3 * 2, theta, w, h); - stim::cpu2image(b2.data(), "data_output/b2_tex.bmp", w, h, stim::cmBrewer); - b3 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r4 * 2, theta, w, h); - stim::cpu2image(b3.data(), "data_output/b3_tex.bmp", w, h, stim::cmBrewer);*/ - - /*for (unsigned i = 0; i #include -//#include #include #include #include -stim::image func_mPb(stim::image lab, unsigned int theta_n, unsigned int w, unsigned int h){ +/// This function evaluates the multicale Pb given an multi-channel image - std::clock_t start; - start = std::clock(); +/// @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){ - //---------------pavel's suggesiton------------------------------------ - std::ostringstream ss; - unsigned int N = w * h; - stim::image mPb_theta(w,h), mPb(w,h); - unsigned size = mPb_theta.size(); - memset ( mPb.data(), 0, size * sizeof(float)); + std::clock_t start; // (optional) set the timer to calculate the total time + start = std::clock(); // (optional) set timer start point - float* ptr; - ptr = (float*) malloc(size * sizeof(float) * theta_n); + + 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 << "_conv2.bmp"; - float theta = 180 * ((float)n/theta_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(lab, theta, w, h); - //mPb_theta.load("101087.bmp"); - float* ptr_n = &ptr[ n * w * h * 1 ]; - mPb_theta.channel(0).data_noninterleaved(ptr_n); + mPb_theta = func_mPb_theta(img, theta, r, alpha, s); // calculate the mPb_theta - double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; - std::cout<<"mPb_theta_"<< theta <<" complished time:"<< duration1 <<"s"<<'\n'; + 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 slice - std::string sss = ss.str(); - //stim::cpu2image(&ptr[idx], sss, w, h, stim::cmBrewer); - + 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 slice + float pixel = ptr[i+idx]; //get the ith pixel in nth mPb_theta - if(pixel > mPb.data()[i]){ + if(pixel > mPb.data()[i]){ //get the maximum value among all mPb_theta for ith pixel mPb.data()[i] = pixel; } @@ -53,83 +65,15 @@ stim::image func_mPb(stim::image lab, unsigned int theta_n, unsign } + ss.str(""); //(optional) clear the space for stream - ss.str(""); } - //stim::cpu2image(mPb.data(), "data_output/mPb_conv2.bmp", w, h, stim::cmBrewer); + 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; - std::cout<<"total time:"<< duration2 <<"s"<<'\n'; - - //getch(); + 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; - //---------------my first method------------------------------------ - /* - std::clock_t start; - start = std::clock(); - - stim::image mPb_stack(w,h,theta_n), mPb(w,h), mPb_theta(w,h), A, B, temp; - float* ptr[8]; - - for (unsigned int n = 0; n < theta_n; n++){ - - //int* x = new int(5); - //int* y = x; - //*y = 1; - - float theta = 180 * ((float)n/theta_n); - mPb_theta = func_mPb_theta(lab, theta, w, h); - mPb_stack.getslice(n) = mPb_theta; - float* ptr[n] = mPb_stack.getslice(n).data(); - - double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; - std::cout<<"mPb_theta, theta = "<< theta <<" time:"<< duration1 <<"s"<<'\n'; - - - for(unsigned long i = 0; i < N; i++){ - - *(ptr[n]+i) = mPb_theta.data()[i]; - - - //float a = mPb_theta.data()[i]; - //float* B = ptr[n]+i; - //A.data()[i] = mPb_theta.data()[i]; - //float* C = ptr[0]+1; - //*C = 1; - - // - } - stim::cpu2image(ptr[0], "data_output/mPb_theta.bmp", w, h, stim::cmBrewer); - } - - for (unsigned long i = 0; i < N; i++){ - - mPb.data()[i] = 0; - for (unsigned int n = 0; n < theta_n; n++){ - - float* ptr2 = ptr[i]+n; - float temp = *ptr2; - - if(temp > mPb.data()[i]){ - mPb.data()[i] = temp; - } - else{ - } - } - } - - stim::cpu2image(mPb.data(), "data_output/cmap_mPb.bmp", w, h, stim::cmBrewer); - - double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; - std::cout<<"total time:"<< duration2 <<"s"<<'\n'; - - getch(); - - return mPb; */ - - - } diff --git a/gauss_derivative_odd.cpp b/gauss_derivative_odd.cpp index 6fe8de7..ca98c41 100644 --- a/gauss_derivative_odd.cpp +++ b/gauss_derivative_odd.cpp @@ -1,7 +1,6 @@ #include #include #include -//#include #define PI 3.1415926 @@ -9,80 +8,56 @@ void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned in void array_abs(float* img, unsigned int N); void array_multiply(float* lhs, float rhs, unsigned int N); -// winsize = 2 * r, side of mask = winsize + 1 -stim::image gaussian_derivative_filter_odd(stim::image image, float sigma, unsigned int sigma_n, unsigned int winsize, float theta, unsigned int w, unsigned int h){ +/// This function evaluates the gaussian derivative gradient given an one-channel image - stim::image mask_x(winsize+1, winsize+1), mask_y(winsize+1, winsize+1), mask_theta(winsize+1, winsize+1), mask_delta(winsize+1, winsize+1, 1), derivative_x, derivative_y, derivative_theta(w, h); - //float* ptr = mask_x.data(); +/// @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){ - //DEBUG calculate a Dirac delta function kernel - memset ( mask_delta.data(), 0, mask_delta.size() * sizeof(float)); - mask_delta.data()[winsize*(winsize+2)/2] = 1; - stim::cpu2image(mask_delta.data(), "data_output/mask_test.bmp", winsize+1, winsize+1, stim::cmBrewer); + 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 - // set parameters - unsigned N = w * h; - float theta_r = (theta * PI)/180; + float theta_r = (theta * PI)/180; //change angle unit from degree to rad - float step = (2*sigma*sigma_n)/winsize; + for (int j = 0; j < winsize; j++){ + for (int i = 0; i< winsize; i++){ - for (unsigned j = 0; j <= winsize; j++){ - for (unsigned i = 0; i<= winsize; i++){ - - float x = (-1)*sigma*sigma_n + i * step; //range of x - float y = (-1)*sigma*sigma_n + j * step; //range of y + 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+1) + i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))); + 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+1) + i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))); + 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+1) + i] = cos(theta_r) * mask_x.data()[j*(winsize+1) + i] + sin(theta_r) * mask_y.data()[j*(winsize+1) + i] ; + 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_x.data(), "data_output/cmapgray_mask_x.bmp", winsize+1, winsize+1, stim::cmBrewer); - - stim::cpu2image(image.data(), "data_output/image.bmp", w, h, stim::cmBrewer); - - - stim::cpu2image(mask_theta.data(), "data_output/mask.bmp", winsize+1, winsize+1, stim::cmBrewer); - - // 2D convolution - //derivative_theta = image.convolve2(mask_theta); - //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta1.bmp", w, h, stim::cmBrewer); - conv2(image.data(), mask_delta.data(), derivative_theta.data(), w, h, winsize+1); - //conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize+1); - stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta_tex1.bmp", w, h, stim::cmBrewer); - - //array_abs(derivative_theta.data(), N); - - /*for (unsigned k = 0; k < w * h; k++){ - - derivative_theta.data()[k] = abs(derivative_theta.data()[k]); - - }*/ - - //stim::cpu2image(derivative_theta.data(), "data_output/derivative_theta2_abs.bmp", w, h, stim::cmBrewer); + //stim::cpu2image(mask_theta.data(), "data_output/mask_0911_2.bmp", winsize, winsize, stim::cmBrewer); // (optional) show the mask result - /*float max = derivative_theta.max(); + // do the 2D convolution with image and mask + conv2(image.data(), mask_theta.data(), derivative_theta.data(), w, h, winsize); - array_multiply(derivative_theta.data(), 1/max, N);*/ - - /*( - for (unsigned k = 0; k < w * h; k++){ - - derivative_theta.data()[k] = derivative_theta.data()[k]/max; + 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 - //float max2 = derivative_theta.max(); - //stim::cpu2image(derivative_theta.data(), "data_output/cmap_colorb_gradient_theta90_r5.bmp", w, h, stim::cmBrewer); - //derivative_x.save("data_output/gradient_x.bmp"); + //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/image_contour_detection.h b/image_contour_detection.h deleted file mode 100644 index 57ce07c..0000000 --- a/image_contour_detection.h +++ /dev/null @@ -1,5 +0,0 @@ -#include -//#include -//#include - -stim::image gaussian_derivative_filter_odd(stim::image image, float sigma, unsigned int sigma_n, unsigned int winsize, float theta, unsigned int w, unsigned int h); \ No newline at end of file diff --git a/test_main.cpp b/test_main.cpp index a8451ac..49cb37c 100644 --- a/test_main.cpp +++ b/test_main.cpp @@ -3,53 +3,26 @@ #include #include #include - +/// calculate the mPb given a multi-channel image int main() { - stim::image rgb,gaussgradient; //generate an image object - - //unsigned int a = 5%5; - //unsigned int b = 5/5; - - rgb.load("101087.bmp"); //load the input image - unsigned int w = rgb.width(); //get the image size - unsigned int h = rgb.height(); - unsigned int s = rgb.size(); - //unsigned a = sizeof(float); - - stim::image lab; //create an image object for a single-channel (grayscale) image - lab = rgb.srgb2lab(); //create the single-channel image - - /* - stim::image pic_light, pic_colora, pic_colorb; - pic_light = lab.channel(0); - pic_light.save("pic_light.bmp"); - - pic_colora = lab.channel(1); - pic_colorb = lab.channel(2); - - float sigma = 2; - unsigned int sigma_n = 3; - unsigned int r = 5; - unsigned int winsize = r * 2; //window size = winsize + 1 - float theta = 90; - - gaussgradient = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, winsize, theta, w, h); - gaussgradient.save("data_output/pic_gray_gradient.bmp"); - */ - - //float theta = 0; - unsigned int theta_n = 8; - - //stim::image mPb_stack(w,h,theta_n); - - //stim::image mPb_theta; - //mPb_theta = func_mPb_theta(lab, theta, w, h); - //mPb_theta.save("data_output/pic_gray_gradient.bmp"); - - stim::image mPb; - mPb = func_mPb(lab, theta_n, w, h); + 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; -- libgit2 0.21.4