From 81e0d2215b30a17309ce7ba7615f94e2eda5c67e Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Fri, 5 Dec 2014 17:53:46 -0600 Subject: [PATCH] separated executable arguments and options in the arglist class --- math/complex.h | 1010 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- math/complexfield.cuh | 274 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------- math/field.cuh | 686 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- math/rect.h | 438 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- math/triangle.h | 376 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- optics/material.h | 306 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------- optics/mirst-1d.cuh | 894 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ui/arguments.h | 111 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------- 8 files changed, 2056 insertions(+), 2039 deletions(-) diff --git a/math/complex.h b/math/complex.h index 94df151..84ba67e 100644 --- a/math/complex.h +++ b/math/complex.h @@ -1,505 +1,505 @@ -/*RTS Complex number class. This class is CUDA compatible, -and can therefore be used in CUDA code and on CUDA devices. -*/ - -#ifndef RTS_COMPLEX -#define RTS_COMPLEX - -#include "../cuda/callable.h" -#include -#include -#include -#include - -namespace stim -{ - -template -struct complex -{ - T r, i; - - //default constructor - CUDA_CALLABLE complex() - { - r = 0; - i = 0; - } - - //constructor when given real and imaginary values - CUDA_CALLABLE complex(T r, T i = 0) - { - this->r = r; - this->i = i; - } - - //access methods - CUDA_CALLABLE T real() - { - return r; - } - - CUDA_CALLABLE T real(T r_val) - { - r = r_val; - return r_val; - } - - CUDA_CALLABLE T imag() - { - return i; - } - CUDA_CALLABLE T imag(T i_val) - { - i = i_val; - return i_val; - } - - - - //return the current value multiplied by i - CUDA_CALLABLE complex imul() - { - complex result; - result.r = -i; - result.i = r; - - return result; - } - - //returns the complex signum (-1, 0, 1) - CUDA_CALLABLE int sgn(){ - if(r > 0) return 1; - else if(r < 0) return -1; - else return (0 < i - i < 0); - } - - //ARITHMETIC OPERATORS-------------------- - - //binary + operator (returns the result of adding two complex values) - CUDA_CALLABLE complex operator+ (const complex rhs) const - { - complex result; - result.r = r + rhs.r; - result.i = i + rhs.i; - return result; - } - - CUDA_CALLABLE complex operator+ (const T rhs) const - { - complex result; - result.r = r + rhs; - result.i = i; - return result; - } - - //binary - operator (returns the result of adding two complex values) - CUDA_CALLABLE complex operator- (const complex rhs) const - { - complex result; - result.r = r - rhs.r; - result.i = i - rhs.i; - return result; - } - - //binary - operator (returns the result of adding two complex values) - CUDA_CALLABLE complex operator- (const T rhs) - { - complex result; - result.r = r - rhs; - result.i = i; - return result; - } - - //binary MULTIPLICATION operators (returns the result of multiplying complex values) - CUDA_CALLABLE complex operator* (const complex rhs) const - { - complex result; - result.r = r * rhs.r - i * rhs.i; - result.i = r * rhs.i + i * rhs.r; - return result; - } - CUDA_CALLABLE complex operator* (const T rhs) - { - return complex(r * rhs, i * rhs); - } - - //binary DIVISION operators (returns the result of dividing complex values) - CUDA_CALLABLE complex operator/ (const complex rhs) const - { - complex result; - T denom = rhs.r * rhs.r + rhs.i * rhs.i; - result.r = (r * rhs.r + i * rhs.i) / denom; - result.i = (- r * rhs.i + i * rhs.r) / denom; - - return result; - } - CUDA_CALLABLE complex operator/ (const T rhs) - { - return complex(r / rhs, i / rhs); - } - - //ASSIGNMENT operators----------------------------------- - CUDA_CALLABLE complex & operator=(const complex &rhs) - { - //check for self-assignment - if(this != &rhs) - { - this->r = rhs.r; - this->i = rhs.i; - } - return *this; - } - CUDA_CALLABLE complex & operator=(const T &rhs) - { - this->r = rhs; - this->i = 0; - - return *this; - } - - //arithmetic assignment operators - CUDA_CALLABLE complex operator+=(const complex &rhs) - { - *this = *this + rhs; - return *this; - } - CUDA_CALLABLE complex operator+=(const T &rhs) - { - *this = *this + rhs; - return *this; - } - - CUDA_CALLABLE complex operator-=(const complex &rhs) - { - *this = *this - rhs; - return *this; - } - CUDA_CALLABLE complex operator-=(const T &rhs) - { - *this = *this - rhs; - return *this; - } - - CUDA_CALLABLE complex operator*=(const complex &rhs) - { - *this = *this * rhs; - return *this; - } - CUDA_CALLABLE complex operator*=(const T &rhs) - { - *this = *this * rhs; - return *this; - } - //divide and assign - CUDA_CALLABLE complex operator/=(const complex &rhs) - { - *this = *this / rhs; - return *this; - } - CUDA_CALLABLE complex operator/=(const T &rhs) - { - *this = *this / rhs; - return *this; - } - - //absolute value operator (returns the absolute value of the complex number) - CUDA_CALLABLE T abs() - { - return std::sqrt(r * r + i * i); - } - - CUDA_CALLABLE complex log() - { - complex result; - result.r = (T)std::log(std::sqrt(r * r + i * i)); - result.i = (T)std::atan2(i, r); - - - return result; - } - - CUDA_CALLABLE complex exp() - { - complex result; - - T e_r = std::exp(r); - result.r = e_r * (T)std::cos(i); - result.i = e_r * (T)std::sin(i); - - return result; - } - - /*CUDA_CALLABLE complex pow(int y) - { - - return pow((double)y); - }*/ - - CUDA_CALLABLE complex pow(T y) - { - complex result; - - result = log() * y; - - return result.exp(); - } - - CUDA_CALLABLE complex sqrt() - { - complex result; - - //convert to polar coordinates - T a = std::sqrt(r*r + i*i); - T theta = std::atan2(i, r); - - //find the square root - T a_p = std::sqrt(a); - T theta_p = theta/2.0f; - - //convert back to cartesian coordinates - result.r = a_p * std::cos(theta_p); - result.i = a_p * std::sin(theta_p); - - return result; - } - - std::string str() - { - std::stringstream ss; - ss<<"("< rhs) - { - if(r == rhs.r && i == rhs.i) - return true; - return false; - } - - CUDA_CALLABLE bool operator==(T rhs) - { - if(r == rhs && i == 0) - return true; - return false; - } - - CUDA_CALLABLE bool operator!=(T rhs) - { - if(r != rhs || i != 0) - return true; - return false; - } - - CUDA_CALLABLE bool operator<(complex rhs){ - return abs() < rhs.abs(); - } - CUDA_CALLABLE bool operator<=(complex rhs){ - return abs() <= rhs.abs(); - } - CUDA_CALLABLE bool operator>(complex rhs){ - return abs() > rhs.abs(); - } - CUDA_CALLABLE bool operator >=(complex rhs){ - return abs() >= rhs.abs(); - } - - //CASTING operators - template < typename otherT > - operator complex() - { - complex result((otherT)r, (otherT)i); - return result; - } - template< typename otherT > - complex( const complex &rhs) - { - r = (T)rhs.r; - i = (T)rhs.i; - } - template< typename otherT > - complex& operator=(const complex &rhs) - { - r = (T)rhs.r; - i = (T)rhs.i; - return *this; - } - -}; - -} //end RTS namespace - -//addition -template -CUDA_CALLABLE static stim::complex operator+(const double a, const stim::complex b) -{ - return stim::complex((T)a + b.r, b.i); -} - -//subtraction with a real value -template -CUDA_CALLABLE static stim::complex operator-(const double a, const stim::complex b) -{ - return stim::complex((T)a - b.r, -b.i); -} - -//minus sign -template -CUDA_CALLABLE static stim::complex operator-(const stim::complex &rhs) -{ - return stim::complex(-rhs.r, -rhs.i); -} - -//multiply a T value by a complex value -template -CUDA_CALLABLE static stim::complex operator*(const double a, const stim::complex b) -{ - return stim::complex((T)a * b.r, (T)a * b.i); -} - -//divide a T value by a complex value -template -CUDA_CALLABLE static stim::complex operator/(const double a, const stim::complex b) -{ - stim::complex result; - - T denom = b.r * b.r + b.i * b.i; - - result.r = ((T)a * b.r) / denom; - result.i = -((T)a * b.i) / denom; - - return result; -} - - -template -CUDA_CALLABLE static stim::complex pow(stim::complex x, T y) -{ - return x.pow(y); -} -template -CUDA_CALLABLE static stim::complex pow(stim::complex x, int y) -{ - return x.pow(y); -} - -//log function -template -CUDA_CALLABLE static stim::complex log(stim::complex x) -{ - return x.log(); -} - -//exp function -template -CUDA_CALLABLE static stim::complex exp(stim::complex x) -{ - return x.exp(); -} - -//sqrt function -template -CUDA_CALLABLE static stim::complex sqrt(stim::complex x) -{ - return x.sqrt(); -} - - -template -CUDA_CALLABLE static T abs(stim::complex a) -{ - return a.abs(); -} - -template -CUDA_CALLABLE static T real(stim::complex a) -{ - return a.r; -} - -//template -CUDA_CALLABLE static float real(float a) -{ - return a; -} - -template -CUDA_CALLABLE static T imag(stim::complex a) -{ - return a.i; -} - -//trigonometric functions -//template -/*CUDA_CALLABLE static stim::complex sinf(const stim::complex x) -{ - stim::complex result; - result.r = sinf(x.r) * coshf(x.i); - result.i = cosf(x.r) * sinhf(x.i); - - return result; -}*/ - -template -CUDA_CALLABLE stim::complex sin(const stim::complex x) -{ - stim::complex result; - result.r = (A)std::sin(x.r) * (A)std::cosh(x.i); - result.i = (A)std::cos(x.r) * (A)std::sinh(x.i); - - return result; -} - -//floating point template -//template -/*CUDA_CALLABLE static stim::complex cosf(const stim::complex x) -{ - stim::complex result; - result.r = cosf(x.r) * coshf(x.i); - result.i = -(sinf(x.r) * sinhf(x.i)); - - return result; -}*/ - -template -CUDA_CALLABLE stim::complex cos(const stim::complex x) -{ - stim::complex result; - result.r = (A)std::cos(x.r) * (A)std::cosh(x.i); - result.i = -((A)std::sin(x.r) * (A)std::sinh(x.i)); - - return result; -} - - -template -std::ostream& operator<<(std::ostream& os, stim::complex x) -{ - os< -std::istream& operator>>(std::istream& is, stim::complex& x) -{ - A r, i; - r = i = 0; //initialize the real and imaginary parts to zero - is>>r; //parse - is>>i; - - x.real(r); //assign the parsed values to x - x.imag(i); - - return is; //return the stream -} - -//#if __GNUC__ > 3 && __GNUC_MINOR__ > 7 -//template using rtsComplex = stim::complex; -//#endif - - - -#endif +/*RTS Complex number class. This class is CUDA compatible, +and can therefore be used in CUDA code and on CUDA devices. +*/ + +#ifndef RTS_COMPLEX +#define RTS_COMPLEX + +#include "../cuda/callable.h" +#include +#include +#include +#include + +namespace stim +{ + +template +struct complex +{ + T r, i; + + //default constructor + CUDA_CALLABLE complex() + { + r = 0; + i = 0; + } + + //constructor when given real and imaginary values + CUDA_CALLABLE complex(T r, T i = 0) + { + this->r = r; + this->i = i; + } + + //access methods + CUDA_CALLABLE T real() + { + return r; + } + + CUDA_CALLABLE T real(T r_val) + { + r = r_val; + return r_val; + } + + CUDA_CALLABLE T imag() + { + return i; + } + CUDA_CALLABLE T imag(T i_val) + { + i = i_val; + return i_val; + } + + + + //return the current value multiplied by i + CUDA_CALLABLE complex imul() + { + complex result; + result.r = -i; + result.i = r; + + return result; + } + + //returns the complex signum (-1, 0, 1) + CUDA_CALLABLE int sgn(){ + if(r > 0) return 1; + else if(r < 0) return -1; + else return (0 < i - i < 0); + } + + //ARITHMETIC OPERATORS-------------------- + + //binary + operator (returns the result of adding two complex values) + CUDA_CALLABLE complex operator+ (const complex rhs) const + { + complex result; + result.r = r + rhs.r; + result.i = i + rhs.i; + return result; + } + + CUDA_CALLABLE complex operator+ (const T rhs) const + { + complex result; + result.r = r + rhs; + result.i = i; + return result; + } + + //binary - operator (returns the result of adding two complex values) + CUDA_CALLABLE complex operator- (const complex rhs) const + { + complex result; + result.r = r - rhs.r; + result.i = i - rhs.i; + return result; + } + + //binary - operator (returns the result of adding two complex values) + CUDA_CALLABLE complex operator- (const T rhs) + { + complex result; + result.r = r - rhs; + result.i = i; + return result; + } + + //binary MULTIPLICATION operators (returns the result of multiplying complex values) + CUDA_CALLABLE complex operator* (const complex rhs) const + { + complex result; + result.r = r * rhs.r - i * rhs.i; + result.i = r * rhs.i + i * rhs.r; + return result; + } + CUDA_CALLABLE complex operator* (const T rhs) + { + return complex(r * rhs, i * rhs); + } + + //binary DIVISION operators (returns the result of dividing complex values) + CUDA_CALLABLE complex operator/ (const complex rhs) const + { + complex result; + T denom = rhs.r * rhs.r + rhs.i * rhs.i; + result.r = (r * rhs.r + i * rhs.i) / denom; + result.i = (- r * rhs.i + i * rhs.r) / denom; + + return result; + } + CUDA_CALLABLE complex operator/ (const T rhs) + { + return complex(r / rhs, i / rhs); + } + + //ASSIGNMENT operators----------------------------------- + CUDA_CALLABLE complex & operator=(const complex &rhs) + { + //check for self-assignment + if(this != &rhs) + { + this->r = rhs.r; + this->i = rhs.i; + } + return *this; + } + CUDA_CALLABLE complex & operator=(const T &rhs) + { + this->r = rhs; + this->i = 0; + + return *this; + } + + //arithmetic assignment operators + CUDA_CALLABLE complex operator+=(const complex &rhs) + { + *this = *this + rhs; + return *this; + } + CUDA_CALLABLE complex operator+=(const T &rhs) + { + *this = *this + rhs; + return *this; + } + + CUDA_CALLABLE complex operator-=(const complex &rhs) + { + *this = *this - rhs; + return *this; + } + CUDA_CALLABLE complex operator-=(const T &rhs) + { + *this = *this - rhs; + return *this; + } + + CUDA_CALLABLE complex operator*=(const complex &rhs) + { + *this = *this * rhs; + return *this; + } + CUDA_CALLABLE complex operator*=(const T &rhs) + { + *this = *this * rhs; + return *this; + } + //divide and assign + CUDA_CALLABLE complex operator/=(const complex &rhs) + { + *this = *this / rhs; + return *this; + } + CUDA_CALLABLE complex operator/=(const T &rhs) + { + *this = *this / rhs; + return *this; + } + + //absolute value operator (returns the absolute value of the complex number) + CUDA_CALLABLE T abs() + { + return std::sqrt(r * r + i * i); + } + + CUDA_CALLABLE complex log() + { + complex result; + result.r = (T)std::log(std::sqrt(r * r + i * i)); + result.i = (T)std::atan2(i, r); + + + return result; + } + + CUDA_CALLABLE complex exp() + { + complex result; + + T e_r = std::exp(r); + result.r = e_r * (T)std::cos(i); + result.i = e_r * (T)std::sin(i); + + return result; + } + + /*CUDA_CALLABLE complex pow(int y) + { + + return pow((double)y); + }*/ + + CUDA_CALLABLE complex pow(T y) + { + complex result; + + result = log() * y; + + return result.exp(); + } + + CUDA_CALLABLE complex sqrt() + { + complex result; + + //convert to polar coordinates + T a = std::sqrt(r*r + i*i); + T theta = std::atan2(i, r); + + //find the square root + T a_p = std::sqrt(a); + T theta_p = theta/2.0f; + + //convert back to cartesian coordinates + result.r = a_p * std::cos(theta_p); + result.i = a_p * std::sin(theta_p); + + return result; + } + + std::string str() + { + std::stringstream ss; + ss<<"("< rhs) + { + if(r == rhs.r && i == rhs.i) + return true; + return false; + } + + CUDA_CALLABLE bool operator==(T rhs) + { + if(r == rhs && i == 0) + return true; + return false; + } + + CUDA_CALLABLE bool operator!=(T rhs) + { + if(r != rhs || i != 0) + return true; + return false; + } + + CUDA_CALLABLE bool operator<(complex rhs){ + return abs() < rhs.abs(); + } + CUDA_CALLABLE bool operator<=(complex rhs){ + return abs() <= rhs.abs(); + } + CUDA_CALLABLE bool operator>(complex rhs){ + return abs() > rhs.abs(); + } + CUDA_CALLABLE bool operator >=(complex rhs){ + return abs() >= rhs.abs(); + } + + //CASTING operators + template < typename otherT > + operator complex() + { + complex result((otherT)r, (otherT)i); + return result; + } + template< typename otherT > + complex( const complex &rhs) + { + r = (T)rhs.r; + i = (T)rhs.i; + } + template< typename otherT > + complex& operator=(const complex &rhs) + { + r = (T)rhs.r; + i = (T)rhs.i; + return *this; + } + +}; + +} //end RTS namespace + +//addition +template +CUDA_CALLABLE static stim::complex operator+(const double a, const stim::complex b) +{ + return stim::complex((T)a + b.r, b.i); +} + +//subtraction with a real value +template +CUDA_CALLABLE static stim::complex operator-(const double a, const stim::complex b) +{ + return stim::complex((T)a - b.r, -b.i); +} + +//minus sign +template +CUDA_CALLABLE static stim::complex operator-(const stim::complex &rhs) +{ + return stim::complex(-rhs.r, -rhs.i); +} + +//multiply a T value by a complex value +template +CUDA_CALLABLE static stim::complex operator*(const double a, const stim::complex b) +{ + return stim::complex((T)a * b.r, (T)a * b.i); +} + +//divide a T value by a complex value +template +CUDA_CALLABLE static stim::complex operator/(const double a, const stim::complex b) +{ + stim::complex result; + + T denom = b.r * b.r + b.i * b.i; + + result.r = ((T)a * b.r) / denom; + result.i = -((T)a * b.i) / denom; + + return result; +} + + +template +CUDA_CALLABLE static stim::complex pow(stim::complex x, T y) +{ + return x.pow(y); +} +template +CUDA_CALLABLE static stim::complex pow(stim::complex x, int y) +{ + return x.pow(y); +} + +//log function +template +CUDA_CALLABLE static stim::complex log(stim::complex x) +{ + return x.log(); +} + +//exp function +template +CUDA_CALLABLE static stim::complex exp(stim::complex x) +{ + return x.exp(); +} + +//sqrt function +template +CUDA_CALLABLE static stim::complex sqrt(stim::complex x) +{ + return x.sqrt(); +} + + +template +CUDA_CALLABLE static T abs(stim::complex a) +{ + return a.abs(); +} + +template +CUDA_CALLABLE static T real(stim::complex a) +{ + return a.r; +} + +//template +CUDA_CALLABLE static float real(float a) +{ + return a; +} + +template +CUDA_CALLABLE static T imag(stim::complex a) +{ + return a.i; +} + +//trigonometric functions +//template +/*CUDA_CALLABLE static stim::complex sinf(const stim::complex x) +{ + stim::complex result; + result.r = sinf(x.r) * coshf(x.i); + result.i = cosf(x.r) * sinhf(x.i); + + return result; +}*/ + +template +CUDA_CALLABLE stim::complex sin(const stim::complex x) +{ + stim::complex result; + result.r = (A)std::sin(x.r) * (A)std::cosh(x.i); + result.i = (A)std::cos(x.r) * (A)std::sinh(x.i); + + return result; +} + +//floating point template +//template +/*CUDA_CALLABLE static stim::complex cosf(const stim::complex x) +{ + stim::complex result; + result.r = cosf(x.r) * coshf(x.i); + result.i = -(sinf(x.r) * sinhf(x.i)); + + return result; +}*/ + +template +CUDA_CALLABLE stim::complex cos(const stim::complex x) +{ + stim::complex result; + result.r = (A)std::cos(x.r) * (A)std::cosh(x.i); + result.i = -((A)std::sin(x.r) * (A)std::sinh(x.i)); + + return result; +} + + +template +std::ostream& operator<<(std::ostream& os, stim::complex x) +{ + os< +std::istream& operator>>(std::istream& is, stim::complex& x) +{ + A r, i; + r = i = 0; //initialize the real and imaginary parts to zero + is>>r; //parse + is>>i; + + x.real(r); //assign the parsed values to x + x.imag(i); + + return is; //return the stream +} + +//#if __GNUC__ > 3 && __GNUC_MINOR__ > 7 +//template using rtsComplex = stim::complex; +//#endif + + + +#endif diff --git a/math/complexfield.cuh b/math/complexfield.cuh index 701e97e..81488cf 100644 --- a/math/complexfield.cuh +++ b/math/complexfield.cuh @@ -1,137 +1,137 @@ -#ifndef RTS_COMPLEXFIELD_H -#define RTS_COMPLEXFIELD_H - -#include "cublas_v2.h" -#include - -#include "../math/field.cuh" -#include "../math/complex.h" -#include "../math/realfield.cuh" - -namespace stim{ - -template -__global__ void gpu_complexfield_mag(T* dest, complex* source, unsigned int r0, unsigned int r1){ - - int iu = blockIdx.x * blockDim.x + threadIdx.x; - int iv = blockIdx.y * blockDim.y + threadIdx.y; - - //make sure that the thread indices are in-bounds - if(iu >= r0 || iv >= r1) return; - - //compute the index into the field - int i = iv*r0 + iu; - - //calculate and store the result - dest[i] = source[i].abs(); -} - -/*This class stores functions for saving images of complex fields -*/ -template -class complexfield : public field< stim::complex, D >{ - using field< stim::complex, D >::R; - using field< stim::complex, D >::X; - using field< stim::complex, D >::shape; - using field< stim::complex, D >::cuda_params; - - - -public: - - //find the maximum value of component n - stim::complex find_max(unsigned int n){ - cublasStatus_t stat; - cublasHandle_t handle; - - //create a CUBLAS handle - stat = cublasCreate(&handle); - if(stat != CUBLAS_STATUS_SUCCESS){ - std::cout<<"CUBLAS Error: initialization failed"< result; - - if(sizeof(T) == 8) - stat = cublasIcamax(handle, L, (const cuComplex*)X[n], 1, &index); - else - stat = cublasIzamax(handle, L, (const cuDoubleComplex*)X[n], 1, &index); - - index -= 1; //adjust for 1-based indexing - - //if there was a GPU error, terminate - if(stat != CUBLAS_STATUS_SUCCESS){ - std::cout<<"CUBLAS Error: failure finding maximum value."<), cudaMemcpyDeviceToHost)); - return result; - } - -public: - - enum attribute {magnitude, real, imaginary}; - - //constructor (no parameters) - complexfield() : field, D>(){}; - - //constructor (resolution specified) - complexfield(unsigned int r0, unsigned int r1) : field, D>(r0, r1){}; - - //assignment from a field of complex values - complexfield & operator=(const field< stim::complex, D > rhs){ - field< complex, D >::operator=(rhs); - return *this; - } - - //assignment operator (scalar value) - complexfield & operator= (const complex rhs){ - - field< complex, D >::operator=(rhs); - return *this; - } - - //assignment operator (vector value) - complexfield & operator= (const vec< complex, D > rhs){ - - field< complex, D >::operator=(rhs); - return *this; - } - - //cropping - complexfield crop(unsigned int width, unsigned int height){ - - complexfield result; - result = field< complex, D>::crop(width, height); - return result; - } - - void toImage(std::string filename, attribute type = magnitude, unsigned int n=0){ - - field rf(R[0], R[1]); - - //get cuda parameters - dim3 blocks, grids; - cuda_params(grids, blocks); - - if(type == magnitude){ - gpu_complexfield_mag <<>> (rf.ptr(), X[n], R[0], R[1]); - rf.toImage(filename, n, true); - } - - } - - -}; - - -} //end namespace rts - - -#endif +#ifndef RTS_COMPLEXFIELD_H +#define RTS_COMPLEXFIELD_H + +#include "cublas_v2.h" +#include + +#include "../math/field.cuh" +#include "../math/complex.h" +#include "../math/realfield.cuh" + +namespace stim{ + +template +__global__ void gpu_complexfield_mag(T* dest, complex* source, unsigned int r0, unsigned int r1){ + + int iu = blockIdx.x * blockDim.x + threadIdx.x; + int iv = blockIdx.y * blockDim.y + threadIdx.y; + + //make sure that the thread indices are in-bounds + if(iu >= r0 || iv >= r1) return; + + //compute the index into the field + int i = iv*r0 + iu; + + //calculate and store the result + dest[i] = source[i].abs(); +} + +/*This class stores functions for saving images of complex fields +*/ +template +class complexfield : public field< stim::complex, D >{ + using field< stim::complex, D >::R; + using field< stim::complex, D >::X; + using field< stim::complex, D >::shape; + using field< stim::complex, D >::cuda_params; + + + +public: + + //find the maximum value of component n + stim::complex find_max(unsigned int n){ + cublasStatus_t stat; + cublasHandle_t handle; + + //create a CUBLAS handle + stat = cublasCreate(&handle); + if(stat != CUBLAS_STATUS_SUCCESS){ + std::cout<<"CUBLAS Error: initialization failed"< result; + + if(sizeof(T) == 8) + stat = cublasIcamax(handle, L, (const cuComplex*)X[n], 1, &index); + else + stat = cublasIzamax(handle, L, (const cuDoubleComplex*)X[n], 1, &index); + + index -= 1; //adjust for 1-based indexing + + //if there was a GPU error, terminate + if(stat != CUBLAS_STATUS_SUCCESS){ + std::cout<<"CUBLAS Error: failure finding maximum value."<), cudaMemcpyDeviceToHost)); + return result; + } + +public: + + enum attribute {magnitude, real, imaginary}; + + //constructor (no parameters) + complexfield() : field, D>(){}; + + //constructor (resolution specified) + complexfield(unsigned int r0, unsigned int r1) : field, D>(r0, r1){}; + + //assignment from a field of complex values + complexfield & operator=(const field< stim::complex, D > rhs){ + field< complex, D >::operator=(rhs); + return *this; + } + + //assignment operator (scalar value) + complexfield & operator= (const complex rhs){ + + field< complex, D >::operator=(rhs); + return *this; + } + + //assignment operator (vector value) + complexfield & operator= (const vec< complex, D > rhs){ + + field< complex, D >::operator=(rhs); + return *this; + } + + //cropping + complexfield crop(unsigned int width, unsigned int height){ + + complexfield result; + result = field< complex, D>::crop(width, height); + return result; + } + + void toImage(std::string filename, attribute type = magnitude, unsigned int n=0){ + + field rf(R[0], R[1]); + + //get cuda parameters + dim3 blocks, grids; + cuda_params(grids, blocks); + + if(type == magnitude){ + gpu_complexfield_mag <<>> (rf.ptr(), X[n], R[0], R[1]); + rf.toImage(filename, n, true); + } + + } + + +}; + + +} //end namespace rts + + +#endif diff --git a/math/field.cuh b/math/field.cuh index aa6b10c..1a3a2ab 100644 --- a/math/field.cuh +++ b/math/field.cuh @@ -1,343 +1,343 @@ -#ifndef RTS_FIELD_CUH -#define RTS_FIELD_CUH - -#include -#include -#include - -#include "cublas_v2.h" -#include - -#include "../math/rect.h" -#include "../cuda/threads.h" -#include "../cuda/error.h" -#include "../cuda/devices.h" -#include "../visualization/colormap.h" - - -namespace stim{ - -//multiply R = X * Y -template -__global__ void gpu_field_multiply(T* R, T* X, T* Y, unsigned int r0, unsigned int r1){ - - int iu = blockIdx.x * blockDim.x + threadIdx.x; - int iv = blockIdx.y * blockDim.y + threadIdx.y; - - //make sure that the thread indices are in-bounds - if(iu >= r0 || iv >= r1) return; - - //compute the index into the field - int i = iv*r0 + iu; - - //calculate and store the result - R[i] = X[i] * Y[i]; -} - -//assign a constant value to all points -template -__global__ void gpu_field_assign(T* ptr, T val, unsigned int r0, unsigned int r1){ - - int iu = blockIdx.x * blockDim.x + threadIdx.x; - int iv = blockIdx.y * blockDim.y + threadIdx.y; - - //make sure that the thread indices are in-bounds - if(iu >= r0 || iv >= r1) return; - - //compute the index into the field - int i = iv*r0 + iu; - - //calculate and store the result - ptr[i] = val; -} - -//crop the field to the new dimensions (width x height) -template -__global__ void gpu_field_crop(T* dest, T* source, - unsigned int r0, unsigned int r1, - unsigned int width, unsigned int height){ - - int iu = blockIdx.x * blockDim.x + threadIdx.x; - int iv = blockIdx.y * blockDim.y + threadIdx.y; - - //make sure that the thread indices are in-bounds - if(iu >= width || iv >= height) return; - - //compute the index into the field - int is = iv*r0 + iu; - int id = iv*width + iu; - - //calculate and store the result - dest[id] = source[is]; -} - -template -class field{ - -protected: - - T* X[D]; //pointer to the field data - unsigned int R[2]; //field resolution - stim::rect shape; //position and shape of the field slice - - //calculates the optimal block and grid sizes using information from the GPU - void cuda_params(dim3& grids, dim3& blocks){ - int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size - int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); - - //create one thread for each detector pixel - blocks = dim3(SQRT_BLOCK, SQRT_BLOCK); - grids = dim3((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); - } - - //find the maximum value of component n - T find_max(unsigned int n){ - cublasStatus_t stat; - cublasHandle_t handle; - - //create a CUBLAS handle - stat = cublasCreate(&handle); - if(stat != CUBLAS_STATUS_SUCCESS){ - std::cout<<"CUBLAS Error: initialization failed"< process_filename(std::string name){ - std::stringstream ss(name); - std::string item; - std::vector elems; - while(std::getline(ss, item, '.')) //split the string at the '.' character (filename and extension) - { - elems.push_back(item); - } - - std::string prefix = elems[0]; //prefix contains the filename (with wildcard '?' characters) - std::string ext = elems[1]; //file extension (ex. .bmp, .png) - ext = std::string(".") + ext; //add a period back into the extension - - size_t i0 = prefix.find_first_of("?"); //find the positions of the first and last wildcard ('?'') - size_t i1 = prefix.find_last_of("?"); - - std::string postfix = prefix.substr(i1+1); - prefix = prefix.substr(0, i0); - - unsigned int digits = i1 - i0 + 1; //compute the number of wildcards - - std::vector flist; //create a vector of file names - //fill the list - for(unsigned int d=0; d>> (X[n], rhs, R[0], R[1]); - - return *this; - } - - //assignment of vector component - field & operator= (const vec rhs){ - - int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size - int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); - - //create one thread for each detector pixel - dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); - dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); - - //assign the constant value to all positions and dimensions - for(unsigned int n=0; n>> (X[n], rhs.v[n], R[0], R[1]); - - return *this; - - } - - //multiply two fields (element-wise multiplication) - field operator* (const field & rhs){ - - int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size - int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); - - //create one thread for each detector pixel - dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); - dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); - - //create a scalar field to store the result - field result(R[0], R[1]); - - for(int n=0; n>> (result.X[n], X[n], rhs.X[n], R[0], R[1]); - - return result; - } - - T* ptr(unsigned int n = 0){ - if(n < D) - return X[n]; - else return NULL; - } - - //return the vector component at position (u, v) - vec get(unsigned int u, unsigned int v){ - - vec result; - for(unsigned int d=0; d crop(unsigned int width, unsigned int height){ - int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size - int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); - - //create one thread for each detector pixel - dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); - dim3 dimGrid((width + SQRT_BLOCK -1)/SQRT_BLOCK, (height + SQRT_BLOCK - 1)/SQRT_BLOCK); - - //create a scalar field to store the result - field result(width, height); - - for(int n=0; n>> (result.X[n], X[n], R[0], R[1], width, height); - - return result; - } - - //save an image representing component n - void toImage(std::string filename, unsigned int n = 0, - bool positive = false, stim::colormapType cmap = stim::cmBrewer){ - T max_val = find_max(n); //find the maximum value - - if(positive) //if the field is positive, use the range [0 max_val] - stim::gpu2image(X[n], filename, R[0], R[1], 0, max_val, cmap); - else - stim::gpu2image(X[n], filename, R[0], R[1], -max_val, max_val, cmap); - } - -}; - -} //end namespace rts -#endif +#ifndef RTS_FIELD_CUH +#define RTS_FIELD_CUH + +#include +#include +#include + +#include "cublas_v2.h" +#include + +#include "../math/rect.h" +#include "../cuda/threads.h" +#include "../cuda/error.h" +#include "../cuda/devices.h" +#include "../visualization/colormap.h" + + +namespace stim{ + +//multiply R = X * Y +template +__global__ void gpu_field_multiply(T* R, T* X, T* Y, unsigned int r0, unsigned int r1){ + + int iu = blockIdx.x * blockDim.x + threadIdx.x; + int iv = blockIdx.y * blockDim.y + threadIdx.y; + + //make sure that the thread indices are in-bounds + if(iu >= r0 || iv >= r1) return; + + //compute the index into the field + int i = iv*r0 + iu; + + //calculate and store the result + R[i] = X[i] * Y[i]; +} + +//assign a constant value to all points +template +__global__ void gpu_field_assign(T* ptr, T val, unsigned int r0, unsigned int r1){ + + int iu = blockIdx.x * blockDim.x + threadIdx.x; + int iv = blockIdx.y * blockDim.y + threadIdx.y; + + //make sure that the thread indices are in-bounds + if(iu >= r0 || iv >= r1) return; + + //compute the index into the field + int i = iv*r0 + iu; + + //calculate and store the result + ptr[i] = val; +} + +//crop the field to the new dimensions (width x height) +template +__global__ void gpu_field_crop(T* dest, T* source, + unsigned int r0, unsigned int r1, + unsigned int width, unsigned int height){ + + int iu = blockIdx.x * blockDim.x + threadIdx.x; + int iv = blockIdx.y * blockDim.y + threadIdx.y; + + //make sure that the thread indices are in-bounds + if(iu >= width || iv >= height) return; + + //compute the index into the field + int is = iv*r0 + iu; + int id = iv*width + iu; + + //calculate and store the result + dest[id] = source[is]; +} + +template +class field{ + +protected: + + T* X[D]; //pointer to the field data + unsigned int R[2]; //field resolution + stim::rect shape; //position and shape of the field slice + + //calculates the optimal block and grid sizes using information from the GPU + void cuda_params(dim3& grids, dim3& blocks){ + int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); + + //create one thread for each detector pixel + blocks = dim3(SQRT_BLOCK, SQRT_BLOCK); + grids = dim3((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); + } + + //find the maximum value of component n + T find_max(unsigned int n){ + cublasStatus_t stat; + cublasHandle_t handle; + + //create a CUBLAS handle + stat = cublasCreate(&handle); + if(stat != CUBLAS_STATUS_SUCCESS){ + std::cout<<"CUBLAS Error: initialization failed"< process_filename(std::string name){ + std::stringstream ss(name); + std::string item; + std::vector elems; + while(std::getline(ss, item, '.')) //split the string at the '.' character (filename and extension) + { + elems.push_back(item); + } + + std::string prefix = elems[0]; //prefix contains the filename (with wildcard '?' characters) + std::string ext = elems[1]; //file extension (ex. .bmp, .png) + ext = std::string(".") + ext; //add a period back into the extension + + size_t i0 = prefix.find_first_of("?"); //find the positions of the first and last wildcard ('?'') + size_t i1 = prefix.find_last_of("?"); + + std::string postfix = prefix.substr(i1+1); + prefix = prefix.substr(0, i0); + + unsigned int digits = i1 - i0 + 1; //compute the number of wildcards + + std::vector flist; //create a vector of file names + //fill the list + for(unsigned int d=0; d>> (X[n], rhs, R[0], R[1]); + + return *this; + } + + //assignment of vector component + field & operator= (const vec rhs){ + + int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); + + //create one thread for each detector pixel + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); + dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); + + //assign the constant value to all positions and dimensions + for(unsigned int n=0; n>> (X[n], rhs.v[n], R[0], R[1]); + + return *this; + + } + + //multiply two fields (element-wise multiplication) + field operator* (const field & rhs){ + + int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); + + //create one thread for each detector pixel + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); + dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); + + //create a scalar field to store the result + field result(R[0], R[1]); + + for(int n=0; n>> (result.X[n], X[n], rhs.X[n], R[0], R[1]); + + return result; + } + + T* ptr(unsigned int n = 0){ + if(n < D) + return X[n]; + else return NULL; + } + + //return the vector component at position (u, v) + vec get(unsigned int u, unsigned int v){ + + vec result; + for(unsigned int d=0; d crop(unsigned int width, unsigned int height){ + int maxThreads = stim::maxThreadsPerBlock(); //compute the optimal block size + int SQRT_BLOCK = (int)std::sqrt((float)maxThreads); + + //create one thread for each detector pixel + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); + dim3 dimGrid((width + SQRT_BLOCK -1)/SQRT_BLOCK, (height + SQRT_BLOCK - 1)/SQRT_BLOCK); + + //create a scalar field to store the result + field result(width, height); + + for(int n=0; n>> (result.X[n], X[n], R[0], R[1], width, height); + + return result; + } + + //save an image representing component n + void toImage(std::string filename, unsigned int n = 0, + bool positive = false, stim::colormapType cmap = stim::cmBrewer){ + T max_val = find_max(n); //find the maximum value + + if(positive) //if the field is positive, use the range [0 max_val] + stim::gpu2image(X[n], filename, R[0], R[1], 0, max_val, cmap); + else + stim::gpu2image(X[n], filename, R[0], R[1], -max_val, max_val, cmap); + } + +}; + +} //end namespace rts +#endif diff --git a/math/rect.h b/math/rect.h index eba4448..7dc48b5 100644 --- a/math/rect.h +++ b/math/rect.h @@ -1,219 +1,219 @@ -#ifndef RTS_RECT_H -#define RTS_RECT_H - -//enable CUDA_CALLABLE macro -#include "../cuda/callable.h" -#include "../math/vector.h" -#include "../math/triangle.h" -#include "../math/quaternion.h" -#include -#include -#include - -namespace stim{ - -//template for a rectangle class in ND space -template -struct rect -{ - /* - ^ O - | - | - Y C - | - | - O---------X---------> - */ - -private: - - stim::vec C; - stim::vec X; - stim::vec Y; - - CUDA_CALLABLE void scale(T factor){ - X *= factor; - Y *= factor; - } - - CUDA_CALLABLE void normal(vec n){ //orient the rectangle along the specified normal - - n = n.norm(); //normalize, just in case - vec n_current = X.cross(Y).norm(); //compute the current normal - quaternion q; //create a quaternion - q.CreateRotation(n_current, n); //initialize a rotation from n_current to n - - //apply the quaternion to the vectors and position - X = q.toMatrix3() * X; - Y = q.toMatrix3() * Y; - } - - CUDA_CALLABLE void init(){ - C = vec(0, 0, 0); - X = vec(1, 0, 0); - Y = vec(0, 1, 0); - } - -public: - - CUDA_CALLABLE rect(){ - init(); - } - - CUDA_CALLABLE rect(T size, T z_pos = (T)0){ - init(); //use the default setup - scale(size); //scale the rectangle - C[2] = z_pos; - } - - CUDA_CALLABLE rect(T size, vec c, vec n = vec(0, 0, 1)){ - init(); //start with the default setting - C = c; - scale(size); //scale the rectangle - normal(n); //orient - - } - - /*CUDA_CALLABLE rect(vec a, vec b, vec c) - { - A = a; - Y = b - a; - X = c - a - Y; - - }*/ - - /******************************************************************* - Constructor - create a rect from a position, normal, and rotation - *******************************************************************/ - /*CUDA_CALLABLE rect(stim::vec c, stim::vec normal, T width, T height, T theta) - { - - //compute the X direction - start along world-space X - Y = stim::vec(0, 1, 0); - if(Y == normal) - Y = stim::vec(0, 0, 1); - - X = Y.cross(normal).norm(); - - std::cout< q; - q.CreateRotation(theta, normal); - X = q.toMatrix3() * X; - Y = normal.cross(X); - - //normalize everything - X = X.norm(); - Y = Y.norm(); - - //scale to match the rect width and height - X = X * width; - Y = Y * height; - - //set the corner of the plane - A = c - X * 0.5f - Y * 0.5f; - - std::cout< & rhs) - { - if(C == rhs.C && X == rhs.X && Y == rhs.Y) - return true; - else - return false; - } - - /******************************************* - Return the normal for the rect - *******************************************/ - CUDA_CALLABLE stim::vec n() - { - return (X.cross(Y)).norm(); - } - - CUDA_CALLABLE stim::vec p(T a, T b) - { - stim::vec result; - //given the two parameters a, b = [0 1], returns the position in world space - vec A = C - X * (T)0.5 - Y * (T)0.5; - result = A + X * a + Y * b; - - return result; - } - - CUDA_CALLABLE stim::vec operator()(T a, T b) - { - return p(a, b); - } - - std::string str() - { - std::stringstream ss; - vec A = C - X * (T)0.5 - Y * (T)0.5; - ss<"<<"C="<"<<"D="< operator*(T rhs) - { - //scales the plane by a scalar value - - //create the new rectangle - rect result = *this; - result.scale(rhs); - - return result; - - } - - CUDA_CALLABLE T dist(vec p) - { - //compute the distance between a point and this rect - - vec A = C - X * (T)0.5 - Y * (T)0.5; - - //first break the rect up into two triangles - triangle T0(A, A+X, A+Y); - triangle T1(A+X+Y, A+X, A+Y); - - - T d0 = T0.dist(p); - T d1 = T1.dist(p); - - if(d0 < d1) - return d0; - else - return d1; - } - - CUDA_CALLABLE T dist_max(vec p) - { - vec A = C - X * (T)0.5 - Y * (T)0.5; - T da = (A - p).len(); - T db = (A+X - p).len(); - T dc = (A+Y - p).len(); - T dd = (A+X+Y - p).len(); - - return std::max( da, std::max(db, std::max(dc, dd) ) ); - } -}; - -} //end namespace rts - -template -std::ostream& operator<<(std::ostream& os, stim::rect R) -{ - os< +#include +#include + +namespace stim{ + +//template for a rectangle class in ND space +template +struct rect +{ + /* + ^ O + | + | + Y C + | + | + O---------X---------> + */ + +private: + + stim::vec C; + stim::vec X; + stim::vec Y; + + CUDA_CALLABLE void scale(T factor){ + X *= factor; + Y *= factor; + } + + CUDA_CALLABLE void normal(vec n){ //orient the rectangle along the specified normal + + n = n.norm(); //normalize, just in case + vec n_current = X.cross(Y).norm(); //compute the current normal + quaternion q; //create a quaternion + q.CreateRotation(n_current, n); //initialize a rotation from n_current to n + + //apply the quaternion to the vectors and position + X = q.toMatrix3() * X; + Y = q.toMatrix3() * Y; + } + + CUDA_CALLABLE void init(){ + C = vec(0, 0, 0); + X = vec(1, 0, 0); + Y = vec(0, 1, 0); + } + +public: + + CUDA_CALLABLE rect(){ + init(); + } + + CUDA_CALLABLE rect(T size, T z_pos = (T)0){ + init(); //use the default setup + scale(size); //scale the rectangle + C[2] = z_pos; + } + + CUDA_CALLABLE rect(T size, vec c, vec n = vec(0, 0, 1)){ + init(); //start with the default setting + C = c; + scale(size); //scale the rectangle + normal(n); //orient + + } + + /*CUDA_CALLABLE rect(vec a, vec b, vec c) + { + A = a; + Y = b - a; + X = c - a - Y; + + }*/ + + /******************************************************************* + Constructor - create a rect from a position, normal, and rotation + *******************************************************************/ + /*CUDA_CALLABLE rect(stim::vec c, stim::vec normal, T width, T height, T theta) + { + + //compute the X direction - start along world-space X + Y = stim::vec(0, 1, 0); + if(Y == normal) + Y = stim::vec(0, 0, 1); + + X = Y.cross(normal).norm(); + + std::cout< q; + q.CreateRotation(theta, normal); + X = q.toMatrix3() * X; + Y = normal.cross(X); + + //normalize everything + X = X.norm(); + Y = Y.norm(); + + //scale to match the rect width and height + X = X * width; + Y = Y * height; + + //set the corner of the plane + A = c - X * 0.5f - Y * 0.5f; + + std::cout< & rhs) + { + if(C == rhs.C && X == rhs.X && Y == rhs.Y) + return true; + else + return false; + } + + /******************************************* + Return the normal for the rect + *******************************************/ + CUDA_CALLABLE stim::vec n() + { + return (X.cross(Y)).norm(); + } + + CUDA_CALLABLE stim::vec p(T a, T b) + { + stim::vec result; + //given the two parameters a, b = [0 1], returns the position in world space + vec A = C - X * (T)0.5 - Y * (T)0.5; + result = A + X * a + Y * b; + + return result; + } + + CUDA_CALLABLE stim::vec operator()(T a, T b) + { + return p(a, b); + } + + std::string str() + { + std::stringstream ss; + vec A = C - X * (T)0.5 - Y * (T)0.5; + ss<"<<"C="<"<<"D="< operator*(T rhs) + { + //scales the plane by a scalar value + + //create the new rectangle + rect result = *this; + result.scale(rhs); + + return result; + + } + + CUDA_CALLABLE T dist(vec p) + { + //compute the distance between a point and this rect + + vec A = C - X * (T)0.5 - Y * (T)0.5; + + //first break the rect up into two triangles + triangle T0(A, A+X, A+Y); + triangle T1(A+X+Y, A+X, A+Y); + + + T d0 = T0.dist(p); + T d1 = T1.dist(p); + + if(d0 < d1) + return d0; + else + return d1; + } + + CUDA_CALLABLE T dist_max(vec p) + { + vec A = C - X * (T)0.5 - Y * (T)0.5; + T da = (A - p).len(); + T db = (A+X - p).len(); + T dc = (A+Y - p).len(); + T dd = (A+X+Y - p).len(); + + return std::max( da, std::max(db, std::max(dc, dd) ) ); + } +}; + +} //end namespace rts + +template +std::ostream& operator<<(std::ostream& os, stim::rect R) +{ + os< - -namespace stim{ - -template -struct triangle -{ - /* - A------>B - | / - | / - | / - | / - | / - | / - C - */ - private: - - vec A; - vec B; - vec C; - - CUDA_CALLABLE vec _p(T s, T t) - { - //This function returns the point specified by p = A + s(B-A) + t(C-A) - vec E0 = B-A; - vec E1 = C-A; - - return A + s*E0 + t*E1; - } - - - public: - - - - CUDA_CALLABLE triangle() - { - - } - - CUDA_CALLABLE triangle(vec a, vec b, vec c) - { - A = a; - B = b; - C = c; - } - - CUDA_CALLABLE stim::vec operator()(T s, T t) - { - return _p(s, t); - } - - CUDA_CALLABLE vec nearest(vec p) - { - //comptue the distance between a point and this triangle - // This code is adapted from: http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf - - vec E0 = B-A; - vec E1 = C-A; - vec D = A - p; - - T a = E0.dot(E0); - T b = E0.dot(E1); - T c = E1.dot(E1); - T d = E0.dot(D); - T e = E1.dot(D); - //T f = D.dot(D); - - T det = a*c - b*b; - T s = b*e - c*d; - T t = b*d - a*e; - - /*std::cout<<"E0: "<= 0 ? 0 : ( -e >= c ? 1 : -e/c ) ); - //done - } - } - else if(t < 0) - { - //region 5 - //std::cout<<"Region 5"<= 0 ? 0 : ( -d >= a ? 1 : -d/a ) ); - t = 0; - //done - } - else - { - //region 0 - //std::cout<<"Region 0"<= denom ? 1 : numer/denom ); - } - t = 1 - s; - //done - } - } - - //std::cout<<"s: "< p) - { - vec n = nearest(p); - - return (p - n).len(); - } -}; - -} - -#endif +#ifndef RTS_TRIANGLE_H +#define RTS_TRIANGLE_H + +//enable CUDA_CALLABLE macro +#include "../cuda/callable.h" +#include "../math/vector.h" +#include + +namespace stim{ + +template +struct triangle +{ + /* + A------>B + | / + | / + | / + | / + | / + | / + C + */ + private: + + vec A; + vec B; + vec C; + + CUDA_CALLABLE vec _p(T s, T t) + { + //This function returns the point specified by p = A + s(B-A) + t(C-A) + vec E0 = B-A; + vec E1 = C-A; + + return A + s*E0 + t*E1; + } + + + public: + + + + CUDA_CALLABLE triangle() + { + + } + + CUDA_CALLABLE triangle(vec a, vec b, vec c) + { + A = a; + B = b; + C = c; + } + + CUDA_CALLABLE stim::vec operator()(T s, T t) + { + return _p(s, t); + } + + CUDA_CALLABLE vec nearest(vec p) + { + //comptue the distance between a point and this triangle + // This code is adapted from: http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf + + vec E0 = B-A; + vec E1 = C-A; + vec D = A - p; + + T a = E0.dot(E0); + T b = E0.dot(E1); + T c = E1.dot(E1); + T d = E0.dot(D); + T e = E1.dot(D); + //T f = D.dot(D); + + T det = a*c - b*b; + T s = b*e - c*d; + T t = b*d - a*e; + + /*std::cout<<"E0: "<= 0 ? 0 : ( -e >= c ? 1 : -e/c ) ); + //done + } + } + else if(t < 0) + { + //region 5 + //std::cout<<"Region 5"<= 0 ? 0 : ( -d >= a ? 1 : -d/a ) ); + t = 0; + //done + } + else + { + //region 0 + //std::cout<<"Region 0"<= denom ? 1 : numer/denom ); + } + t = 1 - s; + //done + } + } + + //std::cout<<"s: "< p) + { + vec n = nearest(p); + + return (p - n).len(); + } +}; + +} + +#endif diff --git a/optics/material.h b/optics/material.h index f22fa0d..fb97cd6 100644 --- a/optics/material.h +++ b/optics/material.h @@ -1,153 +1,153 @@ -#ifndef RTS_MATERIAL_H -#define RTS_MATERIAL_H - -#include -#include -#include -#include -#include -#include -#include -#include "../math/complex.h" -#include "../math/constants.h" -#include "../math/function.h" - -namespace stim{ - -//Material class - default representation for the material property is the refractive index (RI) -template -class material : public function< T, complex >{ - -public: - enum wave_property{microns, inverse_cm}; - enum material_property{ri, absorbance}; - -private: - - using function< T, complex >::X; - using function< T, complex >::Y; - using function< T, complex >::insert; - using function< T, complex >::bounding; - - std::string name; //name for the material (defaults to file name) - - void process_header(std::string str, wave_property& wp, material_property& mp){ - - std::stringstream ss(str); //create a stream from the data string - std::string line; - std::getline(ss, line); //get the first line as a string - while(line[0] == '#'){ //continue looping while the line is a comment - - std::stringstream lstream(line); //create a stream from the line - lstream.ignore(); //ignore the first character ('#') - - std::string prop; //get the property name - lstream>>prop; - - if(prop == "X"){ - std::string wp_name; - lstream>>wp_name; - if(wp_name == "microns") wp = microns; - else if(wp_name == "inverse_cm") wp = inverse_cm; - } - else if(prop == "Y"){ - std::string mp_name; - lstream>>mp_name; - if(mp_name == "ri") mp = ri; - else if(mp_name == "absorbance") mp = absorbance; - } - - std::getline(ss, line); //get the next line - } - - function< T, stim::complex >::process_string(str); - } - - void from_inverse_cm(){ - //convert inverse centimeters to wavelength (in microns) - for(unsigned int i=0; i(1, 0); - } - - -public: - - material(std::string filename, wave_property wp, material_property mp){ - name = filename; - load(filename, wp, mp); - } - - material(std::string filename){ - name = filename; - load(filename); - } - - material(){ - init(); - } - - complex getN(T lambda){ - return function< T, complex >::linear(lambda); - } - - void load(std::string filename, wave_property wp, material_property mp){ - - //load the file as a function - function< T, complex >::load(filename); - } - - void load(std::string filename){ - - wave_property wp = inverse_cm; - material_property mp = ri; - //turn the file into a string - std::ifstream t(filename.c_str()); //open the file as a stream - - if(!t){ - std::cout<<"ERROR: Couldn't open the material file '"<(t)), - std::istreambuf_iterator()); - - //process the header information - process_header(str, wp, mp); - - //convert units - if(wp == inverse_cm) - from_inverse_cm(); - //set the bounding values - bounding[0] = Y[0]; - bounding[1] = Y.back(); - } - std::string str(){ - std::stringstream ss; - ss< >::str(); - return ss.str(); - } - std::string get_name(){ - return name; - } - - void set_name(std::string str){ - name = str; - } - -}; - -} - - - - -#endif +#ifndef RTS_MATERIAL_H +#define RTS_MATERIAL_H + +#include +#include +#include +#include +#include +#include +#include +#include "../math/complex.h" +#include "../math/constants.h" +#include "../math/function.h" + +namespace stim{ + +//Material class - default representation for the material property is the refractive index (RI) +template +class material : public function< T, complex >{ + +public: + enum wave_property{microns, inverse_cm}; + enum material_property{ri, absorbance}; + +private: + + using function< T, complex >::X; + using function< T, complex >::Y; + using function< T, complex >::insert; + using function< T, complex >::bounding; + + std::string name; //name for the material (defaults to file name) + + void process_header(std::string str, wave_property& wp, material_property& mp){ + + std::stringstream ss(str); //create a stream from the data string + std::string line; + std::getline(ss, line); //get the first line as a string + while(line[0] == '#'){ //continue looping while the line is a comment + + std::stringstream lstream(line); //create a stream from the line + lstream.ignore(); //ignore the first character ('#') + + std::string prop; //get the property name + lstream>>prop; + + if(prop == "X"){ + std::string wp_name; + lstream>>wp_name; + if(wp_name == "microns") wp = microns; + else if(wp_name == "inverse_cm") wp = inverse_cm; + } + else if(prop == "Y"){ + std::string mp_name; + lstream>>mp_name; + if(mp_name == "ri") mp = ri; + else if(mp_name == "absorbance") mp = absorbance; + } + + std::getline(ss, line); //get the next line + } + + function< T, stim::complex >::process_string(str); + } + + void from_inverse_cm(){ + //convert inverse centimeters to wavelength (in microns) + for(unsigned int i=0; i(1, 0); + } + + +public: + + material(std::string filename, wave_property wp, material_property mp){ + name = filename; + load(filename, wp, mp); + } + + material(std::string filename){ + name = filename; + load(filename); + } + + material(){ + init(); + } + + complex getN(T lambda){ + return function< T, complex >::linear(lambda); + } + + void load(std::string filename, wave_property wp, material_property mp){ + + //load the file as a function + function< T, complex >::load(filename); + } + + void load(std::string filename){ + + wave_property wp = inverse_cm; + material_property mp = ri; + //turn the file into a string + std::ifstream t(filename.c_str()); //open the file as a stream + + if(!t){ + std::cout<<"ERROR: Couldn't open the material file '"<(t)), + std::istreambuf_iterator()); + + //process the header information + process_header(str, wp, mp); + + //convert units + if(wp == inverse_cm) + from_inverse_cm(); + //set the bounding values + bounding[0] = Y[0]; + bounding[1] = Y.back(); + } + std::string str(){ + std::stringstream ss; + ss< >::str(); + return ss.str(); + } + std::string get_name(){ + return name; + } + + void set_name(std::string str){ + name = str; + } + +}; + +} + + + + +#endif diff --git a/optics/mirst-1d.cuh b/optics/mirst-1d.cuh index e7fde75..9279087 100644 --- a/optics/mirst-1d.cuh +++ b/optics/mirst-1d.cuh @@ -1,447 +1,447 @@ -#include "../optics/material.h" -#include "../math/complexfield.cuh" -#include "../math/constants.h" -//#include "../envi/bil.h" - -#include "cufft.h" - -#include -#include - -namespace stim{ - -//this function writes a sinc function to "dest" such that an iFFT produces a slab -template -__global__ void gpu_mirst1d_layer_fft(complex* dest, complex* ri, - T* src, T* zf, - T w, unsigned int zR, unsigned int nuR){ - //dest = complex field representing the sample - //ri = refractive indices for each wavelength - //src = intensity of the light source for each wavelength - //zf = z position of the slab interface for each wavelength (accounting for optical path length) - //w = width of the slab (in pixels) - //zR = number of z-axis samples - //nuR = number of wavelengths - - //get the current coordinate in the plane slice - int ifz = blockIdx.x * blockDim.x + threadIdx.x; - int inu = blockIdx.y * blockDim.y + threadIdx.y; - - //make sure that the thread indices are in-bounds - if(inu >= nuR || ifz >= zR) return; - - int i = inu * zR + ifz; - - T fz; - if(ifz < zR/2) - fz = ifz / (T)zR; - else - fz = -(zR - ifz) / (T)zR; - - //if the slab starts outside of the simulation domain, just return - if(zf[inu] >= zR) return; - - //fill the array along z with a sinc function representing the Fourier transform of the layer - - T opl = w * ri[inu].real(); //optical path length - - //handle the case where the slab goes outside the simulation domain - if(zf[inu] + opl >= zR) - opl = zR - zf[inu]; - - if(opl == 0) return; - - //T l = w * ri[inu].real(); - //complex e(0.0, -2 * PI * fz * (zf[inu] + zR/2 - l/2.0)); - complex e(0, -2 * stimPI * fz * (zf[inu] + opl/2)); - - complex eta = ri[inu] * ri[inu] - 1; - - //dest[i] = fz;//exp(e) * m[inu] * src[inu] * sin(PI * fz * l) / (PI * fz); - if(ifz == 0) - dest[i] += opl * exp(e) * eta * src[inu]; - else - dest[i] += opl * exp(e) * eta * src[inu] * sin(stimPI * fz * opl) / (stimPI * fz * opl); -} - -template -__global__ void gpu_mirst1d_increment_z(T* zf, complex* ri, T w, unsigned int S){ - //zf = current z depth (optical path length) in pixels - //ri = refractive index of the material - //w = actual width of the layer (in pixels) - - - //compute the index for this thread - int i = blockIdx.x * blockDim.x + threadIdx.x; - if(i >= S) return; - - if(ri == NULL) - zf[i] += w; - else - zf[i] += ri[i].real() * w; -} - -//apply the 1D MIRST filter to an existing sample (overwriting the sample) -template -__global__ void gpu_mirst1d_apply_filter(complex* sampleFFT, T* lambda, - T dFz, - T inNA, T outNA, - unsigned int lambdaR, unsigned int zR, - T sigma = 0){ - //sampleFFT = the sample in the Fourier domain (will be overwritten) - //lambda = list of wavelengths - //dFz = delta along the Fz axis in the frequency domain - //inNA = NA of the internal obscuration - //outNA = NA of the objective - //zR = number of pixels along the Fz axis (same as the z-axis) - //lambdaR = number of wavelengths - //sigma = width of the Gaussian source - int ifz = blockIdx.x * blockDim.x + threadIdx.x; - int inu = blockIdx.y * blockDim.y + threadIdx.y; - - if(inu >= lambdaR || ifz >= zR) return; - - //calculate the index into the sample FT - int i = inu * zR + ifz; - - //compute the frequency (and set all negative spatial frequencies to zero) - T fz; - if(ifz < zR / 2) - fz = ifz * dFz; - //if the spatial frequency is negative, set it to zero and exit - else{ - sampleFFT[i] = 0; - return; - } - - //compute the frequency in inverse microns - T nu = 1/lambda[inu]; - - //determine the radius of the integration circle - T nu_sq = nu * nu; - T fz_sq = (fz * fz) / 4; - - //cut off frequencies above the diffraction limit - T r; - if(fz_sq < nu_sq) - r = sqrt(nu_sq - fz_sq); - else - r = 0; - - //account for the optics - T Q = 0; - if(r > nu * inNA && r < nu * outNA) - Q = 1; - - //account for the source - //T sigma = 30.0; - T s = exp( - (r*r * sigma*sigma) / 2 ); - //T s=1; - - //compute the final filter - T mirst = 0; - if(fz != 0) - mirst = 2 * stimPI * r * s * Q * (1/fz); - - sampleFFT[i] *= mirst; - -} - -/*This object performs a 1-dimensional (layered) MIRST simulation -*/ -template -class mirst1d{ - -private: - unsigned int Z; //z-axis resolution - unsigned int pad; //pixel padding on either side of the sample - - std::vector< material > matlist; //list of materials - std::vector< T > layers; //list of layer thicknesses - - std::vector< T > lambdas; //list of wavelengths that are being simulated - unsigned int S; //number of wavelengths (size of "lambdas") - - T NA[2]; //numerical aperature (central obscuration and outer diameter) - - function source_profile; //profile (spectrum) of the source (expressed in inverse centimeters) - - complexfield scratch; //scratch GPU memory used to build samples, transforms, etc. - - void fft(int direction = CUFFT_FORWARD){ - - unsigned padZ = Z + pad; - - //create cuFFT handles - cufftHandle plan; - cufftResult result; - - if(sizeof(T) == 4) - result = cufftPlan1d(&plan, padZ, CUFFT_C2C, lambdas.size()); //single precision - else - result = cufftPlan1d(&plan, padZ, CUFFT_Z2Z, lambdas.size()); //double precision - - //check for Plan 1D errors - if(result != CUFFT_SUCCESS){ - std::cout<<"Error creating CUFFT plan for computing the FFT:"<(Z + pad , lambdas.size()); - scratch = 0; - } - - //get the list of scattering efficiency (eta) values for a specified layer - std::vector< complex > layer_etas(unsigned int l){ - - std::vector< complex > etas; - - //fill the list of etas - for(unsigned int i=0; i* gpuRi; - HANDLE_ERROR(cudaMalloc( (void**)&gpuRi, sizeof(complex) * S)); - - //allocate memory for the source profile - T* gpuSrc; - HANDLE_ERROR(cudaMalloc( (void**)&gpuSrc, sizeof(T) * S)); - - complex ri; - T source; - //store the refractive index and source profile in a CPU array - for(int inu=0; inu), cudaMemcpyHostToDevice )); - - //save the source profile to the GPU - source = source_profile(10000 / lambdas[inu]); - HANDLE_ERROR(cudaMemcpy( gpuSrc + inu, &source, sizeof(T), cudaMemcpyHostToDevice )); - - } - - //create one thread for each pixel of the field slice - dim3 gridDim, blockDim; - cuda_params(gridDim, blockDim); - stim::gpu_mirst1d_layer_fft<<>>(scratch.ptr(), gpuRi, gpuSrc, zf, wpx, paddedZ, S); - - int linBlock = stim::maxThreadsPerBlock(); //compute the optimal block size - int linGrid = S / linBlock + 1; - stim::gpu_mirst1d_increment_z <<>>(zf, gpuRi, wpx, S); - - //free memory - HANDLE_ERROR(cudaFree(gpuRi)); - HANDLE_ERROR(cudaFree(gpuSrc)); - } - - void build_sample(){ - init_scratch(); //initialize the GPU scratch space - //build_layer(1); - - T* zf; - HANDLE_ERROR(cudaMalloc(&zf, sizeof(T) * S)); - HANDLE_ERROR(cudaMemset(zf, 0, sizeof(T) * S)); - - //render each layer of the sample - for(unsigned int l=0; l>>(scratch.ptr(), gpuLambdas, - dFz, - NA[0], NA[1], - S, Zpad); - } - - //crop the image to the sample thickness - keep in mind that sample thickness != optical path length - void crop(){ - - scratch = scratch.crop(Z, S); - } - - //save the scratch field as a binary file - void to_binary(std::string filename){ - - } - - -public: - - //constructor - mirst1d(unsigned int rZ = 100, - unsigned int padding = 0){ - Z = rZ; - pad = padding; - NA[0] = 0; - NA[1] = 0.8; - S = 0; - source_profile = 1; - } - - //add a layer, thickness = microns - void add_layer(material mat, T thickness){ - matlist.push_back(mat); - layers.push_back(thickness); - } - - void add_layer(std::string filename, T thickness){ - add_layer(material(filename), thickness); - } - - //adds a profile spectrum for the light source - void set_source(std::string filename){ - source_profile.load(filename); - } - - //adds a block of wavenumbers (cm^-1) to the simulation parameters - void add_wavenumbers(unsigned int start, unsigned int stop, unsigned int step){ - unsigned int nu = start; - while(nu <= stop){ - lambdas.push_back((T)10000 / nu); - nu += step; - } - S = lambdas.size(); //increment the number of wavelengths (shorthand for later) - } - - T thickness(){ - T t = 0; - for(unsigned int l=0; l get_source(){ - return source_profile; - } - - void save_sample(std::string filename){ - //create a sample and save the magnitude as an image - build_sample(); - fft(CUFFT_INVERSE); - scratch.toImage(filename, stim::complexfield::magnitude); - } - - void save_mirst(std::string filename, bool binary = true){ - //apply the MIRST filter to a sample and save the image - - //build the sample in the Fourier domain - build_sample(); - - //apply the MIRST filter - apply_filter(); - - //apply an inverse FFT to bring the results back into the spatial domain - fft(CUFFT_INVERSE); - - crop(); - - //save the image - if(binary) - to_binary(filename); - else - scratch.toImage(filename, stim::complexfield::magnitude); - } - - - - - std::string str(){ - - stringstream ss; - ss<<"1D MIRST Simulation========================="< +#include + +namespace stim{ + +//this function writes a sinc function to "dest" such that an iFFT produces a slab +template +__global__ void gpu_mirst1d_layer_fft(complex* dest, complex* ri, + T* src, T* zf, + T w, unsigned int zR, unsigned int nuR){ + //dest = complex field representing the sample + //ri = refractive indices for each wavelength + //src = intensity of the light source for each wavelength + //zf = z position of the slab interface for each wavelength (accounting for optical path length) + //w = width of the slab (in pixels) + //zR = number of z-axis samples + //nuR = number of wavelengths + + //get the current coordinate in the plane slice + int ifz = blockIdx.x * blockDim.x + threadIdx.x; + int inu = blockIdx.y * blockDim.y + threadIdx.y; + + //make sure that the thread indices are in-bounds + if(inu >= nuR || ifz >= zR) return; + + int i = inu * zR + ifz; + + T fz; + if(ifz < zR/2) + fz = ifz / (T)zR; + else + fz = -(zR - ifz) / (T)zR; + + //if the slab starts outside of the simulation domain, just return + if(zf[inu] >= zR) return; + + //fill the array along z with a sinc function representing the Fourier transform of the layer + + T opl = w * ri[inu].real(); //optical path length + + //handle the case where the slab goes outside the simulation domain + if(zf[inu] + opl >= zR) + opl = zR - zf[inu]; + + if(opl == 0) return; + + //T l = w * ri[inu].real(); + //complex e(0.0, -2 * PI * fz * (zf[inu] + zR/2 - l/2.0)); + complex e(0, -2 * stimPI * fz * (zf[inu] + opl/2)); + + complex eta = ri[inu] * ri[inu] - 1; + + //dest[i] = fz;//exp(e) * m[inu] * src[inu] * sin(PI * fz * l) / (PI * fz); + if(ifz == 0) + dest[i] += opl * exp(e) * eta * src[inu]; + else + dest[i] += opl * exp(e) * eta * src[inu] * sin(stimPI * fz * opl) / (stimPI * fz * opl); +} + +template +__global__ void gpu_mirst1d_increment_z(T* zf, complex* ri, T w, unsigned int S){ + //zf = current z depth (optical path length) in pixels + //ri = refractive index of the material + //w = actual width of the layer (in pixels) + + + //compute the index for this thread + int i = blockIdx.x * blockDim.x + threadIdx.x; + if(i >= S) return; + + if(ri == NULL) + zf[i] += w; + else + zf[i] += ri[i].real() * w; +} + +//apply the 1D MIRST filter to an existing sample (overwriting the sample) +template +__global__ void gpu_mirst1d_apply_filter(complex* sampleFFT, T* lambda, + T dFz, + T inNA, T outNA, + unsigned int lambdaR, unsigned int zR, + T sigma = 0){ + //sampleFFT = the sample in the Fourier domain (will be overwritten) + //lambda = list of wavelengths + //dFz = delta along the Fz axis in the frequency domain + //inNA = NA of the internal obscuration + //outNA = NA of the objective + //zR = number of pixels along the Fz axis (same as the z-axis) + //lambdaR = number of wavelengths + //sigma = width of the Gaussian source + int ifz = blockIdx.x * blockDim.x + threadIdx.x; + int inu = blockIdx.y * blockDim.y + threadIdx.y; + + if(inu >= lambdaR || ifz >= zR) return; + + //calculate the index into the sample FT + int i = inu * zR + ifz; + + //compute the frequency (and set all negative spatial frequencies to zero) + T fz; + if(ifz < zR / 2) + fz = ifz * dFz; + //if the spatial frequency is negative, set it to zero and exit + else{ + sampleFFT[i] = 0; + return; + } + + //compute the frequency in inverse microns + T nu = 1/lambda[inu]; + + //determine the radius of the integration circle + T nu_sq = nu * nu; + T fz_sq = (fz * fz) / 4; + + //cut off frequencies above the diffraction limit + T r; + if(fz_sq < nu_sq) + r = sqrt(nu_sq - fz_sq); + else + r = 0; + + //account for the optics + T Q = 0; + if(r > nu * inNA && r < nu * outNA) + Q = 1; + + //account for the source + //T sigma = 30.0; + T s = exp( - (r*r * sigma*sigma) / 2 ); + //T s=1; + + //compute the final filter + T mirst = 0; + if(fz != 0) + mirst = 2 * stimPI * r * s * Q * (1/fz); + + sampleFFT[i] *= mirst; + +} + +/*This object performs a 1-dimensional (layered) MIRST simulation +*/ +template +class mirst1d{ + +private: + unsigned int Z; //z-axis resolution + unsigned int pad; //pixel padding on either side of the sample + + std::vector< material > matlist; //list of materials + std::vector< T > layers; //list of layer thicknesses + + std::vector< T > lambdas; //list of wavelengths that are being simulated + unsigned int S; //number of wavelengths (size of "lambdas") + + T NA[2]; //numerical aperature (central obscuration and outer diameter) + + function source_profile; //profile (spectrum) of the source (expressed in inverse centimeters) + + complexfield scratch; //scratch GPU memory used to build samples, transforms, etc. + + void fft(int direction = CUFFT_FORWARD){ + + unsigned padZ = Z + pad; + + //create cuFFT handles + cufftHandle plan; + cufftResult result; + + if(sizeof(T) == 4) + result = cufftPlan1d(&plan, padZ, CUFFT_C2C, lambdas.size()); //single precision + else + result = cufftPlan1d(&plan, padZ, CUFFT_Z2Z, lambdas.size()); //double precision + + //check for Plan 1D errors + if(result != CUFFT_SUCCESS){ + std::cout<<"Error creating CUFFT plan for computing the FFT:"<(Z + pad , lambdas.size()); + scratch = 0; + } + + //get the list of scattering efficiency (eta) values for a specified layer + std::vector< complex > layer_etas(unsigned int l){ + + std::vector< complex > etas; + + //fill the list of etas + for(unsigned int i=0; i* gpuRi; + HANDLE_ERROR(cudaMalloc( (void**)&gpuRi, sizeof(complex) * S)); + + //allocate memory for the source profile + T* gpuSrc; + HANDLE_ERROR(cudaMalloc( (void**)&gpuSrc, sizeof(T) * S)); + + complex ri; + T source; + //store the refractive index and source profile in a CPU array + for(int inu=0; inu), cudaMemcpyHostToDevice )); + + //save the source profile to the GPU + source = source_profile(10000 / lambdas[inu]); + HANDLE_ERROR(cudaMemcpy( gpuSrc + inu, &source, sizeof(T), cudaMemcpyHostToDevice )); + + } + + //create one thread for each pixel of the field slice + dim3 gridDim, blockDim; + cuda_params(gridDim, blockDim); + stim::gpu_mirst1d_layer_fft<<>>(scratch.ptr(), gpuRi, gpuSrc, zf, wpx, paddedZ, S); + + int linBlock = stim::maxThreadsPerBlock(); //compute the optimal block size + int linGrid = S / linBlock + 1; + stim::gpu_mirst1d_increment_z <<>>(zf, gpuRi, wpx, S); + + //free memory + HANDLE_ERROR(cudaFree(gpuRi)); + HANDLE_ERROR(cudaFree(gpuSrc)); + } + + void build_sample(){ + init_scratch(); //initialize the GPU scratch space + //build_layer(1); + + T* zf; + HANDLE_ERROR(cudaMalloc(&zf, sizeof(T) * S)); + HANDLE_ERROR(cudaMemset(zf, 0, sizeof(T) * S)); + + //render each layer of the sample + for(unsigned int l=0; l>>(scratch.ptr(), gpuLambdas, + dFz, + NA[0], NA[1], + S, Zpad); + } + + //crop the image to the sample thickness - keep in mind that sample thickness != optical path length + void crop(){ + + scratch = scratch.crop(Z, S); + } + + //save the scratch field as a binary file + void to_binary(std::string filename){ + + } + + +public: + + //constructor + mirst1d(unsigned int rZ = 100, + unsigned int padding = 0){ + Z = rZ; + pad = padding; + NA[0] = 0; + NA[1] = 0.8; + S = 0; + source_profile = 1; + } + + //add a layer, thickness = microns + void add_layer(material mat, T thickness){ + matlist.push_back(mat); + layers.push_back(thickness); + } + + void add_layer(std::string filename, T thickness){ + add_layer(material(filename), thickness); + } + + //adds a profile spectrum for the light source + void set_source(std::string filename){ + source_profile.load(filename); + } + + //adds a block of wavenumbers (cm^-1) to the simulation parameters + void add_wavenumbers(unsigned int start, unsigned int stop, unsigned int step){ + unsigned int nu = start; + while(nu <= stop){ + lambdas.push_back((T)10000 / nu); + nu += step; + } + S = lambdas.size(); //increment the number of wavelengths (shorthand for later) + } + + T thickness(){ + T t = 0; + for(unsigned int l=0; l get_source(){ + return source_profile; + } + + void save_sample(std::string filename){ + //create a sample and save the magnitude as an image + build_sample(); + fft(CUFFT_INVERSE); + scratch.toImage(filename, stim::complexfield::magnitude); + } + + void save_mirst(std::string filename, bool binary = true){ + //apply the MIRST filter to a sample and save the image + + //build the sample in the Fourier domain + build_sample(); + + //apply the MIRST filter + apply_filter(); + + //apply an inverse FFT to bring the results back into the spatial domain + fft(CUFFT_INVERSE); + + crop(); + + //save the image + if(binary) + to_binary(filename); + else + scratch.toImage(filename, stim::complexfield::magnitude); + } + + + + + std::string str(){ + + stringstream ss; + ss<<"1D MIRST Simulation========================="< args; + //vector of options + std::vector opts; + std::vector args; - //column width of the longest argument + //column width of the longest option int col_width; //list of sections @@ -261,28 +262,28 @@ namespace stim{ void set_ansi(bool b) { ansi = b; - for(unsigned int i=0; i(col_width, arg.col_width()); + col_width = std::max(col_width, opt.col_width()); } void section(std::string _name) { argsection s; s.name = _name; - s.index = args.size(); + s.index = opts.size(); sections.push_back(s); } - //output the arguments (generally in response to --help) + //output the options (generally in response to --help) std::string str() { std::stringstream ss; @@ -292,8 +293,8 @@ namespace stim{ if(sections.size() > 0) si = 0; - //for each argument - for(unsigned int a=0; a= args.size()) + if(i >= opts.size()) return -1; return (int)i; @@ -327,52 +328,57 @@ namespace stim{ if(i != -1) { - args[i].set(_value); + opts[i].set(_value); //adjust the column width if necessary - col_width = (std::max)(col_width, args[i].col_width()); + col_width = (std::max)(col_width, opts[i].col_width()); } else - std::cout<<"ERROR - Argument not recognized: "<<_name<