Commit 8d8c3ebd609c8bc1175c0ce9ff1c92b558529b51

Authored by David Mayerich
0 parents

initial commit transferring bimsym to bimsim

Showing 4 changed files with 520 additions and 0 deletions   Show diff stats
CMakeLists.txt 0 → 100644
  1 +++ a/CMakeLists.txt
  1 +#Specify the version being used aswell as the language
  2 +cmake_minimum_required(VERSION 2.8)
  3 +#Name your project here
  4 +project(bimsym)
  5 +
  6 +
  7 +
  8 +#set the module directory
  9 +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}")
  10 +
  11 +#set up CUDA
  12 +find_package(CUDA)
  13 +
  14 +#find OpenCV
  15 +find_package(OpenCV REQUIRED )
  16 +
  17 +find_package(Boost REQUIRED)
  18 +
  19 +#find the STIM library
  20 +find_package(STIM)
  21 +
  22 +#build the executable in the binary directory on MS Visual Studio
  23 +if ( MSVC )
  24 + SET( CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG "${OUTPUT_DIRECTORY}")
  25 + SET( CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE "${OUTPUT_DIRECTORY}")
  26 +endif ( MSVC )
  27 +
  28 +#set the include directories
  29 +include_directories(
  30 + ${CMAKE_CURRENT_BINARY_DIR}
  31 + ${CMAKE_CURRENT_SOURCE_DIR}
  32 + ${STIM_INCLUDE_DIRS}
  33 + ${Boost_INCLUDE_DIRS}
  34 +)
  35 +
  36 +#SET(CUDA_NVCC_FLAGS_RELEASE ${CUDA_NVCC_FLAGS};/maxregcount=60)
  37 +
  38 +#create an executable
  39 +cuda_add_executable(bimsym
  40 + bimsym.cu
  41 + sources.h
  42 +)
  43 +
  44 +#set the link libraries
  45 +target_link_libraries(bimsym
  46 + ${CUDA_cufft_LIBRARY}
  47 + ${CUDA_cublas_LIBRARY}
  48 + ${OpenCV_LIBS}
  49 + )
  50 +
  51 +
  52 +
FindSTIM.cmake 0 → 100644
  1 +++ a/FindSTIM.cmake
  1 +include(FindPackageHandleStandardArgs)
  2 +
  3 +set(STIM_INCLUDE_DIR $ENV{STIMLIB_PATH})
  4 +
  5 +find_package_handle_standard_args(STIM DEFAULT_MSG STIM_INCLUDE_DIR)
  6 +
  7 +if(STIM_FOUND)
  8 + set(STIM_INCLUDE_DIRS ${STIM_INCLUDE_DIR})
  9 +endif()
bimsym.cu 0 → 100644
  1 +++ a/bimsym.cu
  1 +#define CUDA_FOUND
  2 +#define USING_OPENCV
  3 +
  4 +#include <iostream>
  5 +
  6 +#include <stim/visualization/colormap.h>
  7 +
  8 +#include <stim/parser/arguments.h>
  9 +#include <stim/math/vec3.h>
  10 +#include <stim/optics/scalarbeam.h>
  11 +#include <stim/optics/scalarmie.h>
  12 +#include <stim/parser/filename.h>
  13 +#include <stim/ui/progressbar.h>
  14 +#include <stim/parser/table.h>
  15 +
  16 +#include <stim/cuda/cudatools.h>
  17 +
  18 +#include <sources.h> //default source types
  19 +
  20 +#include <chrono>
  21 +#include <thread>
  22 +
  23 +//field parameters
  24 +size_t g_R_nf; //near-field resolution along X and Y
  25 +double g_S_nf; //size of the near field simulation
  26 +double g_Z_nf; //z position of the field slice
  27 +std::vector<stim::vec3<float>> g_sources; //point source position list
  28 +
  29 +bool use_cuda = true;
  30 +int cuda_device;
  31 +int cuda_major = 3;
  32 +int cuda_minor = 0;
  33 +bool use_backprop = true;
  34 +
  35 +//incident beam parameters
  36 +double g_NA_cond[2]; //condenser numerical aperture, [internal, external]
  37 +double g_lambda; //incident wavelength
  38 +size_t g_Nl; //number of orders used to calculate the incident field
  39 +double g_A; //amplitude of the incident field
  40 +
  41 +//scattering parameters
  42 +stim::scalarcluster<float> spheres;
  43 +//std::vector<double> g_radius(1); //radius of the scattering sphere
  44 +//std::vector< stim::complex<double> > g_n(1); //refractive index of the scattering sphere
  45 +size_t g_MC; //number of monte-carlo samples used to calculate the scattered field
  46 +//std::vector< stim::vec3<float> > g_c(1); //center of the sphere
  47 +
  48 +//far field parameters
  49 +double g_NA_obj[2]; //objective numerical aperture
  50 +size_t g_padding; //padding applied to the near-field calculation in order to improve the quality of the band-pass
  51 +double g_backprop = 0; //number of units that the field has to be backpropagated by in order to approximate a near-field slice passing through a sphere
  52 +size_t g_R_ff; //number of samples along X and Y in the far field image
  53 +double g_S_ff; //size (in units) of the far field image
  54 +double g_Z_ff; //position of the near-field plane projected into the far field (position of the back-propagated plane)
  55 +
  56 +//debugging
  57 +bool g_debug_images = false;
  58 +bool g_verbose = false;
  59 +
  60 +stim::filename outfile("out.bmp"); //output file name
  61 +
  62 +stim::arglist args; //input arguments
  63 +
  64 +template<typename T>
  65 +__global__ void cuda_absorbance(T* out, T* I, T* I0, size_t N){
  66 +
  67 + size_t i = blockIdx.x * blockDim.x + threadIdx.x;
  68 +
  69 + if(i >= N) return;
  70 +
  71 + out[i] = -log10(I[i] / I0[i]);
  72 +}
  73 +
  74 +/// Crops a 2D image composed of elements of type T
  75 +/// @param dest is a device pointer to memory of size dx*dy that will store the cropped image
  76 +/// @param src is a device pointer to memory of size sx*sy that stores the original image
  77 +/// @param sx is the size of the source image along x
  78 +/// @param sy is the size of the source image along y
  79 +/// @param x is the x-coordinate of the start position of the cropped region within the source image
  80 +/// @param y is the y-coordinate of the start position of the cropped region within the source image
  81 +/// @param dx is the size of the destination image along x
  82 +/// @param dy is the size of the destination image along y
  83 +template<typename T>
  84 +void gpu_absorbance(T* A, T* I, T* I0, size_t N){
  85 + int threads = stim::maxThreadsPerBlock(); //get the maximum number of threads per block for the CUDA device
  86 +
  87 + dim3 blocks( N / threads + 1 ); //calculate the optimal number of blocks
  88 + cuda_absorbance<T> <<< blocks, threads >>>(A, I, I0, N);
  89 +}
  90 +
  91 +//function to display a progress bar
  92 +void progressbar_thread(double* e){
  93 +
  94 + unsigned int p = 0;
  95 + while(p != 100){
  96 + if(*e > p){
  97 + p = (unsigned int)(*e);
  98 + rtsProgressBar(p);
  99 + }
  100 + }
  101 + std::cout<<std::endl; //put a newline after the completed progress bar
  102 +}
  103 +
  104 +void advertise(){
  105 + //output advertisement
  106 + std::cout<<std::endl<<std::endl;
  107 + std::cout<<"========================================================================="<<std::endl;
  108 +
  109 + std::cout<<"========================================================================="<<std::endl<<std::endl;
  110 +
  111 + std::cout<<std::endl<<std::endl;
  112 +
  113 + std::cout<<args.str();
  114 +}
  115 +
  116 +void set_arguments(int argc, char* argv[]){
  117 + args.add("help", "prints this help text");
  118 +
  119 + args.section("Field Slice Parameters");
  120 + args.add("res", "field resolution (number of samples along X and Y)", "256", "nonzero positive integer");
  121 + args.add("size", "size of the calculated field", "10", "nonzero positive value");
  122 + args.add("z", "position of the field slice along the z axis", "0", "any real value");
  123 +
  124 + args.section("Incoherent Sources");
  125 + args.add("simage", "image of the source at the focal plane (z=0)", "", "any standard image format");
  126 + args.add("sgrid", "grid of n x n point sources", "", "n");
  127 +
  128 + args.section("Incident Field Parameters");
  129 + args.add("f", "position of a single focal point", "0 0", "x y position (z = 0)");
  130 + args.add("condenser", "condenser numerical aperture (NA) and center obscuration", "1", "'sin(a1)' or 'sin(a1) sin(a0)', where a0 = center obscuration");
  131 + args.add("lambda", "incident wavelength", "1", "any real positive value");
  132 + args.add("wavenumber", "incident wavelength specified in cm^(-1)", "", "any real positive value, sets all units to microns");
  133 + args.add("order", "order used to calculate the incident field", "20", "any nonzero integer");
  134 +
  135 + args.section("Mie Scattering Parameters");
  136 + args.add("c", "center of the scattering sphere", "0 0 0", "any real 3D coordinate");
  137 + args.add("radius", "radius of the scattering sphere", "0.5", "any real positive value");
  138 + args.add("n", "complex refractive index of the scattering sphere", "1.4 0.1", "[r] or [r i]");
  139 + args.add("samples", "number of Monte-Carlo samples used to calculate the scattered field", "1000", "any real positive integer");
  140 + args.add("spheres", "specify a file defining several spheres", "", "(1/line: radius, real RI, imag RI, x, y, z)");
  141 +
  142 + args.section("Far Field Parameters");
  143 + args.add("objective", "objective numerical aperture (NA) and center obscuration", "1", "'sin(b1)' or 'sin(b1) sin(b0)', where b0 = center obscuration");
  144 + args.add("padding", "padding for the near field calculation (improves far-field calculation)", "1", "any real positive integer");
  145 + args.add("nobackprop", "don't use back-propagation to deal with sphere intersections");
  146 +
  147 + args.section("Debugging");
  148 + args.add("cuda", "CUDA device ID of device to use (-1 disables CUDA)", "0", "integer value");
  149 + args.add("debug", "outputs a debug image at several important intermediate steps of the Mie scattering calculation");
  150 + args.add("verbose", "outputs additional information regarding the Mie calculations as they are performed");
  151 + args.parse(argc, argv);
  152 +
  153 + if(args["help"].is_set()){
  154 + advertise();
  155 + exit(1);
  156 + }
  157 +}
  158 +
  159 +void output_simulation(){
  160 + std::cout<<"far field slice-------"<<std::endl;
  161 + std::cout<<"R: "<<g_R_ff<<" x "<<g_R_ff<<" pixels"<<std::endl;
  162 + std::cout<<"S: "<<g_S_ff<<" x "<<g_S_ff<<" units"<<std::endl;
  163 + std::cout<<"Z: "<<g_Z_ff<<" units"<<std::endl;
  164 +
  165 + std::cout<<"objective NA: "<<g_NA_obj[1];
  166 + if(g_NA_obj[0] != 0) std::cout<<" (internal = "<<g_NA_obj[0]<<")";
  167 + std::cout<<std::endl;
  168 + std::cout<<"padding: "<<g_padding<<std::endl;
  169 + std::cout<<"backprop: "<<g_backprop<<" units"<<std::endl;
  170 + std::cout<<std::endl;
  171 +
  172 + std::cout<<"near field slice------"<<std::endl;
  173 + std::cout<<"R: "<<g_R_nf<<" x "<<g_R_nf<<" pixels"<<std::endl;
  174 + std::cout<<"S: "<<g_S_nf<<" x "<<g_S_nf<<" units"<<std::endl;
  175 + std::cout<<"Z: "<<g_Z_nf<<" units"<<std::endl;
  176 + std::cout<<std::endl;
  177 +
  178 + std::cout<<"incident field--------"<<std::endl;
  179 + std::cout<<"condenser NA: "<<g_NA_cond[1];
  180 + if(g_NA_cond[0] != 0) std::cout<<" (internal = "<<g_NA_cond[0]<<")";
  181 + std::cout<<std::endl;
  182 + std::cout<<"wavelength: "<<g_lambda<<" units"<<std::endl;
  183 + std::cout<<"order: "<<g_Nl<<std::endl;
  184 + std::cout << "# sources: " << g_sources.size()<<" ";
  185 + if (g_sources.size() == 1)
  186 + std::cout << g_sources[0].str();
  187 + std::cout << std::endl;
  188 + std::cout<<std::endl;
  189 +
  190 + std::cout<<"mie scattering-------"<<std::endl;
  191 + if (spheres.size() == 1) {
  192 + std::cout << "center: " << spheres[0].c << " units" << std::endl;
  193 + std::cout << "radius: " << spheres[0].radius << " units" << std::endl;
  194 + std::cout << "ref. index: " << spheres[0].n << std::endl;
  195 + }
  196 + else {
  197 + for (size_t si = 0; si < spheres.size(); si++) {
  198 + std::cout << "sphere [" << si << "]: r = " << spheres[si].radius << " n = " << spheres[si].n << " c = " << spheres[si].c << std::endl;
  199 + }
  200 + }
  201 + std::cout << "MC samples: " << g_MC << std::endl;
  202 + std::cout<<std::endl;
  203 +}
  204 +
  205 +void load_spheres(std::string filename) {
  206 + stim::table t;
  207 + t.read_ascii(filename, ' ');
  208 +
  209 + std::cout << t.str() << std::endl;
  210 + size_t ns = t.rows();
  211 + spheres.resize(ns);
  212 +
  213 + std::vector< std::vector< double > > s = t.get_vector< double >();
  214 +
  215 +
  216 + for (size_t si = 0; si < ns; si++) { //for each sphere in the file
  217 + spheres[si] = stim::scalarmie<float>(s[si][0], //generate a corresponding scalar mie object in the cluster
  218 + stim::complex<float>(s[si][1], s[si][2]),
  219 + stim::vec3<float>(s[si][3], s[si][4], s[si][5]));
  220 + }
  221 +}
  222 +
  223 +void set_simulation(){
  224 +
  225 + if(args.nargs() > 0) outfile = args.arg(0); //the first argument is the filename (if any)
  226 +
  227 + if (args["cuda"].as_int() == -1) use_cuda = false;
  228 + if (use_cuda) cuda_device = args["cuda"].as_int();
  229 + if (!stim::testDevice(cuda_device, cuda_major, cuda_minor)) { //make sure the device supports the necessary compute capability
  230 + std::cout << "ERROR: CUDA device doesn't support compute capability " << cuda_major << "." << cuda_minor << std::endl;
  231 + exit(1);
  232 + }
  233 + if(args["debug"].is_set()) g_debug_images = true; //set debug flags
  234 + if(args["verbose"].is_set()) g_verbose = true;
  235 +
  236 + //set the condenser NA
  237 + g_NA_cond[1] = args["condenser"].as_float(0); //set the external condenser NA
  238 + if(args["condenser"].nargs() == 2) g_NA_cond[0] = args["condenser"].as_float(1); //if the internal NA is specified, use it
  239 + else g_NA_cond[0] = 0; //otherwise set the internal NA to zero
  240 +
  241 + //set the incident wavelength
  242 + if(args["wavenumber"].is_set()) //if the wavelength is given in terms of wavenumber
  243 + g_lambda = 10000.0 / args["wavenumber"].as_float(0); //convert to wavelength (assume microns)
  244 + else
  245 + g_lambda = args["lambda"].as_float(0); //otherwise assume the wavelength is specified in microns
  246 +
  247 + g_Nl = args["order"].as_int(0); //set the incident field order
  248 +
  249 + //set the mie scattering parameters
  250 + g_MC = args["samples"].as_int(); //get the number of monte-carlo samples
  251 + if (args["spheres"]) {
  252 + load_spheres(args["spheres"].as_string());
  253 + }
  254 + else { //otherwise store the sphere specified at the command line
  255 + float radius = args["radius"].as_float(0); //set the radius of the sphere
  256 + stim::complex<float> n;
  257 + n.real(args["n"].as_float(0)); //set the real part of the refractive index
  258 + if (args["n"].nargs() == 2) //if an imaginary part is specified
  259 + n.imag(args["n"].as_float(1)); //set the imaginary part
  260 + else n.imag(0); //otherwise assume the imaginary part is zero
  261 + stim::vec3<float> c = stim::vec3<float>(args["c"].as_float(0), args["c"].as_float(1), args["c"].as_float(2));
  262 + spheres.resize(1);
  263 + spheres[0] = stim::scalarmie<float>(radius, n, c);
  264 + }
  265 +
  266 + g_padding = args["padding"].as_int(); //get the near field padding factor
  267 + g_R_ff = args["res"].as_int(0);
  268 + g_R_nf = g_R_ff * (2 * g_padding + 1); //set the field resolution to the first resolution argument
  269 + g_S_ff = args["size"].as_int(0);
  270 + g_S_nf = g_S_ff * (2 * g_padding + 1); //get the size of the field slice (assume square at first)
  271 + g_Z_ff = args["z"].as_float(); //get the z-axis position for the desired far-field image
  272 +
  273 + if (args["nobackprop"]) use_backprop = false;
  274 + if(use_backprop && (abs(g_Z_ff) < spheres[0].radius)){ //if the near field slice is cutting through the sphere
  275 + g_backprop = g_Z_ff - spheres[0].radius; //calculate the number of units that the field has to be back-propagated
  276 + g_Z_nf = spheres[0].radius;
  277 + }
  278 + else {
  279 + g_Z_nf = g_Z_ff;
  280 + }
  281 +
  282 + g_NA_obj[1] = args["objective"].as_float(0); //set the external condenser NA
  283 + if(args["objective"].nargs() == 2) g_NA_obj[0] = args["objective"].as_float(1); //if the internal NA is specified, use it
  284 + else g_NA_obj[0] = 0; //otherwise set the internal NA to zero
  285 +
  286 + if (args["simage"].is_set()) //if a source image is provided
  287 + g_sources = simage<float>(g_S_ff, args["simage"].as_string()); //convert the image to a list of sources
  288 + else if (args["sgrid"]) {
  289 + if (args["sgrid"].nargs() == 1)
  290 + g_sources = sgrid<float>(g_S_ff, args["sgrid"].as_int());
  291 + else
  292 + g_sources = sgrid<float>(g_S_ff, args["sgrid"].as_int(0), args["sgrid"].as_int(1));
  293 + }
  294 + else
  295 + g_sources.push_back( stim::vec3<float>(args["f"].as_float(0), args["f"].as_float(1), 0) ); //otherwise just put a single point source at the origin
  296 +}
  297 +
  298 +int main(int argc, char* argv[]){
  299 + set_arguments(argc, argv); //set the input arguments
  300 + set_simulation(); //set all simulation parameters
  301 +
  302 + if(g_verbose) output_simulation(); //output all simulation parameters
  303 +
  304 + stim::vec3<float> d(0, 0, 1); //beam direction
  305 +
  306 + size_t I_bytes = sizeof(float) * g_R_ff * g_R_ff; //number of bytes in the final intensity images
  307 +
  308 + //background intensity image
  309 + float* I0; //allocate space for the background intensity image
  310 + float* I;
  311 +
  312 + stim::scalarfield<float> E0(g_R_nf, g_R_nf, (float)g_S_nf, (float)g_Z_nf + (float)g_backprop); //incident nearfield
  313 + stim::scalarfield<float> E(g_R_nf, g_R_nf, (float)g_S_nf, (float)g_Z_nf); //create the scattered field
  314 + stim::scalarfield<float> Eff(g_R_ff, g_R_ff, (float)g_S_ff, (float)g_Z_ff); //create the far field
  315 +
  316 + float* X = (float*) malloc( E0.size() * sizeof(float) );
  317 + float* Y = (float*) malloc( E0.size() * sizeof(float) );
  318 + float* Z = (float*) malloc( E0.size() * sizeof(float) );
  319 + E.meshgrid(X, Y, Z, stim::CPUmem); //create a meshgrid for the scattering plane (different from E0 b/c we can't propagate the internal field)
  320 +
  321 +
  322 + E0.meshgrid(); //create meshgrids in the field object for field calculations
  323 + E.meshgrid();
  324 +
  325 + if (use_cuda) {
  326 + cudaSetDevice(cuda_device);
  327 + HANDLE_ERROR(cudaMalloc(&I0, I_bytes));
  328 + HANDLE_ERROR(cudaMemset(I0, 0, I_bytes)); //create an array to store the image intensity
  329 + HANDLE_ERROR(cudaMalloc(&I, I_bytes));
  330 + HANDLE_ERROR(cudaMemset(I, 0, I_bytes));
  331 + E0.to_gpu();
  332 + E.to_gpu();
  333 + }
  334 + else {
  335 + I0 = (float*)calloc(I_bytes, 1);
  336 + I = (float*)calloc(I_bytes, 1);
  337 + }
  338 +
  339 + std::chrono::high_resolution_clock::time_point t_begin = std::chrono::high_resolution_clock::now();
  340 + double P = 0;
  341 + std::thread t1(progressbar_thread, &P); //start the progress bar thread
  342 + for(size_t s = 0; s < g_sources.size(); s++){
  343 + stim::scalarbeam<float> beam((float)g_lambda, 1, g_sources[s], d, (float)g_NA_cond[1], (float)g_NA_cond[0]); //create a focused beam
  344 +
  345 + beam.eval(E0, g_Nl); //evaluate the incident field
  346 +
  347 + if(g_debug_images && s == 0){
  348 + stim::filename incident_file = outfile.prefix(outfile.prefix() + "_0_nfi");
  349 + E0.image(incident_file.str(), stim::complexMag);
  350 + }
  351 + E0.crop(g_padding, Eff);
  352 + Eff.intensity(I0);
  353 +
  354 + spheres.eval(E, beam, g_Nl, g_MC);
  355 +
  356 + if(g_debug_images && s == 0){
  357 + stim::filename nf_file = outfile.prefix(outfile.prefix() + "_1_nfs");
  358 + E.image(nf_file.str(), stim::complexMag);
  359 + }
  360 +
  361 + E.propagate(g_backprop, stim::TAU / g_lambda); //propagate the field back to the origin
  362 +
  363 + if(g_debug_images && s == 0){
  364 + stim::filename propagate_file = outfile.prefix(outfile.prefix()+"_2_nfprop");
  365 + E.image(propagate_file.str(), stim::complexMag);
  366 + }
  367 +
  368 + E.crop(g_padding, Eff); //crop out the far field slice
  369 + Eff.intensity(I); //calculate the far field intensity and add it to the single-beam image
  370 +
  371 + if(g_debug_images && s == 0){
  372 + stim::filename cropped_file = outfile.prefix(outfile.prefix() + "_3_cropped");
  373 + Eff.image(cropped_file.str(), stim::complexMag);
  374 + }
  375 + P = (double)(s + 1) / (double)g_sources.size() * 100;
  376 +
  377 + }
  378 + t1.join();
  379 +
  380 + std::chrono::high_resolution_clock::time_point t_end = std::chrono::high_resolution_clock::now();
  381 + std::chrono::duration<double> t_total = (t_end - t_begin);
  382 + std::cout<<"Time: "<<t_total.count()<<" s"<<std::endl;
  383 +
  384 + //save the near-field image
  385 + if (g_debug_images) {
  386 + stim::filename nearfield_real_file = outfile.prefix(outfile.prefix() + "_3_detector_real");
  387 + E.image(nearfield_real_file.str(), stim::complexReal);
  388 + stim::filename nearfield_imag_file = outfile.prefix(outfile.prefix() + "_3_detector_imag");
  389 + E.image(nearfield_imag_file.str(), stim::complexImaginary);
  390 + }
  391 +
  392 + //save the background and single-beam images
  393 + stim::filename background_file = outfile.prefix(outfile.prefix() + "background");
  394 + stim::filename intensity_file = outfile.prefix(outfile.prefix() + "intensity");
  395 +
  396 + if (use_cuda) {
  397 + stim::gpu2image<float>(I0, background_file.str(), g_R_ff, g_R_ff, stim::cmBrewer);
  398 + stim::gpu2image<float>(I, intensity_file.str(), g_R_ff, g_R_ff, stim::cmBrewer);
  399 + }
  400 + else {
  401 + stim::cpu2image<float>(I0, background_file.str(), g_R_ff, g_R_ff, stim::cmBrewer);
  402 + stim::cpu2image<float>(I, intensity_file.str(), g_R_ff, g_R_ff, stim::cmBrewer);
  403 + }
  404 +
  405 +
  406 + float* A;
  407 + HANDLE_ERROR( cudaMalloc(&A, sizeof(float) * g_R_ff * g_R_ff) );
  408 + gpu_absorbance(A, I, I0, g_R_ff * g_R_ff);
  409 + stim::filename absorbance_file = outfile.prefix(outfile.prefix() + "absorbance");
  410 + stim::gpu2image<float>(A, absorbance_file.str(), g_R_ff, g_R_ff, stim::cmBrewer);
  411 +}
0 \ No newline at end of file 412 \ No newline at end of file
sources.h 0 → 100644
  1 +++ a/sources.h
  1 +template<typename T>
  2 +std::vector< stim::vec3<T> > simage(T size, std::string imagefile) {
  3 + std::vector< stim::vec3<T> > sources;
  4 +
  5 + stim::image<unsigned char> src(imagefile);
  6 +
  7 + float xmin = -size / 2;
  8 + float ymin = -size / 2;
  9 +
  10 + size_t xi, yi;
  11 + float xa, ya;
  12 + for (yi = 0; yi < src.height(); yi++) {
  13 + for (xi = 0; xi < src.width(); xi++) {
  14 + if (src(xi, yi)) {
  15 + xa = (float)xi / (float)(src.width() - 1);
  16 + ya = (float)yi / (float)(src.height() - 1);
  17 +
  18 + sources.push_back(stim::vec3<float>(xa * size + xmin, ya * size + ymin, 0));
  19 +
  20 + }
  21 + }
  22 + }
  23 +
  24 + return sources;
  25 +}
  26 +
  27 +template<typename T>
  28 +std::vector< stim::vec3<T> > sgrid(T size, unsigned int n, unsigned int m = 0) {
  29 + std::vector< stim::vec3<T> > sources;
  30 +
  31 + if (m == 0) m = n;
  32 +
  33 + float xmin = -size / 2 + size/(n + 1);
  34 + float ymin = -size / 2 + size/(m + 1);
  35 +
  36 + size_t xi, yi;
  37 + float xa, ya;
  38 + for (yi = 0; yi < m; yi++) {
  39 + for (xi = 0; xi < n; xi++) {
  40 + xa = (float)xi / (float)(n + 1);
  41 + ya = (float)yi / (float)(m + 1);
  42 +
  43 + sources.push_back(stim::vec3<float>(xa * size + xmin, ya * size + ymin, 0));
  44 + }
  45 + }
  46 +
  47 + return sources;
  48 +}
0 \ No newline at end of file 49 \ No newline at end of file