Commit 883e39e475893ee4ba64f0536035c978ead82562

Authored by Pavel Govyadinov
2 parents 0df6a853 d24662d2

Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib into Graph

matlab/images2ENVI.m 0 → 100644
  1 +function images2ENVI( image_struct, filename )
  2 +%HSIWRITE : This function generates a hyperspectral binary image and header
  3 +% based on multi-dimention image
  4 +%
  5 +%
  6 +% This function accepts a srutcture which contains several channels.
  7 +% It creates a 'bsq' interleaved image and a ENVI header filer correspoing
  8 +% to the information of the image
  9 +
  10 +fields = fieldnames(image_struct);
  11 +
  12 +% check all images to have the same size
  13 +for i = 1 : size (fields, 1)
  14 + for j = 1 : size (fields, 1)
  15 + if size(image_struct.(fields{i}))~= size(image_struct.(fields{i}))
  16 + fprintf ('Dimension of input images do not Match !');
  17 + exit;
  18 + end
  19 + end
  20 +end
  21 +
  22 +
  23 +
  24 +
  25 +
  26 +% Extracting the needed information
  27 +base_image = image_struct.(fields{1});
  28 +
  29 +rows = size(base_image, 1);
  30 +columns = size(base_image, 2);
  31 +bands = size(fields,1);
  32 +multiDimImage = base_image'; % Matlab stores column major, transposing make it row major to be saved on memory
  33 +
  34 +% Create multi-dimentional image
  35 +for j = 1 : size(fields,1)
  36 + multiDimImage(:,:,j) = (image_struct.(fields{j})'); % Matlab stores column major, transposing make it row major to be saved on memory
  37 +end
  38 +
  39 +% create the binary image
  40 +HSImage = reshape (multiDimImage, [rows*columns*bands, 1]);
  41 +if (~isfloat(HSImage))
  42 + HSImage_float = single(HSImage); % casting to float
  43 +end
  44 +multibandwrite(HSImage_float,filename,'bsq');
  45 +
  46 +
  47 +% create the header file
  48 +info.samples = columns;
  49 +info.lines = rows;
  50 +info.bands = bands;
  51 +info.header_offset = 0;
  52 +info.file_type = 'ENVI Standard';
  53 +info.data_type = 4;
  54 +info.interleave = 'bsq';
  55 +
  56 +
  57 +
  58 +
  59 +fid = fopen(strcat(filename, '.hdr'),'w'); % open a file to write the header file
  60 +fprintf(fid,'%s\n','ENVI'); % first line of the header file showing it is an ENVI file
  61 +
  62 +info_fields=fieldnames(info);
  63 +for i = 1:length(info_fields)
  64 + parameter = info_fields{i}; % get the parameter of each line
  65 + value = info.(info_fields{i}); % get the corresponding value of the parameter
  66 + parameter(strfind(parameter,'_')) = ' '; % change '_' to ' '
  67 + line=[parameter,' = ',num2str(value)]; % create the line of the hdr file
  68 + fprintf(fid,'%s\n',line); % write it to the hdr file
  69 +end
  70 +
  71 +fclose(fid); % clode the header file
  72 +
  73 +end
  74 +
... ...
stim/biomodels/cellset.h
... ... @@ -167,6 +167,26 @@ public:
167 167 cells.push_back(newcell); //push the new memory entry into the cell array
168 168 }
169 169  
  170 +void save(std::string filename){
  171 +
  172 +
  173 + size_t ncells = cells.size(); // get the number of cells
  174 + std::ofstream file(filename); //open a file to store the cell's coordinates
  175 + if (file.is_open()) {
  176 +
  177 + file << "x y z radius\n"; //write the file header
  178 + for (size_t c=0; c < ncells; c++){ //for each cell
  179 + if (cells[c][ip[3]] != NULL) //check if for the current cell, radius has been assigned
  180 + file << cells[c][ip[0]] << delim << cells[c][ip[1]] << delim << cells[c][ip[2]] << delim << cells[c][ip[3]] << '\n' ;
  181 + else //check if for the current cell, radius has not been assigned, set radius to zero
  182 + file << cells[c][ip[0]] << delim << cells[c][ip[1]] << delim << cells[c][ip[2]] << delim << 0 << '\n' ;
  183 + }
  184 +
  185 + }
  186 + file.close();
  187 +
  188 + }
  189 +
170 190  
171 191 }; //end class cellset
172 192 }; //end namespace stim
... ...
stim/cuda/sharedmem.cuh
... ... @@ -48,7 +48,11 @@ namespace stim{
48 48  
49 49 /// Threaded copying of 2D data on a CUDA device
50 50 /// @param dest is a linear destination array of size nx * ny
  51 + /// @param nx is the size of the region to be copied along the X dimension
  52 + /// @param ny is the size of the region to be copied along the Y dimension
51 53 /// @param src is a 2D image stored as a linear array with a pitch of X
  54 + /// @param x is the x position in the source image where the copy is started
  55 + /// @param y is the y position in the source image where the copy is started
52 56 /// @param X is the number of bytes in a row of src
53 57 /// @param tid is a 1D id for the current thread
54 58 /// @param nt is the number of threads in the block
... ...
stim/image/image.h
1 1 #ifndef STIM_IMAGE_H
2 2 #define STIM_IMAGE_H
3 3  
4   -#include <opencv2/core/core.hpp>
5   -#include <opencv2/highgui/highgui.hpp>
  4 +#ifdef USING_OPENCV
  5 + #include <opencv2/core/core.hpp>
  6 + #include <opencv2/highgui/highgui.hpp>
  7 +#endif
6 8 #include <vector>
7 9 #include <iostream>
8 10 #include <limits>
9 11 #include <typeinfo>
  12 +#include <fstream>
10 13  
11 14 namespace stim{
12 15 /// This static class provides the STIM interface for loading, saving, and storing 2D images.
... ... @@ -17,13 +20,12 @@ namespace stim{
17 20 template <class T>
18 21 class image{
19 22  
20   - //cimg_library::CImg<T> img;
21   - T* img; //pointer to the image data (assumes RGB for loading/saving)
  23 + T* img; //pointer to the image data (interleaved RGB for color)
22 24 size_t R[3];
23 25  
24   - size_t X() const { return R[1]; }
25   - size_t Y() const { return R[2]; }
26   - size_t C() const { return R[0]; }
  26 + inline size_t X() const { return R[1]; }
  27 + inline size_t Y() const { return R[2]; }
  28 + inline size_t C() const { return R[0]; }
27 29  
28 30 void init(){ //initializes all variables, assumes no memory is allocated
29 31 memset(R, 0, sizeof(size_t) * 3); //set the resolution and number of channels to zero
... ... @@ -52,22 +54,12 @@ class image{
52 54  
53 55 size_t bytes(){ return size() * sizeof(T); }
54 56  
55   - size_t idx(size_t x, size_t y, size_t c = 0){
56   - return y * C() * X() + x * C() + c;
  57 + inline size_t idx(size_t x, size_t y, size_t c = 0) const {
  58 + return y * R[0] * R[1] + x * R[0] + c;
57 59 }
58 60  
59   -
  61 +#ifdef USING_OPENCV
60 62 int cv_type(){
61   - // The following is C++ 11 code, but causes problems on some compilers (ex. nvcc). Below is my best approximation to a solution
62   -
63   - //if(std::is_same<T, unsigned char>::value) return CV_MAKETYPE(CV_8U, (int)C());
64   - //if(std::is_same<T, char>::value) return CV_MAKETYPE(CV_8S, (int)C());
65   - //if(std::is_same<T, unsigned short>::value) return CV_MAKETYPE(CV_16U, (int)C());
66   - //if(std::is_same<T, short>::value) return CV_MAKETYPE(CV_16S, (int)C());
67   - //if(std::is_same<T, int>::value) return CV_MAKETYPE(CV_32S, (int)C());
68   - //if(std::is_same<T, float>::value) return CV_MAKETYPE(CV_32F, (int)C());
69   - //if(std::is_same<T, double>::value) return CV_MAKETYPE(CV_64F, (int)C());
70   -
71 63 if(typeid(T) == typeid(unsigned char)) return CV_MAKETYPE(CV_8U, (int)C());
72 64 if(typeid(T) == typeid(char)) return CV_MAKETYPE(CV_8S, (int)C());
73 65 if(typeid(T) == typeid(unsigned short)) return CV_MAKETYPE(CV_16U, (int)C());
... ... @@ -79,18 +71,9 @@ class image{
79 71 std::cout<<"ERROR in stim::image::cv_type - no valid data type found"<<std::endl;
80 72 exit(1);
81 73 }
82   -
  74 +#endif
83 75 /// Returns the value for "white" based on the dynamic range (assumes white is 1.0 for floating point images)
84 76 T white(){
85   - // The following is C++ 11 code, but causes problems on some compilers (ex. nvcc). Below is my best approximation to a solution
86   -
87   - //if(std::is_same<T, unsigned char>::value) return UCHAR_MAX;
88   - //if(std::is_same<T, unsigned short>::value) return SHRT_MAX;
89   - //if(std::is_same<T, unsigned>::value) return UINT_MAX;
90   - //if(std::is_same<T, unsigned long>::value) return ULONG_MAX;
91   - //if(std::is_same<T, unsigned long long>::value) return ULLONG_MAX;
92   - //if(std::is_same<T, float>::value) return 1.0f;
93   - //if(std::is_same<T, double>::value) return 1.0;
94 77  
95 78 if(typeid(T) == typeid(unsigned char)) return UCHAR_MAX;
96 79 if(typeid(T) == typeid(unsigned short)) return SHRT_MAX;
... ... @@ -142,6 +125,11 @@ public:
142 125 free(img);
143 126 }
144 127  
  128 + ///Resize an image
  129 + void resize(size_t x, size_t y, size_t c = 1) {
  130 + allocate(x, y, c);
  131 + }
  132 +
145 133 stim::image<T>& operator=(const stim::image<T>& I){
146 134 init();
147 135 if(&I == this) //handle self-assignment
... ... @@ -151,9 +139,88 @@ public:
151 139 return *this;
152 140 }
153 141  
  142 + //save a Netpbm file
  143 + void load_netpbm(std::string filename) {
  144 + std::ifstream infile(filename, std::ios::in | std::ios::binary); //open an output file
  145 + if (!infile) {
  146 + std::cout << "Error opening input file in image::load_netpbm()" << std::endl;
  147 + exit(1);
  148 + }
  149 + if (sizeof(T) != 1) {
  150 + std::cout << "Error in image::load_netpbm() - data type must be 8-bit integer." << std::endl;
  151 + exit(1);
  152 + }
  153 +
  154 + size_t nc; //allocate space for the number of channels
  155 + char format[2]; //allocate space to hold the image format tag
  156 + infile.read(format, 2); //read the image format tag
  157 + infile.seekg(1, std::ios::cur); //skip the newline character
  158 +
  159 + if (format[0] != 'P') {
  160 + std::cout << "Error in image::load_netpbm() - file format tag is invalid: " << format[0] << format[1] << std::endl;
  161 + exit(1);
  162 + }
  163 + if (format[1] == '5') nc = 1; //get the number of channels from the format flag
  164 + else if (format[1] == '6') nc = 3;
  165 + else {
  166 + std::cout << "Error in image::load_netpbm() - file format tag is invalid: " << format[0] << format[1] << std::endl;
  167 + exit(1);
  168 + }
  169 +
  170 + unsigned char c; //stores a character
  171 + while (infile.peek() == '#') { //if the next character indicates the start of a comment
  172 + while (true) {
  173 + c = infile.get();
  174 + if (c == 0x0A) break;
  175 + }
  176 + }
  177 +
  178 + std::string sw; //create a string to store the width of the image
  179 + while(true){
  180 + c = infile.get(); //get a single character
  181 + if (c == ' ') break; //exit if we've encountered a space
  182 + sw.push_back(c); //push the character on to the string
  183 + }
  184 + size_t w = atoi(sw.c_str()); //convert the string into an integer
  185 +
  186 + std::string sh;
  187 + while (true) {
  188 + c = infile.get();
  189 + if (c == 0x0A) break;
  190 + sh.push_back(c);
  191 + }
  192 +
  193 + while (true) { //skip the maximum value
  194 + c = infile.get();
  195 + if (c == 0x0A) break;
  196 + }
  197 + size_t h = atoi(sh.c_str()); //convert the string into an integer
  198 +
  199 + allocate(w, h, nc); //allocate space for the image
  200 + infile.read((char*)img, size()); //copy the binary data from the file to the image
  201 + infile.close();
  202 + }
  203 +
  204 +
  205 +#ifdef USING_OPENCV
  206 + void from_opencv(unsigned char* buffer, size_t width, size_t height) {
  207 + allocate(width, height, 3);
  208 + T value;
  209 + size_t i;
  210 + for (size_t c = 0; c < C(); c++) { //copy directly
  211 + for (size_t y = 0; y < Y(); y++) {
  212 + for (size_t x = 0; x < X(); x++) {
  213 + i = y * X() * C() + x * C() + (2 - c);
  214 + value = buffer[i];
  215 + img[idx(x, y, c)] = value;
  216 + }
  217 + }
  218 + }
  219 + }
  220 +#endif
154 221 /// Load an image from a file
155 222 void load(std::string filename){
156   -
  223 +#ifdef USING_OPENCV
157 224 cv::Mat cvImage = cv::imread(filename, CV_LOAD_IMAGE_UNCHANGED); //use OpenCV to open the image file
158 225 if(!cvImage.data){
159 226 std::cout<<"ERROR stim::image::load() - unable to find image "<<filename<<std::endl;
... ... @@ -168,25 +235,42 @@ public:
168 235 memcpy(img, cv_ptr, bytes());
169 236 if(C() == 3) //if this is a 3-color image, OpenCV uses BGR interleaving
170 237 from_opencv(cv_ptr, X(), Y());
  238 +#else
  239 + load_netpbm(filename);
  240 +#endif
171 241 }
172 242  
173   - void from_opencv(unsigned char* buffer, size_t width, size_t height){
174   - allocate(width, height, 3);
175   - T value;
176   - size_t i;
177   - for(size_t c = 0; c < C(); c++){ //copy directly
178   - for(size_t y = 0; y < Y(); y++){
179   - for(size_t x = 0; x < X(); x++){
180   - i = y * X() * C() + x * C() + (2-c);
181   - value = buffer[i];
182   - img[idx(x, y, c)] = value;
183   - }
184   - }
  243 +
  244 +
  245 + //save a Netpbm file
  246 + void save_netpbm(std::string filename) {
  247 + std::ofstream outfile(filename, std::ios::out | std::ios::binary); //open an output file
  248 + if(!outfile) {
  249 + std::cout << "Error generating output file in image::save_netpbm()" << std::endl;
  250 + exit(1);
  251 + }
  252 + if (sizeof(T) != 1) {
  253 + std::cout << "Error in image::save_netpbm() - data type must be 8-bit integer." << std::endl;
  254 + exit(1);
185 255 }
  256 + std::string format;
  257 + if (channels() == 1) outfile << "P5" << (char)0x0A; //output P5 if the file is grayscale
  258 + else if (channels() == 3) outfile << "P6" << (char)0x0A; //output P6 if the file is color
  259 + else {
  260 + std::cout << "Error in image::save_netpbm() - data must be grayscale or RGB." << std::endl;
  261 + exit(1);
  262 + }
  263 + size_t w = width();
  264 + size_t h = height();
  265 + outfile << w << " " << h << (char)0x0A; //save the width and height
  266 + outfile << "255" << (char)0x0A; //output the maximum value
  267 + outfile.write((const char*)img, size()); //write the binary data
  268 + outfile.close();
186 269 }
187 270  
188 271 //save a file
189 272 void save(std::string filename){
  273 +#ifdef USING_OPENCV
190 274 //OpenCV uses an interleaved format, so convert first and then output
191 275 T* buffer = (T*) malloc(bytes());
192 276  
... ... @@ -197,6 +281,9 @@ public:
197 281 cv::Mat cvImage((int)Y(), (int)X(), cv_type(), buffer);
198 282 cv::imwrite(filename, cvImage);
199 283 free(buffer);
  284 +#else
  285 + save_netpbm(filename);
  286 +#endif
200 287 }
201 288  
202 289 void set_interleaved(T* buffer, size_t width, size_t height, size_t channels){
... ... @@ -256,20 +343,36 @@ public:
256 343 }
257 344 }
258 345  
259   -
260   - image<T> channel(size_t c){
261   -
262   - //create a new image
263   - image<T> r(X(), Y(), 1);
264   -
  346 + /// Return an image representing a specified channel
  347 + /// @param c is the channel to be returned
  348 + image<T> channel(size_t c) const {
  349 + image<T> r(X(), Y(), 1); //create a new image
265 350 for(size_t x = 0; x < X(); x++){
266 351 for(size_t y = 0; y < Y(); y++){
267 352 r.img[r.idx(x, y, 0)] = img[idx(x, y, c)];
268 353 }
269 354 }
  355 + return r;
  356 + }
  357 +
  358 + /// Returns an std::vector containing each channel as a separate image
  359 + std::vector<image<T>> split() const {
  360 + std::vector<image<T>> r; //create an image array
  361 + r.resize(C()); //create images for each channel
270 362  
  363 + for (size_t c = 0; c < C(); c++) { //for each channel
  364 + r[c] = channel(c); //copy the channel image to the array
  365 + }
271 366 return r;
  367 + }
272 368  
  369 + /// Merge a series of single-channel images into a multi-channel image
  370 + void merge(std::vector<image<T>>& list) {
  371 + size_t x = list[0].width(); //calculate the size of the image
  372 + size_t y = list[0].height();
  373 + allocate(x, y, list.size()); //re-allocate the image
  374 + for (size_t c = 0; c < list.size(); c++) //for each channel
  375 + set_channel(list[c].channel(0).data(), c); //insert the channel into the output image
273 376 }
274 377  
275 378 T& operator()(size_t x, size_t y, size_t c = 0){
... ... @@ -300,20 +403,20 @@ public:
300 403 size_t x, y;
301 404 for(y = 0; y < Y(); y++){
302 405 for(x = 0; x < X(); x++){
303   - img[idx(x, y, c)] = buffer[c];
  406 + img[idx(x, y, c)] = buffer[y * X() + x];
304 407 }
305 408 }
306 409 }
307 410  
308   - size_t channels(){
  411 + size_t channels() const{
309 412 return C();
310 413 }
311 414  
312   - size_t width(){
  415 + size_t width() const{
313 416 return X();
314 417 }
315 418  
316   - size_t height(){
  419 + size_t height() const{
317 420 return Y();
318 421 }
319 422  
... ... @@ -423,6 +526,14 @@ public:
423 526 exit(1);
424 527 }
425 528  
  529 + /// Casting operator, casts every value in an image to a different data type V
  530 + template<typename V>
  531 + operator image<V>() {
  532 + image<V> r(X(), Y(), C()); //create a new image
  533 + std::copy(img, img + size(), r.data()); //copy and cast the data
  534 + return r; //return the new image
  535 + }
  536 +
426 537 };
427 538  
428 539 }; //end namespace stim
... ...
stim/math/filters/conv2.h 0 → 100644
  1 +#ifndef STIM_CUDA_CONV2_H
  2 +#define STIM_CUDA_CONV2_H
  3 +//#define __CUDACC__
  4 +
  5 +#ifdef __CUDACC__
  6 +#include <stim/cuda/cudatools.h>
  7 +#include <stim/cuda/sharedmem.cuh>
  8 +#endif
  9 +
  10 +namespace stim {
  11 +#ifdef __CUDACC__
  12 + //Kernel function that performs the 2D convolution.
  13 + template<typename T, typename K = T>
  14 + __global__ void kernel_conv2(T* out, T* in, K* kernel, size_t sx, size_t sy, size_t kx, size_t ky) {
  15 + extern __shared__ T s[]; //declare a shared memory array
  16 + size_t xi = blockIdx.x * blockDim.x + threadIdx.x; //threads correspond to indices into the output image
  17 + size_t yi = blockIdx.y * blockDim.y + threadIdx.y;
  18 + size_t tid = threadIdx.y * blockDim.x + threadIdx.x;
  19 + size_t nt = blockDim.x * blockDim.y;
  20 +
  21 + size_t cx = blockIdx.x * blockDim.x; //find the upper left corner of the input region
  22 + size_t cy = blockIdx.y * blockDim.y;
  23 +
  24 + size_t X = sx - kx + 1; //calculate the size of the output image
  25 + size_t Y = sy - ky + 1;
  26 +
  27 + if (cx >= X || cy >= Y) return; //return if the entire block is outside the image
  28 + size_t smx = min(blockDim.x + kx - 1, sx - cx); //size of the shared copy of the input image
  29 + size_t smy = min(blockDim.y + ky - 1, sy - cy); // min function is used to deal with boundary blocks
  30 + stim::cuda::threadedMemcpy2D<T>(s, smx, smy, in, cx, cy, sx, sy, tid, nt); //copy the input region to shared memory
  31 + __syncthreads();
  32 +
  33 + if (xi >= X || yi >= Y) return; //returns if the thread is outside of the output image
  34 +
  35 + //loop through the kernel
  36 + size_t kxi, kyi;
  37 + K v = 0;
  38 + for (kyi = 0; kyi < ky; kyi++) {
  39 + for (kxi = 0; kxi < kx; kxi++) {
  40 + v += s[(threadIdx.y + kyi) * smx + threadIdx.x + kxi] * kernel[kyi * kx + kxi];
  41 + //v += in[(yi + kyi) * sx + xi + kxi] * kernel[kyi * kx + kxi];
  42 + }
  43 + }
  44 + out[yi * X + xi] = (T)v; //write the result to global memory
  45 +
  46 + }
  47 +
  48 + //Performs a convolution of a 2D image using the GPU. All pointers are assumed to be to memory on the current device.
  49 + //@param out is a pointer to the output image
  50 + //@param in is a pointer to the input image
  51 + //@param sx is the size of the input image along X
  52 + //@param sy is the size of the input image along Y
  53 + //@param kx is the size of the kernel along X
  54 + //@param ky is the size of the kernel along Y
  55 + template<typename T, typename K = T>
  56 + void gpu_conv2(T* out, T* in, K* kernel, size_t sx, size_t sy, size_t kx, size_t ky) {
  57 + cudaDeviceProp p;
  58 + HANDLE_ERROR(cudaGetDeviceProperties(&p, 0));
  59 + size_t tmax = p.maxThreadsPerBlock;
  60 + dim3 nt(sqrt(tmax), sqrt(tmax)); //calculate the block dimensions
  61 + size_t X = sx - kx + 1; //calculate the size of the output image
  62 + size_t Y = sy - ky + 1;
  63 + dim3 nb(X / nt.x + 1, Y / nt.y + 1); //calculate the grid dimensions
  64 + size_t sm = (nt.x + kx - 1) * (nt.y + ky - 1) * sizeof(T); //shared memory bytes required to store block data
  65 + if (sm > p.sharedMemPerBlock) {
  66 + std::cout << "Error in stim::gpu_conv2() - insufficient shared memory for this kernel." << std::endl;
  67 + exit(1);
  68 + }
  69 + kernel_conv2 <<<nb, nt, sm>>> (out, in, kernel, sx, sy, kx, ky); //launch the kernel
  70 + }
  71 +#endif
  72 + //Performs a convolution of a 2D image. Only valid pixels based on the kernel are returned.
  73 + // As a result, the output image will be smaller than the input image by (kx-1, ky-1)
  74 + //@param out is a pointer to the output image
  75 + //@param in is a pointer to the input image
  76 + //@param sx is the size of the input image along X
  77 + //@param sy is the size of the input image along Y
  78 + //@param kx is the size of the kernel along X
  79 + //@param ky is the size of the kernel along Y
  80 + template<typename T, typename K = T>
  81 + void cpu_conv2(T* out, T* in, K* kernel, size_t sx, size_t sy, size_t kx, size_t ky) {
  82 + size_t X = sx - kx + 1; //x size of the output image
  83 + size_t Y = sy - ky + 1; //y size of the output image
  84 +
  85 +#ifdef __CUDACC__
  86 + //allocate memory and copy everything to the GPU
  87 + T* gpu_in;
  88 + HANDLE_ERROR(cudaMalloc(&gpu_in, sx * sy * sizeof(T)));
  89 + HANDLE_ERROR(cudaMemcpy(gpu_in, in, sx * sy * sizeof(T), cudaMemcpyHostToDevice));
  90 + K* gpu_kernel;
  91 + HANDLE_ERROR(cudaMalloc(&gpu_kernel, kx * ky * sizeof(K)));
  92 + HANDLE_ERROR(cudaMemcpy(gpu_kernel, kernel, kx * ky * sizeof(K), cudaMemcpyHostToDevice));
  93 + T* gpu_out;
  94 + HANDLE_ERROR(cudaMalloc(&gpu_out, X * Y * sizeof(T)));
  95 + gpu_conv2(gpu_out, gpu_in, gpu_kernel, sx, sy, kx, ky); //execute the GPU kernel
  96 + HANDLE_ERROR(cudaMemcpy(out, gpu_out, X * Y * sizeof(T), cudaMemcpyDeviceToHost)); //copy the result to the host
  97 + HANDLE_ERROR(cudaFree(gpu_in));
  98 + HANDLE_ERROR(cudaFree(gpu_kernel));
  99 + HANDLE_ERROR(cudaFree(gpu_out));
  100 +#else
  101 + K v; //register stores the integral of the current pixel value
  102 + size_t yi, xi, kyi, kxi, yi_kyi_sx;
  103 + for (yi = 0; yi < Y; yi++) { //for each pixel in the output image
  104 + for (xi = 0; xi < X; xi++) {
  105 + v = 0;
  106 + for (kyi = 0; kyi < ky; kyi++) { //for each pixel in the kernel
  107 + yi_kyi_sx = (yi + kyi) * sx;
  108 + for (kxi = 0; kxi < kx; kxi++) {
  109 + v += in[yi_kyi_sx + xi + kxi] * kernel[kyi * kx + kxi];
  110 + }
  111 + }
  112 + out[yi * X + xi] = v; //save the result to the output array
  113 + }
  114 + }
  115 +
  116 +#endif
  117 + }
  118 +
  119 +
  120 +}
  121 +
  122 +
  123 +#endif
0 124 \ No newline at end of file
... ...
stim/math/filters/gauss2.h 0 → 100644
  1 +#ifndef STIM_CUDA_GAUSS2_H
  2 +#define STIM_CUDA_GAUSS2_H
  3 +
  4 +#include <stim/image/image.h>
  5 +#include <stim/math/filters/sepconv2.h>
  6 +#include <stim/math/constants.h>
  7 +
  8 +namespace stim {
  9 +
  10 + template<typename T>
  11 + T gauss1d(T x, T mu, T std) {
  12 + return ((T)1 / (T)sqrt(stim::TAU * std * std) * exp(-(x - mu)*(x - mu) / (2 * std*std)));
  13 + }
  14 + ///Perform a 2D gaussian convolution on an input image.
  15 + ///@param in is a pointer to the input image
  16 + ///@param stdx is the standard deviation (in pixels) along the x axis
  17 + ///@param stdy is the standard deviation (in pixels) along the y axis
  18 + ///@param nstds specifies the number of standard deviations of the Gaussian that will be kept in the kernel
  19 + template<typename T, typename K = T>
  20 + stim::image<T> cpu_gauss2(const stim::image<T>& in, K stdx, K stdy, size_t nstds = 3) {
  21 + size_t kx = stdx * nstds * 2; //calculate the kernel sizes
  22 + size_t ky = stdy * nstds * 2;
  23 + size_t X = in.width() - kx + 1; //calculate the size of the output image
  24 + size_t Y = in.height() - ky + 1;
  25 + stim::image<T> r(X, Y, in.channels()); //create an output image
  26 +
  27 + K* gx = (K*)malloc(kx * sizeof(K)); //allocate space for the gaussian kernels
  28 + K* gy = (K*)malloc(ky * sizeof(K));
  29 + K mux = (K)kx / (K)2; //calculate the mean of the gaussian (center of the kernel)
  30 + K muy = (K)ky / (K)2;
  31 + for (size_t xi = 0; xi < kx; xi++) //evaluate the kernels
  32 + gx[xi] = gauss1d((K)xi, mux, stdx);
  33 + for (size_t yi = 0; yi < ky; yi++)
  34 + gy[yi] = gauss1d((K)yi, muy, stdy);
  35 +
  36 + std::vector<stim::image<T>> IN = in.split(); //split the input image into channels
  37 + std::vector<stim::image<T>> R = r.split(); //split the output image into channels
  38 + for (size_t c = 0; c < in.channels(); c++) //for each channel
  39 + cpu_sepconv2(R[c].data(), IN[c].data(), gx, gy, IN[c].width(), IN[c].height(), kx, ky);
  40 +
  41 + r.merge(R); //merge the blurred channels into the final image
  42 + return r;
  43 + }
  44 +}
  45 +
  46 +
  47 +#endif
0 48 \ No newline at end of file
... ...
stim/math/filters/sepconv2.h 0 → 100644
  1 +#ifndef STIM_CUDA_SEPCONV2_H
  2 +#define STIM_CUDA_SEPCONV2_H
  3 +#include <stim/math/filters/conv2.h>
  4 +#ifdef __CUDACC__
  5 +#include <stim/cuda/cudatools.h>
  6 +#include <stim/cuda/sharedmem.cuh>
  7 +#endif
  8 +
  9 +namespace stim {
  10 +#ifdef __CUDACC__
  11 + //Performs a convolution of a 2D image using the GPU. All pointers are assumed to be to memory on the current device.
  12 + //@param out is a pointer to the output image
  13 + //@param in is a pointer to the input image
  14 + //@param sx is the size of the input image along X
  15 + //@param sy is the size of the input image along Y
  16 + //@param kx is the size of the kernel along X
  17 + //@param ky is the size of the kernel along Y
  18 + template<typename T, typename K = T>
  19 + void gpu_sepconv2(T* out, T* in, K* k0, K* k1, size_t sx, size_t sy, size_t kx, size_t ky) {
  20 + cudaDeviceProp p;
  21 + HANDLE_ERROR(cudaGetDeviceProperties(&p, 0));
  22 + size_t tmax = p.maxThreadsPerBlock;
  23 + dim3 nt(sqrt(tmax), sqrt(tmax)); //calculate the block dimensions
  24 + size_t X = sx - kx + 1; //calculate the x size of the output image
  25 + T* temp; //declare a temporary variable to store the intermediate image
  26 + HANDLE_ERROR(cudaMalloc(&temp, X * sy * sizeof(T))); //allocate memory for the intermediate image
  27 +
  28 + dim3 nb(X / nt.x + 1, sy / nt.y + 1); //calculate the grid dimensions
  29 + size_t sm = (nt.x + kx - 1) * nt.y * sizeof(T); //shared memory bytes required to store block data
  30 + if (sm > p.sharedMemPerBlock) {
  31 + std::cout << "Error in stim::gpu_conv2() - insufficient shared memory for this kernel." << std::endl;
  32 + exit(1);
  33 + }
  34 + kernel_conv2 <<<nb, nt, sm>>> (temp, in, k0, sx, sy, kx, 1); //launch the kernel to compute the intermediate image
  35 +
  36 + size_t Y = sy - ky + 1; //calculate the y size of the output image
  37 + nb.y = Y / nt.y + 1; //update the grid dimensions to reflect the Y-axis size of the output image
  38 + sm = nt.x * (nt.y + ky - 1) * sizeof(T); //calculate the amount of shared memory needed for the second pass
  39 + if (sm > p.sharedMemPerBlock) {
  40 + std::cout << "Error in stim::gpu_conv2() - insufficient shared memory for this kernel." << std::endl;
  41 + exit(1);
  42 + }
  43 + kernel_conv2 <<<nb, nt, sm>>> (out, temp, k1, X, sy, 1, ky); //launch the kernel to compute the final image
  44 + HANDLE_ERROR(cudaFree(temp)); //free memory allocated for the intermediate image
  45 + }
  46 +#endif
  47 + //Performs a separable convolution of a 2D image. Only valid pixels based on the kernel are returned.
  48 + // As a result, the output image will be smaller than the input image by (kx-1, ky-1)
  49 + //@param out is a pointer to the output image
  50 + //@param in is a pointer to the input image
  51 + //@param k0 is the x-axis convolution filter
  52 + //@param k1 is the y-axis convolution filter
  53 + //@param sx is the size of the input image along X
  54 + //@param sy is the size of the input image along Y
  55 + //@param kx is the size of the kernel along X
  56 + //@param ky is the size of the kernel along Y
  57 + template<typename T, typename K = T>
  58 + void cpu_sepconv2(T* out, T* in, K* k0, K* k1, size_t sx, size_t sy, size_t kx, size_t ky) {
  59 + size_t X = sx - kx + 1; //x size of the output image
  60 + size_t Y = sy - ky + 1;
  61 +#ifdef __CUDACC__
  62 + //allocate memory and copy everything to the GPU
  63 + T* gpu_in;
  64 + HANDLE_ERROR(cudaMalloc(&gpu_in, sx * sy * sizeof(T)));
  65 + HANDLE_ERROR(cudaMemcpy(gpu_in, in, sx * sy * sizeof(T), cudaMemcpyHostToDevice));
  66 + K* gpu_k0;
  67 + HANDLE_ERROR(cudaMalloc(&gpu_k0, kx * sizeof(K)));
  68 + HANDLE_ERROR(cudaMemcpy(gpu_k0, k0, kx * sizeof(K), cudaMemcpyHostToDevice));
  69 + K* gpu_k1;
  70 + HANDLE_ERROR(cudaMalloc(&gpu_k1, ky * sizeof(K)));
  71 + HANDLE_ERROR(cudaMemcpy(gpu_k1, k1, ky * sizeof(K), cudaMemcpyHostToDevice));
  72 + T* gpu_out;
  73 + HANDLE_ERROR(cudaMalloc(&gpu_out, X * Y * sizeof(T)));
  74 + gpu_sepconv2(gpu_out, gpu_in, gpu_k0, gpu_k1, sx, sy, kx, ky); //execute the GPU kernel
  75 + HANDLE_ERROR(cudaMemcpy(out, gpu_out, X * Y * sizeof(T), cudaMemcpyDeviceToHost)); //copy the result to the host
  76 + HANDLE_ERROR(cudaFree(gpu_in));
  77 + HANDLE_ERROR(cudaFree(gpu_k0));
  78 + HANDLE_ERROR(cudaFree(gpu_k1));
  79 + HANDLE_ERROR(cudaFree(gpu_out));
  80 +#else
  81 + T* temp = (T*)malloc(X * sy * sizeof(T)); //allocate space for the intermediate image
  82 + cpu_conv2(temp, in, k0, sx, sy, kx, 1); //evaluate the intermediate image
  83 + cpu_conv2(out, temp, k1, X, sy, 1, ky); //evaluate the final image
  84 + free(temp); //free the memory for the intermediate image
  85 +#endif
  86 + }
  87 +}
  88 +
  89 +#endif
0 90 \ No newline at end of file
... ...
stim/parser/filename.h
... ... @@ -302,6 +302,15 @@ public:
302 302 return insert(ss.str());
303 303 }
304 304  
  305 + ///This method returns true if any characters in the filename contain '*' or '?'
  306 + bool wildcards() {
  307 + if (_prefix.find('*') != std::string::npos) return true;
  308 + if (_prefix.find('?') != std::string::npos) return true;
  309 + if (_extension.find('*') != std::string::npos) return true;
  310 + if (_extension.find('?') != std::string::npos) return true;
  311 + return false;
  312 + }
  313 +
305 314  
306 315 /// Returns a list of files using the current filename as a template.
307 316 /// For example:
... ...