Pb.cpp 3 KB
#include <stim/image/image.h>
//#include <cmath>
#include <stim/visualization/colormap.h>
#include <stim/image/image_contour_detection.h>

#define SIGMA_N 3

void array_abs(float* img, unsigned int N);
void array_multiply(float* lhs, float rhs, unsigned int N);
void array_cos(float* ptr1, float* cpu_out, unsigned int N);
void array_sin(float* ptr1, float* cpu_out, unsigned int N);
void array_atan(float* ptr1, float* cpu_out, unsigned int N);
void array_divide(float* ptr1, float* ptr2,float* cpu_quotient, unsigned int N);
void array_multiply(float* ptr1, float* ptr2, float* product, unsigned int N);
void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);

/// This function uses odd-symmetric gaussian derivative filter to evaluate 
/// the max probability of a contour on one scale, given an one-channel image 

/// @param img is an one-channel image
/// @param r is an array of radii for different scaled discs(filters)
/// @param sigma_n is the number of standard deviations used to define the sigma

stim::image<float> Pb(stim::image<float> image, int sigma){

	unsigned int w = image.width();    // get the width of picture
	unsigned int h = image.height();   // get the height of picture
	unsigned N = w * h;				   // get the number of pixels of picture
	int r = SIGMA_N * sigma;           // set the radius of filter
	int winsize = 2 * r + 1;           // set the winsdow size of filter

	stim::image<float> I(w, h, 1, 2);       // allocate space for return image of dG1_conv2
	stim::image<float> theta(w, h);       // allocate space for theta matrix
	stim::image<float> cos(w, h);       // allocate space for cos(theta)
	stim::image<float> sin(w, h);       // allocate space for sin(theta)
	stim::image<float> temp(w, h);       // allocate space for temp
	stim::image<float> Ix(w, h);       // allocate space for Ix
	stim::image<float> Iy(w, h);       // allocate space for Iy
	stim::image<float> Pb(w, h);       // allocate space for Pb

	I = dG1_conv2(image, sigma);  // calculate the Ix, Iy
	Ix = I.channel(0);
	array_abs(Ix.data(), N);	  //get |Ix|;
	//stim::cpu2image(Ix.data(), "data_output/Pb_Ix_0924.bmp", w, h, stim::cmBrewer); 
	Iy = I.channel(1); 
	array_abs(Iy.data(), N);	  //get |Iy|;
	//stim::cpu2image(Iy.data(), "data_output/Pb_Iy_0924.bmp", w, h, stim::cmBrewer); 

	array_divide(Iy.data(), Ix.data(), temp.data(), N);                //temp = Iy./Ix
	array_atan(temp.data(), theta.data(), N);                   //theta = atan(temp)
	array_cos(theta.data(), cos.data(), N);						//cos = cos(theta)
	array_sin(theta.data(), sin.data(), N);						//sin = sin(theta)
	array_multiply(Ix.data(), cos.data(), Ix.data(), N);		//Ix = Ix.*cos
	array_multiply(Iy.data(), sin.data(), Iy.data(), N);		//Iy = Iy.*sin
	array_add(Ix.data(), Iy.data(), Pb.data(), N);				//Pb = Ix + Iy;

	float max = Pb.maxv();						// get the maximum of Pb used for normalization
	array_multiply(Pb.data(), 1/max, N);		// normalize the Pb

	//stim::cpu2image(Pb.data(), "data_output/Pb_0924.bmp", w, h, stim::cmBrewer); show the Pb(optional)

	return Pb;