diff --git a/stim/cuda/imageproc/gaussian_blur.cuh b/stim/cuda/imageproc/gaussian_blur.cuh deleted file mode 100644 index 99056af..0000000 --- a/stim/cuda/imageproc/gaussian_blur.cuh +++ /dev/null @@ -1,141 +0,0 @@ -#ifndef STIM_CUDA_GAUSSIAN_BLUR_H -#define STIM_CUDA_GAUSSIAN_BLUR_H - -#include -#include -#include "gaussian_blur.cuh" -#include - -#define pi 3.14159 - -// This CUDA kernel applies a Gaussian blur along the x axis -template -__global__ void gaussian_blur_x(T* out, T* in, T sigma, unsigned int x, unsigned int y){ - - //calculate the 1D image index for this thread - int i = blockIdx.x * blockDim.x + threadIdx.x; - - //convert this to 2D pixel coordinates - int yi = i / x; - int xi = i - (yi * x); - - int k_size = sigma * 4; //calculate the kernel size - - int gx, gy, gi; //variables to store global coordinates - T sum = 0; //variable stores the integral of the kernel times the image - T G; //stores the value of the gaussian kernel - - out[i] = 0; - //for each element of the kernel - for(int ki = -k_size; ki <= k_size; ki++){ - - //calculate the gaussian value - G = 1.0 / (sigma * sqrt(2 * pi)) * exp(-(ki*ki) / (2*sigma*sigma)); - - //calculate the global coordinates for this point in the kernel - gx = xi + ki; - gy = yi; - - //make sure we are inside the bounds of the image - if(gx >= 0 && gy >= 0 && gx < x && gy < y){ - gi = gy * x + gx; - - //perform the integration - sum += G * in[gi]; - } - } - - //set the output value to the integral of the function times the kernel - out[i] = sum; -} - -// This CUDA kernel applies a Gaussian blur along the y axis. -template -__global__ void gaussian_blur_y(T* out, T* in, T sigma, unsigned int x, unsigned int y){ - - //calculate the 1D image index for this thread - int i = blockIdx.x * blockDim.x + threadIdx.x; - - //convert this to 2D pixel coordinates - int yi = i / x; - int xi = i - (yi * x); - - int k_size = sigma * 4; //calculate the kernel size - - int gx, gy, gi; //variables to store global coordinates - T sum = 0; //variable stores the integral of the kernel times the image - T G; //stores the value of the gaussian kernel - - out[i] = 0; - //for each element of the kernel - for(int ki = -k_size; ki <= k_size; ki++){ - - //calculate the gaussian value - G = 1.0 / (sigma * sqrt(2 * pi)) * exp(-(ki*ki) / (2*sigma*sigma)); - - //calculate the global coordinates for this point in the kernel - gx = xi; - gy = yi + ki; - - //make sure we are inside the bounds of the image - if(gx >= 0 && gy >= 0 && gx < x && gy < y){ - gi = gy * x + gx; - - //perform the integration - sum += G * in[gi]; - } - } - - //set the output value to the integral of the function times the kernel - out[i] = sum; -} - -/// Applies a Gaussian blur to a 2D image stored on the GPU -template -void gpu_gaussian_blur_2d(T* image, T sigma, unsigned int x, unsigned int y){ - - //get the number of pixels in the image - unsigned int pixels = x * y; - unsigned int bytes = sizeof(T) * pixels; - - //allocate temporary space on the GPU - T* gpuI1; - cudaMalloc(&gpuI1, bytes); - - //get the maximum number of threads per block for the CUDA device - int threads = stim::maxThreadsPerBlock(); - - //calculate the number of blocks - int blocks = pixels / threads + (pixels%threads == 0 ? 0:1); - - //run the kernel for the x dimension - gaussian_blur_x <<< blocks, threads >>>(gpuI1, image, sigma, x, y); - - //run the kernel for the y dimension - gaussian_blur_y <<< blocks, threads >>>(image, gpuI1, sigma, x, y); - -} - -/// Applies a Gaussian blur to a 2D image stored on the CPU -template -void cpu_gaussian_blur_2d(T* image, T sigma, unsigned int x, unsigned int y){ - - //get the number of pixels in the image - unsigned int pixels = x * y; - unsigned int bytes = sizeof(T) * pixels; - - //allocate space on the GPU - T* gpuI0; - cudaMalloc(&gpuI0, bytes); - - //copy the image data to the GPU - cudaMemcpy(gpuI0, image, bytes, cudaMemcpyHostToDevice); - - //run the GPU-based version of the algorithm - gpu_gaussian_blur_2d(gpuI0, sigma, x, y); - - //copy the image data from the device - cudaMemcpy(image, gpuI0, bytes, cudaMemcpyDeviceToHost); -} - -#endif diff --git a/stim/envi/binary.h b/stim/envi/binary.h index ff68201..16ff0b9 100644 --- a/stim/envi/binary.h +++ b/stim/envi/binary.h @@ -100,7 +100,7 @@ public: /// @param filename is the name of the binary file /// @param r is a STIM vector specifying the size of the binary file along each dimension /// @param h is the length (in bytes) of any header file (default zero) - bool open(std::string filename, vec r, unsigned int h = 0){ + bool open(std::string filename, vec r, unsigned int h = 0){ for(unsigned int i = 0; i < D; i++) //set the dimensions of the binary file object R[i] = r[i]; @@ -117,7 +117,7 @@ public: /// @param filename is the name of the binary file to be created /// @param r is a STIM vector specifying the size of the file along each dimension /// @offset specifies how many bytes to offset the file (used to leave room for a header) - bool create(std::string filename, vec r, unsigned int offset = 0){ + bool create(std::string filename, vec r, unsigned int offset = 0){ std::ofstream target(filename.c_str(), std::ios::binary); diff --git a/stim/math/mathvec.h b/stim/math/mathvec.h deleted file mode 100644 index ae2ff3a..0000000 --- a/stim/math/mathvec.h +++ /dev/null @@ -1,376 +0,0 @@ -#ifndef RTS_VECTOR_H -#define RTS_VECTOR_H - -#include -#include -#include -#include -#include "../cuda/callable.h" - -namespace stim -{ - - - -template -struct vec : public std::vector -{ - using std::vector::size; - using std::vector::at; - using std::vector::resize; - using std::vector::push_back; - - vec(){ - - } - - /// Create a vector with a set dimension d - vec(int d) - { - resize(d,0); - } - - -// //efficiency constructors, makes construction easier for 1D-4D vectors - vec(T x, T y) - { - resize(2, 0); - at(0) = x; - at(1) = y; - } - vec(T x, T y, T z) - { - resize(3, 0); - at(0) = x; - at(1) = y; - at(2) = z; - } - vec(T x, T y, T z, T w) - { - resize(4, 0); - at(0) = x; - at(1) = y; - at(2) = z; - at(3) = w; - } - - - - //copy constructor - vec( const vec& other){ - unsigned int N = other.size(); - for(int i=0; i push(T x) - { - push_back(x); - return *this; - } - - vec push(T x, T y) - { - push_back(x); - push_back(y); - return *this; - } - vec push(T x, T y, T z) - { - push_back(x); - push_back(y); - push_back(z); - return *this; - } - vec push(T x, T y, T z, T w) - { - push_back(x); - push_back(y); - push_back(z); - push_back(w); - return *this; - } - - /// Casting operator. Creates a new vector with a new type U. - template< typename U > - operator vec(){ - unsigned int N = size(); - vec result; - for(int i=0; i r, theta, phi where theta = [0, 2*pi]) - vec cart2sph() const - { - - - vec sph; - sph.push_back(std::sqrt(at(0)*at(0) + at(1)*at(1) + at(2)*at(2))); - sph.push_back(std::atan2(at(1), at(0))); - - if(sph[0] == 0) - sph.push_back(0); - else - sph.push_back(std::acos(at(2) / sph[0])); - - return sph; - } - - /// Convert the vector from cartesian to spherical coordinates (r, theta, phi -> x, y, z where theta = [0, 2*pi]) - vec sph2cart() const - { - vec cart; - cart.push_back(at(0) * std::cos(at(1)) * std::sin(at(2))); - cart.push_back(at(0) * std::sin(at(1)) * std::sin(at(2))); - cart.push_back(at(0) * std::cos(at(2))); - - return cart; - } - - /// Computes the normalized vector (where each coordinate is divided by the L2 norm) - vec norm() const - { - unsigned int N = size(); - - //compute and return the unit vector - vec result; - - //compute the vector length - T l = len(); - - //normalize - for(int i=0; i cross(const vec rhs) const - { - - vec result(3); - - //compute the cross product (only valid for 3D vectors) - result[0] = (at(1) * rhs[2] - at(2) * rhs[1]); - result[1] = (at(2) * rhs[0] - at(0) * rhs[2]); - result[2] = (at(0) * rhs[1] - at(1) * rhs[0]); - - return result; - } - - /// Compute the Euclidean inner (dot) product - T dot(vec rhs) const - { - T result = (T)0; - unsigned int N = size(); - for(int i=0; i operator+(vec rhs) const - { - unsigned int N = size(); - vec result(N); - - for(int i=0; i operator+(T rhs) const - { - unsigned int N = size(); - - vec result(N); - for(int i=0; i operator-(vec rhs) const - { - unsigned int N = size(); - - vec result(N); - - for(int i=0; i operator-(T rhs) const - { - unsigned int N = size(); - - vec result(N); - for(int i=0; i operator*(T rhs) const - { - unsigned int N = size(); - - vec result(N); - - for(int i=0; i operator/(T rhs) const - { - unsigned int N = size(); - - vec result(N); - - for(int i=0; i operator*=(T rhs){ - - unsigned int N = size(); - for(int i=0; i operator+=(vec rhs){ - unsigned int N = size(); - for(int i=0; i & operator=(T rhs){ - - unsigned int N = size(); - for(int i=0; i - vec & operator=(vec rhs){ - unsigned int N = rhs.size(); - resize(N); - - for(int i=0; i operator-() const{ - - unsigned int N = size(); - - vec r(N); - - //negate the vector - for(int i=0; i -std::ostream& operator<<(std::ostream& os, stim::vec v) -{ - os< -stim::vec operator*(T lhs, stim::vec rhs) -{ - stim::vec r; - - return rhs * lhs; -} - -#endif diff --git a/stim/math/vector.h b/stim/math/vector.h new file mode 100644 index 0000000..ae2ff3a --- /dev/null +++ b/stim/math/vector.h @@ -0,0 +1,376 @@ +#ifndef RTS_VECTOR_H +#define RTS_VECTOR_H + +#include +#include +#include +#include +#include "../cuda/callable.h" + +namespace stim +{ + + + +template +struct vec : public std::vector +{ + using std::vector::size; + using std::vector::at; + using std::vector::resize; + using std::vector::push_back; + + vec(){ + + } + + /// Create a vector with a set dimension d + vec(int d) + { + resize(d,0); + } + + +// //efficiency constructors, makes construction easier for 1D-4D vectors + vec(T x, T y) + { + resize(2, 0); + at(0) = x; + at(1) = y; + } + vec(T x, T y, T z) + { + resize(3, 0); + at(0) = x; + at(1) = y; + at(2) = z; + } + vec(T x, T y, T z, T w) + { + resize(4, 0); + at(0) = x; + at(1) = y; + at(2) = z; + at(3) = w; + } + + + + //copy constructor + vec( const vec& other){ + unsigned int N = other.size(); + for(int i=0; i push(T x) + { + push_back(x); + return *this; + } + + vec push(T x, T y) + { + push_back(x); + push_back(y); + return *this; + } + vec push(T x, T y, T z) + { + push_back(x); + push_back(y); + push_back(z); + return *this; + } + vec push(T x, T y, T z, T w) + { + push_back(x); + push_back(y); + push_back(z); + push_back(w); + return *this; + } + + /// Casting operator. Creates a new vector with a new type U. + template< typename U > + operator vec(){ + unsigned int N = size(); + vec result; + for(int i=0; i r, theta, phi where theta = [0, 2*pi]) + vec cart2sph() const + { + + + vec sph; + sph.push_back(std::sqrt(at(0)*at(0) + at(1)*at(1) + at(2)*at(2))); + sph.push_back(std::atan2(at(1), at(0))); + + if(sph[0] == 0) + sph.push_back(0); + else + sph.push_back(std::acos(at(2) / sph[0])); + + return sph; + } + + /// Convert the vector from cartesian to spherical coordinates (r, theta, phi -> x, y, z where theta = [0, 2*pi]) + vec sph2cart() const + { + vec cart; + cart.push_back(at(0) * std::cos(at(1)) * std::sin(at(2))); + cart.push_back(at(0) * std::sin(at(1)) * std::sin(at(2))); + cart.push_back(at(0) * std::cos(at(2))); + + return cart; + } + + /// Computes the normalized vector (where each coordinate is divided by the L2 norm) + vec norm() const + { + unsigned int N = size(); + + //compute and return the unit vector + vec result; + + //compute the vector length + T l = len(); + + //normalize + for(int i=0; i cross(const vec rhs) const + { + + vec result(3); + + //compute the cross product (only valid for 3D vectors) + result[0] = (at(1) * rhs[2] - at(2) * rhs[1]); + result[1] = (at(2) * rhs[0] - at(0) * rhs[2]); + result[2] = (at(0) * rhs[1] - at(1) * rhs[0]); + + return result; + } + + /// Compute the Euclidean inner (dot) product + T dot(vec rhs) const + { + T result = (T)0; + unsigned int N = size(); + for(int i=0; i operator+(vec rhs) const + { + unsigned int N = size(); + vec result(N); + + for(int i=0; i operator+(T rhs) const + { + unsigned int N = size(); + + vec result(N); + for(int i=0; i operator-(vec rhs) const + { + unsigned int N = size(); + + vec result(N); + + for(int i=0; i operator-(T rhs) const + { + unsigned int N = size(); + + vec result(N); + for(int i=0; i operator*(T rhs) const + { + unsigned int N = size(); + + vec result(N); + + for(int i=0; i operator/(T rhs) const + { + unsigned int N = size(); + + vec result(N); + + for(int i=0; i operator*=(T rhs){ + + unsigned int N = size(); + for(int i=0; i operator+=(vec rhs){ + unsigned int N = size(); + for(int i=0; i & operator=(T rhs){ + + unsigned int N = size(); + for(int i=0; i + vec & operator=(vec rhs){ + unsigned int N = rhs.size(); + resize(N); + + for(int i=0; i operator-() const{ + + unsigned int N = size(); + + vec r(N); + + //negate the vector + for(int i=0; i +std::ostream& operator<<(std::ostream& os, stim::vec v) +{ + os< +stim::vec operator*(T lhs, stim::vec rhs) +{ + stim::vec r; + + return rhs * lhs; +} + +#endif diff --git a/stim/visualization/colormap.h b/stim/visualization/colormap.h index 01077fd..14fbde9 100644 --- a/stim/visualization/colormap.h +++ b/stim/visualization/colormap.h @@ -287,7 +287,7 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T } template -static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale)//, bool positive = false) +static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale) { //computes the max and min range automatically @@ -329,13 +329,13 @@ static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, u } template -static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, colormapType cm = cmGrayscale, bool positive = false) +static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, colormapType cm = cmGrayscale) { //allocate a color buffer unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size); //do the mapping - cpu2cpu(cpuSource, cpuBuffer, x_size * y_size, cm, positive); + cpu2cpu(cpuSource, cpuBuffer, x_size * y_size, cm); //copy the buffer to an image buffer2image(cpuBuffer, fileDest, x_size, y_size); -- libgit2 0.21.4