diff --git a/CMakeLists.txt b/CMakeLists.txt index 28fb22a..03f6820 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,11 @@ project(bimsim) #set the module directory set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}") +set(CMAKE_AUTOMOC ON) + +# As moc files are generated in the binary dir, tell CMake + # to always look for includes there: +set(CMAKE_INCLUDE_CURRENT_DIR ON) #find BOOST set(Boost_USE_STATIC_LIBS ON) @@ -12,44 +17,34 @@ set(Boost_USE_MULTITHREADED ON) set(Boost_USE_STATIC_RUNTIME OFF) find_package( Boost 1.46.0 COMPONENTS program_options ) -#find the Qt4 -find_package(Qt4 REQUIRED) +#find the Qt5 +find_package(Qt5Widgets REQUIRED) +find_package(Qt5Core REQUIRED) +find_package(Qt5Gui REQUIRED) +#find_package(Qt5OpenGL REQUIRED) include_directories(${QT_INCLUDE_DIRECTORY}) -include(${QT_USE_FILE}) #set up CUDA find_package(CUDA) -#find OpenGL -#find_package(OpenGL REQUIRED) - -#find GLUT -#set(GLUT_ROOT_PATH $ENV{GLUT_ROOT_PATH}) -#find_package(GLUT REQUIRED) - -#find GLEW -#find_package(GLEW REQUIRED) - -#add Qt OpenGL stuff -#set(QT_USE_QTOPENGL TRUE) - #ask the user for the RTS location -set(RTS_ROOT_PATH $ENV{RTS_ROOT_PATH}) find_package(RTS REQUIRED) #set the include directories include_directories( ${CMAKE_CURRENT_BINARY_DIR} - ${QT_INCLUDES} - ${QT_QTOPENGL_INCLUDE_DIR} -# ${OPENGL_INCLUDE_DIR} -# ${GLEW_INCLUDE_PATH} -# ${GLUT_INCLUDE_DIR} + ${Qt5Widgets_INCLUDE_DIRS} + ${Qt5Core_INCLUDE_DIRS} + ${Qt5Gui_INCLUDE_DIRS} +# ${Qt5OpenGL_INCLUDE_DIRS} ${RTS_INCLUDE_DIR} - ${Boost_INCLUDE_DIR} + ${Boost_INCLUDE_DIR} ) -#enable warnings +#build position independent code for Qt (-fPIC) +#set(CMAKE_CXX_FLAGS "-fPIC") + +#enable warnings (-Wall) if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) add_definitions(-Wall) endif() @@ -57,39 +52,39 @@ endif() #Assign source files to the appropriate variables file(GLOB SRC_CPP "*.cpp") file(GLOB SRC_H "*.h") -#file(GLOB SRC_UI "*.ui") -#file(GLOB SRC_QRC "*.qrc") +file(GLOB SRC_UI "*.ui") file(GLOB SRC_CU "*.cu") -#set up copying data files -#configure_file(rtsSU8_k.txt ${CMAKE_CURRENT_BINARY_DIR}/rtsSU8_k.txt @ONLY) -#configure_file(rtsSU8_n.txt ${CMAKE_CURRENT_BINARY_DIR}/rtsSU8_n.txt @ONLY) -#configure_file(SurfaceMagnitude.glsl ${CMAKE_CURRENT_BINARY_DIR}/SurfaceMagnitude.glsl @ONLY) -#configure_file(SurfaceDeform.glsl ${CMAKE_CURRENT_BINARY_DIR}/SurfaceDeform.glsl @ONLY) - #determine which source files have to be moc'd -#Qt4_wrap_cpp(UI_MOC ${SRC_H}) -#Qt4_wrap_ui(UI_H ${SRC_UI}) -#Qt4_add_resources(ALL_RCC ${ALL_QRC}) +Qt5_wrap_cpp(UI_MOC ${SRC_H}) +Qt5_wrap_ui(UI_H ${SRC_UI}) +Qt5_add_resources(ALL_RCC ${ALL_QRC}) -#moc the necessary files -#Qt4_automoc(${ALL_CPP}) - -#source_group(QtMoc FILES ${UI_MOC}) -#source_group(QtUI FILES ${SRC_UI}) +set(CMAKE_CXX_FLAGS "${Qt5Widgets_EXECUTABLE_COMPILE_FLAGS}") #create an executable cuda_add_executable(bimsim ${SRC_CPP} ${SRC_H} -# ${UI_H} -# ${UI_MOC} -# ${ALL_RCC} + ${RTS_SOURCE} + ${UI_H} + ${SRC_UI} ${SRC_CU} ) +#specify which qt5 modules to use +qt5_use_modules(bimsim Core Widgets OpenGL Gui) + #set the link libraries -target_link_libraries(bimsim ${QT_LIBRARIES} ${QT_QTOPENGL_LIBRARY} ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ${GLEW_LIBRARY} ${CUDA_cufft_LIBRARY} ${CUDA_cublas_LIBRARY} ${GLUT_glut_LIBRARY} ${Boost_LIBRARIES}) +target_link_libraries(bimsim + ${Qt5Widgets_LIBRARIES} + ${Qt5Core_LIBRARIES} + ${Qt5Gui_LIBRARIES} +# ${Qt5OpenGL_LIBRARIES} + ${CUDA_cufft_LIBRARY} + ${CUDA_cublas_LIBRARY} + ${Boost_LIBRARIES} + ) diff --git a/FindGLEW.cmake b/FindGLEW.cmake index 9fe32e8..07b152c 100644 --- a/FindGLEW.cmake +++ b/FindGLEW.cmake @@ -10,15 +10,14 @@ IF (WIN32) FIND_PATH( GLEW_INCLUDE_PATH GL/glew.h $ENV{PROGRAMFILES}/GLEW/include - ${PROJECT_SOURCE_DIR}/src/nvgl/glew/include - DOC "The directory where GL/glew.h resides") + $ENV{GLEW_PATH}/include + ) FIND_LIBRARY( GLEW_LIBRARY NAMES glew GLEW glew32 glew32s PATHS $ENV{PROGRAMFILES}/GLEW/lib - ${PROJECT_SOURCE_DIR}/src/nvgl/glew/bin - ${PROJECT_SOURCE_DIR}/src/nvgl/glew/lib - DOC "The GLEW library") + $ENV{GLEW_PATH}/lib/Release/Win32 + ) ELSE (WIN32) FIND_PATH( GLEW_INCLUDE_PATH GL/glew.h /usr/include @@ -42,10 +41,4 @@ IF (GLEW_INCLUDE_PATH) SET( GLEW_FOUND 1 CACHE STRING "Set to 1 if GLEW is found, 0 otherwise") ELSE (GLEW_INCLUDE_PATH) SET( GLEW_FOUND 0 CACHE STRING "Set to 1 if GLEW is found, 0 otherwise") -ENDIF (GLEW_INCLUDE_PATH) - -MARK_AS_ADVANCED( - GLEW_FOUND - GLEW_INCLUDE_PATH - GLEW_LIBRARY -) \ No newline at end of file +ENDIF (GLEW_INCLUDE_PATH) \ No newline at end of file diff --git a/FindRTS.cmake b/FindRTS.cmake index fbfcc44..0353e71 100644 --- a/FindRTS.cmake +++ b/FindRTS.cmake @@ -1,14 +1,8 @@ # Tries to find the RTS include directory - FIND_PATH( RTS_INCLUDE_DIR NAMES rts_glShaderProgram.h - PATHS - ${CMAKE_CURRENT_SOURCE_DIR}/rts - ${RTS_ROOT_PATH} -) +find_path( RTS_INCLUDE_DIR rts_glShaderProgram.h HINTS $ENV{RTS_PATH}) -IF (RTS_FOUND) - #The following deprecated settings are for backwards compatibility with CMake1.4 - SET (RTS_INCLUDE_PATH ${RTS_INCLUDE_DIR}) -ENDIF(RTS_FOUND) +set (RTS_INCLUDE_DIRS ${RTS_INCLUDE_DIR}) +file (GLOB RTS_SOURCE ${RTS_INCLUDE_DIR}/rts/source/*.cpp) FIND_PACKAGE_HANDLE_STANDARD_ARGS(RTS REQUIRED_VARS TRUE RTS_INCLUDE_DIR) diff --git a/dataTypes.h b/dataTypes.h index ee2efdf..f30bc88 100644 --- a/dataTypes.h +++ b/dataTypes.h @@ -32,10 +32,10 @@ extern bool verbose; #include "rts/math/point.h" #include "rts/math/quad.h" -typedef rts::rtsComplex bsComplex; -typedef rts::rtsVector bsVector; -typedef rts::rtsPoint bsPoint; -typedef rts::rtsQuad bsRect; +typedef rts::complex bsComplex; +typedef rts::vector bsVector; +typedef rts::point bsPoint; +typedef rts::quad bsRect; #endif diff --git a/fileout.cu b/fileout.cu index 37e86b3..a0f935c 100644 --- a/fileout.cu +++ b/fileout.cu @@ -1,66 +1,4 @@ #include "fileout.h" -//#include "scalarfield.h" - - -/*void fileoutStruct::saveMag(fieldslice* U, std::string filename, rts::colormap::colormapType cmap) -{ - int Rx = U->R[0]; - int Ry = U->R[1]; - - //allocate space for a scalar field on the GPU - ptype* gpuScalar; - int memsize = sizeof(ptype) * Rx * Ry; - HANDLE_ERROR(cudaMalloc((void**) &gpuScalar, memsize)); - HANDLE_ERROR(cudaMemset(gpuScalar, 0, memsize)); - U->Mag(gpuScalar); - - - rts::colormap::gpu2image(gpuScalar, filename, Rx, Ry, 0, colorMax, cmap); - - HANDLE_ERROR(cudaFree(gpuScalar)); -} - -void fileoutStruct::saveReal_scalar(fieldslice* U, std::string filename, rts::colormap::colormapType cmap) -{ - //returns the real component - scalarslice sf = U->Real(); - sf.toImage(filename, false, cmap); - -} - -void fileoutStruct::saveImag_scalar(fieldslice* U, std::string filename, rts::colormap::colormapType cmap) -{ - //returns the imaginary component of a field (assumed scalar) - scalarslice sf = U->Imag(); - sf.toImage(filename, false, cmap); -} - -void fileoutStruct::saveIntensity(fieldslice* U, std::string filename, rts::colormap::colormapType cmap) -{ - //get the intensity of the field - scalarslice sf = U->Intensity(); - sf.toImage(filename, true, cmap); -} - -void fileoutStruct::saveAngularSpectrum(fieldslice* U, std::string filename, rts::colormap::colormapType cmap) -{ - ptype* gpuScalar; - int memsize = sizeof(ptype) * U->R[0] * U->R[1]; - HANDLE_ERROR(cudaMalloc((void**) &gpuScalar, memsize)); - HANDLE_ERROR(cudaMemset(gpuScalar, 0, memsize)); - - //convert the field slice to its angular spectrum - U->toAngularSpectrum(); - - //convert the angular spectrum to a scalar field - U->Mag(gpuScalar); - - //save the color image - rts::colormap::gpu2image(gpuScalar, filename, U->R[0], U->R[1], 0, colorMax, cmap); - - HANDLE_ERROR(cudaFree(gpuScalar)); - -}*/ void fileoutStruct::saveNearField(nearfieldStruct* nf) { @@ -221,7 +159,7 @@ void fileoutStruct::Save(microscopeStruct* scope) } - + } diff --git a/fileout.h b/fileout.h index 96d0c70..b1a4fd3 100644 --- a/fileout.h +++ b/fileout.h @@ -5,7 +5,7 @@ //#include "defaults.h" #include "dataTypes.h" -#include "rts/graphics/colormap.h" +#include "rts/visualization/colormap.h" #include "fieldslice.h" #include "nearfield.h" #include "microscope.h" diff --git a/main.cpp b/main.cpp index 9440b8f..dfb4038 100644 --- a/main.cpp +++ b/main.cpp @@ -23,6 +23,12 @@ microscopeStruct* SCOPE; #include "warnings.h" +#include "planewave.h" + +//user interface +#include "qtMainDialog.h" +bool gui = false; + fileoutStruct gFileOut; bool verbose = false; using namespace std; @@ -33,6 +39,15 @@ int cbessjyva(double v,complex z,double &vm,complex*cjv, int main(int argc, char *argv[]) { + //benchtest planewave class + rts::vector k(1, 0, 0); + rts::vector E(2, 2, 0); + planewave P(k, E); + + std::cout<init(); - gFileOut.Save(SCOPE); + if(gui) + { + QApplication app(argc, argv); + qtMainDialog bsDialog; + + //populate the user interface with the default and command-line values + bsDialog.populateUi(); + bsDialog.show(); + return app.exec(); + } + else + { + + gFileOut.Save(SCOPE); - if(verbose) - OutputOptions(); + if(verbose) + OutputOptions(); - SCOPE->destroy(); + SCOPE->destroy(); + } diff --git a/microscope.cu b/microscope.cu index 3c2dfb2..02e1979 100644 --- a/microscope.cu +++ b/microscope.cu @@ -4,7 +4,7 @@ #include "rts/tools/progressbar.h" #include "rts/cuda/timer.h" #include "dataTypes.h" -#include "rts/graphics/colormap.h" +#include "rts/visualization/colormap.h" #include @@ -49,6 +49,8 @@ microscopeStruct::microscopeStruct() void microscopeStruct::init() { + + nf.scalarSim = scalarSim; //Ud.scalarField = scalarSim; //Ufd.scalarField = scalarSim; @@ -137,7 +139,7 @@ void microscopeStruct::clearDetector() //flag for a vector simulation void microscopeStruct::setPos(bsPoint pMin, bsPoint pMax, bsVector normal) { - pos = rts::rtsQuad(pMin, pMax, normal); + pos = rts::quad(pMin, pMax, normal); } void microscopeStruct::setRes(int x_res, int y_res, int pad, int supersampling) diff --git a/microscope.h b/microscope.h index b36317c..3380663 100644 --- a/microscope.h +++ b/microscope.h @@ -31,7 +31,7 @@ struct microscopeStruct ptype objective[2]; //image position and orientation in world space - rts::rtsQuad pos; + rts::quad pos; vector focalPoints; diff --git a/montecarlo.cpp b/montecarlo.cpp index 0690d20..7609d7f 100644 --- a/montecarlo.cpp +++ b/montecarlo.cpp @@ -20,15 +20,15 @@ void mcSampleNA(bsVector* samples, int N, bsVector k, ptype NAin, ptype NAout) //get the axis of rotation for transforming (0, 0, 1) to k //k = -k; ptype cos_angle = k.dot(bsVector(0, 0, 1)); - rts::rtsMatrix rotation; + rts::matrix rotation; if(cos_angle != 1.0) { bsVector axis = bsVector(0, 0, 1).cross(k).norm(); ptype angle = acos(cos_angle); - rts::rtsQuaternion quat; + rts::quaternion quat; quat.CreateRotation(angle, axis); - rotation = quat.toMatrix(); + rotation = quat.toMatrix3(); } //find the phi values associated with the cassegrain ring diff --git a/nearfield.cpp b/nearfield.cpp index 2a32fbb..668e0e4 100644 --- a/nearfield.cpp +++ b/nearfield.cpp @@ -39,7 +39,7 @@ void nearfieldStruct::destroy() void nearfieldStruct::setPos(bsPoint pMin, bsPoint pMax, bsVector normal) { - pos = rts::rtsQuad(pMin, pMax, normal); + pos = rts::quad(pMin, pMax, normal); } void nearfieldStruct::setRes(int x_res, int y_res) @@ -115,7 +115,7 @@ void nearfieldStruct::calcSpheres() //set the refractive index for the sphere int imat = sVector[i].iMaterial; - rts::rtsComplex n = mVector[imat](lambda); + rts::complex n = mVector[imat](lambda); //calculate the scattering coefficients sVector[i].calcCoeff(lambda, n); @@ -137,19 +137,47 @@ void nearfieldStruct::calcSpheres() void nearfieldStruct::calcUs() { - - if(lut_us) - scalarUpLut(); + if(scalarSim) + { + if(lut_us) + scalarUpLut(); + else + scalarUs(); + + if(lut_us) + std::cout<<"Using LUT for Us"< pos; + rts::quad pos; //slices for the focused field fieldslice Uf; @@ -89,9 +92,15 @@ struct nearfieldStruct //this function re-computes the focused field void calcUf(); + + //scalar functions void scalarUf(); void scalarUfLut(); + //vector functions + void vectorUf(); + void vectorUfLut(); + void calcBesselLut(ptype* j, ptype d_min, ptype d_max, int dR); //compute the field scattered by all of the materials diff --git a/nfScalarUf.cu b/nfScalarUf.cu index adc61ab..261c3d8 100644 --- a/nfScalarUf.cu +++ b/nfScalarUf.cu @@ -178,6 +178,8 @@ void nearfieldStruct::scalarUf() else { //pre-compute the cosine of the obscuration and objective angles + //cout<<"Condenser angle in: "< +#include "rts/cuda/error.h" +#include "rts/cuda/timer.h" + +//Incident field for a single plane wave +__global__ void gpuVectorUfp(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR) +{ + /*Compute the scalar focused field using Debye focusing + k = direction of focused light, where |k| = 2*pi/lambda + P = rect struct describing the field slice + rX, rY = resolution of the field slice + cNAin = inner NA of the condenser + cNAout = outer NA of the condenser + */ + + //get the current coordinate in the plane slice + 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 >= uR || iv >= vR) return; + + //compute the index (easier access to the scalar field array) + int i = iv*uR + iu; + + //compute the parameters for u and v + ptype u = (ptype)iu / uR; + ptype v = (ptype)iv / vR; + + //get the rtsPoint in world space and then the r vector + bsPoint p = ABCD(u, v); + bsVector r = p - f; + //ptype d = r.len(); + + ptype k_dot_r = kmag * k.dot(r); + bsComplex d(0, k_dot_r); + + Uf[i] = exp(d) * A; + +} + +//Incident field for a focused point source +__global__ void gpuVectorUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR, ptype cosAlpha, ptype cosBeta, int nl, ptype j_conv = 1.4) +{ + //Compute the scalar focused field using Debye focusing + // k = direction of focused light, where |k| = 2*pi/lambda + // P = rect struct describing the field slice + // rX, rY = resolution of the field slice + // cNAin = inner NA of the condenser + // cNAout = outer NA of the condenser + + + //get the current coordinate in the plane slice + 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 >= uR || iv >= vR) return; + + //compute the index (easier access to the scalar field array) + int i = iv*uR + iu; + + //compute the parameters for u and v + ptype u = (ptype)iu / (uR); + ptype v = (ptype)iv / (vR); + + //get the rtsPoint in world space and then the r vector + bsPoint p = ABCD(u, v); + bsVector r = p - f; + ptype d = r.len(); + if(d < EPSILON_FLOAT) + { + Uf[i] = A * 2 * PI * (cosAlpha - cosBeta); + return; + } + + //get info for the light direction and frequency + //k = k.norm(); + r = r.norm(); + + //compute the imaginary factor i^l + bsComplex im = bsComplex(0, 1); + bsComplex il = bsComplex(1, 0); + + //Bessel and Legendre functions are computed dynamically to save memory + //initialize the Bessel and Legendre functions + ptype j[2]; + ptype kd = kmag * d; + rts::init_sbesselj(kd, j); + + ptype P[2]; + //get the angle between k and r (light direction and position vector) + ptype cosTheta; + cosTheta = k.dot(r); + + //deal with the degenerate case where r == 0 + //if(isnan(cosTheta)) + // cosTheta = 0; + rts::init_legendre(cosTheta, P[0], P[1]); + + //initialize legendre functions for the cassegrain angles + ptype Palpha[3]; + //ptype cosAlpha = cos(asin(cNAin)); + rts::init_legendre(cosAlpha, Palpha[0], Palpha[1]); + Palpha[2] = 1; + + ptype Pbeta[3]; + //ptype cosBeta = cos(asin(cNAout)); + rts::init_legendre(cosBeta, Pbeta[0], Pbeta[1]); + Pbeta[2] = 1; + + //for each order l + bsComplex sumUf(0.0, 0.0); + ptype jl = 0.0; + ptype Pl; + for(int l = 0; l<=nl; l++) + { + + if(l==0) + { + + jl = j[0]; + Pl = P[0]; + } + else if(l==1) + { + jl = j[1]; + Pl = P[1]; + + //adjust the cassegrain Legendre function + Palpha[2] = Palpha[0]; + rts::shift_legendre(l+1, cosAlpha, Palpha[0], Palpha[1]); + Pbeta[2] = Pbeta[0]; + rts::shift_legendre(l+1, cosBeta, Pbeta[0], Pbeta[1]); + } + else + { + rts::shift_sbesselj(l, kd, j);//, j_conv); + rts::shift_legendre(l, cosTheta, P[0], P[1]); + + jl = j[1]; + Pl = P[1]; + + //adjust the cassegrain outer Legendre function + Palpha[2] = Palpha[0]; + rts::shift_legendre(l+1, cosAlpha, Palpha[0], Palpha[1]); + Pbeta[2] = Pbeta[0]; + rts::shift_legendre(l+1, cosBeta, Pbeta[0], Pbeta[1]); + } + + sumUf += il * jl * Pl * (Palpha[1] - Palpha[2] - Pbeta[1] + Pbeta[2]); + + il *= im; + } + + Uf[i] = sumUf * 2 * PI * A; + +} + + +void nearfieldStruct::vectorUf() +{ + + + gpuStartTimer(); + + //create one thread for each pixel of the field slice + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); + dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); + + //if we are computing a plane wave, call the gpuScalarUfp function + if(planeWave) + { + std::cout<<"Calculating vector plane wave..."<>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1]); + } + //otherwise compute the condenser info and create a focused field + else + { + //pre-compute the cosine of the obscuration and objective angles + ptype cosAlpha = cos(asin(condenser[0])); + ptype cosBeta = cos(asin(condenser[1])); + //compute the scalar Uf field (this will be in the x_hat channel of Uf) + gpuVectorUf<<>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1], cosAlpha, cosBeta, m); + } + + t_Uf = gpuStopTimer(); +} diff --git a/nfVectorUfLut.cu b/nfVectorUfLut.cu new file mode 100644 index 0000000..76fbbc2 --- /dev/null +++ b/nfVectorUfLut.cu @@ -0,0 +1,10 @@ +#include "nearfield.h" + +#include "rts/math/legendre.h" +#include "rts/cuda/error.h" +#include "rts/cuda/timer.h" + +void nearfieldStruct::vectorUfLut() +{ + +} diff --git a/options.h b/options.h index 586fabf..3aa53ea 100644 --- a/options.h +++ b/options.h @@ -5,7 +5,7 @@ #include "nearfield.h" #include "microscope.h" -#include "rts/graphics/colormap.h" +#include "rts/visualization/colormap.h" #include "fileout.h" //extern nearfieldStruct* NF; extern microscopeStruct* SCOPE; @@ -24,11 +24,18 @@ using namespace std; namespace po = boost::program_options; extern bool verbose; +extern bool gui; static void lNearfield(po::variables_map vm) { + //test to see if we are running a vector field simulation + bool vectorField = false; + if(vm.count("vector")) + vectorField = true; + SCOPE->scalarSim = !vectorField; + //test to see if we are simulating a plane wave bool planeWave = DEFAULT_PLANEWAVE; if(vm.count("plane-wave")) @@ -176,6 +183,10 @@ void lFlags(po::variables_map vm, po::options_description desc) { SCOPE->nf.lut_uf = true; } + + //gui + if(vm.count("gui")) + gui = true; } void lWavelength(po::variables_map vm) @@ -315,6 +326,8 @@ void lSpheres(po::variables_map vm) { //otherwise output an error cout<<"BIMSIM Error - A material is not loaded for sphere "<nf.sVector[s].iMaterial + 1<nf.mVector.size()<(ifs)), std::istreambuf_iterator()); + //std::string instr((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); //load the list of spheres from a string - rts::material newM; - newM.fromStr(instr, ""); + rts::material newM(filenames[i].c_str()); + //newM.fromStr(instr, ""); SCOPE->nf.mVector.push_back(newM); } } @@ -497,19 +510,24 @@ static void SetOptions(po::options_description &desc) { desc.add_options() ("help", "prints this help") - ("verbose", "verbose output\n") + ("gui", "run using the Qt GUI") + ("verbose", "verbose output\n\nOutput Parameters\n--------------------------") + ("vector", "run a vector field simulation") ("intensity", po::value()->default_value(DEFAULT_INTENSITY_FILE), "output measured intensity (filename)") ("absorbance", po::value()->default_value(DEFAULT_ABSORBANCE_FILE), "output measured absorbance (filename)") ("transmittance", po::value()->default_value(DEFAULT_TRANSMITTANCE_FILE), "output measured transmittance (filename)") ("far-field", po::value()->default_value(DEFAULT_FAR_FILE), "output far-field at detector (filename)") ("near-field", po::value()->default_value(DEFAULT_NEAR_FILE), "output field at focal plane (filename)") - ("extended-source", po::value()->default_value(DEFAULT_EXTENDED_SOURCE), "image of source at focus (filename)\n") + ("extended-source", po::value()->default_value(DEFAULT_EXTENDED_SOURCE), "image of source at focus (filename)") + ("output-type", po::value()->default_value(DEFAULT_FIELD_TYPE), "output field value:\n magnitude, polarization, real, imaginary, angular-spectrum") + ("colormap", po::value()->default_value(DEFAULT_COLORMAP), "colormap: gray, brewer") + ("append", "append result to an existing file\n (binary files only)\n\nSphere Parameters\n--------------------------") ("spheres", po::value< vector >()->multitoken(), "sphere position: x y z a m") ("sphere-file", po::value< vector >()->multitoken(), "sphere file:\n [x y z radius material]") ("materials", po::value< vector >()->multitoken(), "refractive indices as n, k pairs:\n ex. -m n0 k0 n1 k1 n2 k2") - ("material-file", po::value< vector >()->multitoken(), "material file:\n [lambda n k]\n") + ("material-file", po::value< vector >()->multitoken(), "material file:\n [lambda n k]\n\nOptics\n--------------------------") ("lambda", po::value()->default_value(DEFAULT_LAMBDA), "incident wavelength") ("nu", po::value(), "incident frequency (in cm^-1)\n(if specified, lambda is ignored)") @@ -518,7 +536,7 @@ static void SetOptions(po::options_description &desc) ("condenser", po::value< vector >()->multitoken(), "condenser numerical aperature\nA pair of values can be used to specify an inner obscuration: -c NAin NAout") ("objective", po::value< vector >()->multitoken(), "objective numerical aperature\nA pair of values can be used to specify an inner obscuration: -c NAin NAout") ("focus", po::value< vector >()->multitoken(), "focal position for the incident point source\n (default = --focus 0 0 0)") - ("plane-wave", "simulates an incident plane wave\n") + ("plane-wave", "simulates an incident plane wave\n\n\nImaging Parameters\n--------------------------") ("resolution", po::value()->default_value(DEFAULT_SLICE_RES), "resolution of the detector") ("plane-lower-left", po::value< vector >()->multitoken(), "lower-left position of the image plane") @@ -526,7 +544,7 @@ static void SetOptions(po::options_description &desc) ("plane-normal", po::value< vector >()->multitoken(), "normal for the image plane") ("xy", po::value< vector >()->multitoken(), "specify an x-y image plane\n (standard microscope)") ("xz", po::value< vector >()->multitoken(), "specify a x-z image plane\n (cross-section of the focal volume)") - ("yz", po::value< vector >()->multitoken(), "specify a y-z image plane\n (cross-section of the focal volume)\n") + ("yz", po::value< vector >()->multitoken(), "specify a y-z image plane\n (cross-section of the focal volume)\n\nSampling Parameters\n--------------------------") ("samples", po::value()->default_value(DEFAULT_SAMPLES), "Monte-Carlo samples used to compute Us") ("padding", po::value()->default_value(DEFAULT_PADDING), "FFT padding for the objective bandpass") @@ -536,10 +554,6 @@ static void SetOptions(po::options_description &desc) ("recursive", "evaluate all Bessel functions recursively\n") ("recursive-us", "evaluate scattered-field Bessel functions recursively\n") ("lut-uf", "evaluate the focused-field using a look-up table\n") - - ("output-type", po::value()->default_value(DEFAULT_FIELD_TYPE), "output field value:\n magnitude, polarization, real, imaginary, angular-spectrum") - ("colormap", po::value()->default_value(DEFAULT_COLORMAP), "colormap: gray, brewer") - ("append", "append result to an existing file\n (binary files only)") ; } @@ -563,7 +577,6 @@ static void LoadParameters(int argc, char *argv[]) lWavelength(vm); //load materials - //loadMaterials(vm); lMaterials(vm); //load the sphere data @@ -575,19 +588,11 @@ static void LoadParameters(int argc, char *argv[]) //load the position and orientation of the image plane lImagePlane(vm); - //load spheres - //loadSpheres(vm); - - lNearfield(vm); loadOutputParams(vm); - //loadMicroscopeParams(vm); - - //loadSliceParams(vm); - //if an extended source will be used if(vm["extended-source"].as() != "") { diff --git a/planewave.h b/planewave.h new file mode 100644 index 0000000..57437f0 --- /dev/null +++ b/planewave.h @@ -0,0 +1,72 @@ +#ifndef __PLANEWAVE__ +#define __PLANEWAVE__ + +#include +#include + +#include "rts/math/vector.h" + +template +class planewave +{ + rts::vector k; + rts::vector E; + + public: + + //constructor, initialize to an x-polarized wave propagating along z + planewave() + { + k = rts::vector(0, 0, 1); + E = rts::vector(1, 0, 0); + } + + planewave(rts::vector k_vec, rts::vector E0) + { + k = k_vec; + + //enforce k \dot E = 0 + rts::vector s = E0.cross(k); + rts::vector E_hat; + + if(s.len() == 0) + E_hat = rts::vector(0, 0, 0); + else + E_hat = (s.cross(k)).norm(); + + E = E_hat * (E_hat.dot(E0)); + } + + // refract will bend the wave vector k to correspond to the normalized vector v + refract(rts::vector v) + { + //make sure that v is normalized + v = v.norm(); + + + } + + std::string toStr() + { + std::stringstream ss; + + ss<<"k = "< + +qtMainDialog::qtMainDialog(QWidget *parent, Qt::WindowFlags flags) + : QMainWindow(parent, flags) +{ + ui.setupUi(this); + + outfile = "ui-out.bmp"; +} + +qtMainDialog::~qtMainDialog() +{ + updating = false; +} + +void qtMainDialog::closeEvent(QCloseEvent *event) +{ + std::cout<<"Exiting"<nf.pos.n(); + + ui.spinNx->setValue(n[0]); + ui.spinNy->setValue(n[1]); + ui.spinNz->setValue(n[2]); + + //get the center point for the image plane + bsPoint c = SCOPE->nf.pos.p(0.5, 0.5); + + ui.spinCx->setValue(c[0]); + ui.spinCy->setValue(c[1]); + ui.spinCz->setValue(c[2]); + + //get the plane size (in microns) + ptype S = SCOPE->nf.pos.X.len(); + ptype pad_div = SCOPE->padding * 2 + 1; + S /= pad_div; + + ui.spinS->setValue(S); + + //get the detector resolution + ui.spinR->setValue(SCOPE->Ud.R[0]); + + updating = false; +} + +/********************************* +Change the image plane position +*********************************/ +void qtMainDialog::positionImage() +{ + if(updating) return; + + //get the plane normal + bsVector n(ui.spinNx->value(), ui.spinNy->value(), ui.spinNz->value()); + + //get the plane center point + bsPoint c(ui.spinCx->value(), ui.spinCy->value(), ui.spinCz->value()); + + //get the plane orientation + ptype theta = ui.spinTheta->value(); + + //get the plane size + ptype S = ui.spinS->value() * (2 * SCOPE->padding + 1); + + //create a new image plane + SCOPE->nf.pos = rts::quad(c, n, S, S, theta); + + +} + +/********************************* +Render an image +*********************************/ +void qtMainDialog::renderImage() +{ + //run the near-field simulation + SCOPE->SimulateScattering(); + + + //determine the colormap type + rts::colormapType cmap = rts::cmGrayscale; + + if( ui.cmbColormap->currentText() == tr("brewer") ) + cmap = rts::cmBrewer; + else + cmap = rts::cmGrayscale; + + //near field rendering + if( ui.radDisplayNearfield->isChecked() ) + { + scalarslice S; + bool positive_vals = false; + + if( ui.cmbDisplayNearfield->currentText() == tr("magnitude") ) + { + std::cout<<"magnitude"<nf.U.Mag(); + positive_vals = true; + //S.toImage(outfile.toStdString(), positive_vals, cmap); + } + else if( ui.cmbDisplayNearfield->currentText() == tr("real") ) + S = SCOPE->nf.U.Real(); + else if( ui.cmbDisplayNearfield->currentText() == tr("imaginary") ) + S = SCOPE->nf.U.Imag(); + + S.toImage(outfile.toStdString(), positive_vals, cmap); + } + + //run the far-field simulation + SCOPE->SimulateImaging(); + + //far field rendering + if( ui.radDisplayFarfield->isChecked() ) + { + scalarslice S; + bool positive_vals = false; + + if( ui.cmbDisplayFarfield->currentText() == tr("magnitude") ) + { + S = SCOPE->Ud.Mag(); + positive_vals = true; + } + else if( ui.cmbDisplayFarfield->currentText() == tr("real") ) + S = SCOPE->Ud.Real(); + else if( ui.cmbDisplayFarfield->currentText() == tr("imaginary") ) + S = SCOPE->Ud.Imag(); + + S.toImage(outfile.toStdString(), positive_vals, cmap); + } + + //detector rendering + if( ui.radDisplayDetector->isChecked() ) + { + scalarslice I; + bool positive_vals = true; + + if( ui.cmbDisplayDetector->currentText() == tr("intensity") ) + I = SCOPE->getIntensity(); + + else if( ui.cmbDisplayDetector->currentText() == tr("absorbance") ) + I = SCOPE->getAbsorbance(); + + I.toImage(outfile.toStdString(), positive_vals, cmap); + } + +} diff --git a/qtMainDialog.h b/qtMainDialog.h new file mode 100755 index 0000000..9a2af07 --- /dev/null +++ b/qtMainDialog.h @@ -0,0 +1,98 @@ +#ifndef VolumeSpiderDialog_H +#define VolumeSpiderDialog_H + +#include +#include +#include +#include +#include +#include "ui_qtMainDialog.h" + +//simulation parameters +#include "microscope.h" +extern microscopeStruct* SCOPE; + +//#include "fileout.h" +//extern fileoutStruct gFileOut; +#include "rts/visualization/colormap.h" + + +class qtMainDialog : public QMainWindow +{ +Q_OBJECT + +public: +qtMainDialog(QWidget *parent = 0, Qt::WindowFlags flags = 0); +~qtMainDialog(); + +void closeEvent(QCloseEvent *event); +bool updating; + +QString outfile; + +void refreshUI() +{ + updating = false; +} + +void populateUi(); +void renderImage(); +void positionImage(); + +private: +Ui::qtMainDialogUI ui; + +public slots: + +/************* +Buttons +*************/ +void on_btnRender_pressed() +{ + renderImage(); +} +void on_btnClearDetector_pressed() +{ + SCOPE->clearDetector(); +} + +/*************** +Spinners +***************/ +void on_spinCx_valueChanged(double d) +{ + positionImage(); +} +void on_spinCy_valueChanged(double d) +{ + positionImage(); +} +void on_spinCz_valueChanged(double d) +{ + positionImage(); +} +void on_spinNx_valueChanged(double d) +{ + positionImage(); +} +void on_spinNy_valueChanged(double d) +{ + positionImage(); +} +void on_spinNz_valueChanged(double d) +{ + positionImage(); +} +void on_spinS_valueChanged(double d) +{ + positionImage(); +} +void on_spinTheta_valueChanged(double d) +{ + positionImage(); +} + + +}; + +#endif // INTERACTIVEMIE_H diff --git a/qtMainDialog.ui b/qtMainDialog.ui new file mode 100644 index 0000000..2f26af4 --- /dev/null +++ b/qtMainDialog.ui @@ -0,0 +1,424 @@ + + + qtMainDialogUI + + + + 0 + 0 + 800 + 600 + + + + MainWindow + + + + + + 80 + 40 + 81 + 27 + + + + -99.989999999999995 + + + 0.010000000000000 + + + + + + 170 + 40 + 81 + 27 + + + + -99.989999999999995 + + + 0.010000000000000 + + + + + + 260 + 40 + 81 + 27 + + + + -99.989999999999995 + + + 0.010000000000000 + + + + + + 170 + 70 + 81 + 27 + + + + -99.989999999999995 + + + 0.010000000000000 + + + + + + 260 + 70 + 81 + 27 + + + + -99.989999999999995 + + + 0.010000000000000 + + + + + + 80 + 70 + 81 + 27 + + + + -99.989999999999995 + + + 0.010000000000000 + + + + + + 10 + 40 + 66 + 17 + + + + pc = + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 10 + 70 + 66 + 17 + + + + n = + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 80 + 130 + 81 + 27 + + + + 0.010000000000000 + + + + + + 10 + 130 + 66 + 17 + + + + S = + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 70 + 210 + 331 + 161 + + + + Display + + + + + 113 + 30 + 116 + 22 + + + + near field + + + true + + + + + + 113 + 60 + 116 + 22 + + + + far field + + + + + + 10 + 30 + 101 + 27 + + + + + magnitude + + + + + real + + + + + imaginary + + + + + + + 10 + 90 + 101 + 27 + + + + + intensity + + + + + absorbance + + + + + + + 113 + 90 + 116 + 22 + + + + detector + + + + + + 10 + 60 + 101 + 27 + + + + + magnitude + + + + + real + + + + + imaginary + + + + + + + 220 + 30 + 91 + 27 + + + + + brewer + + + + + grayscale + + + + + + + 10 + 130 + 98 + 27 + + + + Render + + + + + + 220 + 90 + 51 + 27 + + + + Clear + + + + + + + 260 + 130 + 81 + 27 + + + + 999999 + + + 0 + + + + + + 190 + 130 + 66 + 17 + + + + R = + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + 80 + 100 + 81 + 27 + + + + 6.280000000000000 + + + 0.010000000000000 + + + + + + 10 + 100 + 66 + 17 + + + + theta = + + + Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter + + + + + + + 0 + 0 + 800 + 25 + + + + + + + + diff --git a/scalarslice.cu b/scalarslice.cu index c6f215a..8c129e5 100644 --- a/scalarslice.cu +++ b/scalarslice.cu @@ -11,20 +11,50 @@ scalarslice::scalarslice(int x, int y) R[1] = y; //allocate memory on the GPU - HANDLE_ERROR(cudaMalloc( (void**)&S, sizeof(ptype) * x * y )); + HANDLE_ERROR(cudaMalloc( (void**)&S, sizeof(ptype) * x * y )); + //std::cout<<"Scalerslice created."< #include -using namespace rts; +//using namespace rts; using namespace std; int cbessjyva(double v,complex z,double &vm,complex*cjv, @@ -18,7 +18,7 @@ int cbessjyva_sph(int v,complex z,double &vm,complex*cjv, int bessjyv_sph(int v, double z, double &vm, double* cjv, double* cyv, double* cjvp, double* cyvp); -void sphere::calcCoeff(ptype lambda, rtsComplex ri) +void sphere::calcCoeff(ptype lambda, bsComplex ri) { /* These calculations are done at high-precision on the CPU since they are only required once for each wavelength. @@ -219,7 +219,7 @@ void sphere::calcLut(bsComplex* j, bsComplex* h, ptype lambda, bsComplex n, int calcHankelLut(h, k, rR); } -void sphere::calcUp(ptype lambda, bsComplex n, rts::rtsQuad nfPlane, unsigned int R) +void sphere::calcUp(ptype lambda, bsComplex n, rts::quad nfPlane, unsigned int R) { //calculate the parameters of the lookup table diff --git a/sphere.h b/sphere.h index e3aebe7..7a3bf8a 100644 --- a/sphere.h +++ b/sphere.h @@ -52,7 +52,7 @@ struct sphere //assignment operator sphere & operator=(const sphere &rhs); - + //copy constructor sphere(const sphere &rhs); @@ -81,7 +81,7 @@ struct sphere void calcHankelLut(bsComplex* h, ptype k, int rR); //calculate the scattering domain Us(theta, r) - void calcUp(ptype lambda, bsComplex n, rts::rtsQuad nfPlane, unsigned int R); + void calcUp(ptype lambda, bsComplex n, rts::quad nfPlane, unsigned int R); void scalarUsp(bsComplex* h, int rR, int thetaR); void scalarUip(bsComplex* j, int aR, int thetaR); -- libgit2 0.21.4