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