segmentation / bsds500

  #include <stim/image/image.h>
  #include <cmath>
  #include <stim/visualization/colormap.h>

  
  #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<float> gaussian_derivative_filter_odd(stim::image<float> 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<float> mask_x(winsize, winsize);    // allocate space for x-axis-oriented filter
  	stim::image<float> mask_y(winsize, winsize);    // allocate space for y-axis-oriented filter
  	stim::image<float> mask_theta(winsize, winsize);// allocate space for theta-oriented filter
  	stim::image<float> derivative_theta(w, h);      // allocate space for theta-oriented gradient

  	stim::image<float> 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;
  

  }