#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; }