diff --git a/ext/ghoul b/ext/ghoul index b922d28993..7fca620ada 160000 --- a/ext/ghoul +++ b/ext/ghoul @@ -1 +1 @@ -Subproject commit b922d289938a14c4be9e64e3955c955541de6e5f +Subproject commit 7fca620ada5b2a8ab0f6ead12e8a00d87837a939 diff --git a/include/openspace/rendering/renderablevolumeexpert.h b/include/openspace/rendering/renderablevolumeexpert.h index 1b274f293a..5d445a8abb 100644 --- a/include/openspace/rendering/renderablevolumeexpert.h +++ b/include/openspace/rendering/renderablevolumeexpert.h @@ -97,6 +97,8 @@ private: ghoul::opencl::CLWorkSize* _ws; ghoul::filesystem::File* _kernelSourceFile; std::vector > _kernelOptions; + std::vector _kernelIncludes; + std::vector > _kernelDefinitions; bool _programUpdateOnSave; // mutexes to prevent inconsistencies @@ -107,6 +109,7 @@ private: GLuint _screenQuad; VolumeRaycasterBox* _colorBoxRenderer; + glm::vec3 _boxScaling; }; diff --git a/include/openspace/util/kameleonwrapper.h b/include/openspace/util/kameleonwrapper.h index 467373ac25..3ba0afbbb0 100644 --- a/include/openspace/util/kameleonwrapper.h +++ b/include/openspace/util/kameleonwrapper.h @@ -49,6 +49,7 @@ public: private: ccmc::Model* _model; + Model _type; ccmc::Interpolator* _interpolator; }; diff --git a/kernels/volume_helpers.cl b/kernels/volume_helpers.cl new file mode 100644 index 0000000000..90ad1f4d9a --- /dev/null +++ b/kernels/volume_helpers.cl @@ -0,0 +1,26 @@ + +float3 CartesianToSpherical(float3 _cartesian); + +float intensityNormalizer(float intensity, float iMin, float iMax); + +float3 CartesianToSpherical(float3 _cartesian) { + // Put cartesian in [-1..1] range first + _cartesian = (float3)(-1.0) + 2.0f* _cartesian; + float r = length(_cartesian); + float theta, phi; + if (r == 0.0) { + theta = phi = 0.0; + } else { + theta = acospi(_cartesian.z/r); + phi = (M_PI + atan2(_cartesian.y, _cartesian.x)) / (2.0*M_PI); + } + r = r / native_sqrt(3.0f); + // Sampler ignores w component + return (float3)(r, theta, phi); +} + +float intensityNormalizer(float intensity, float iMin, float iMax) { + float i = clamp(intensity, iMin, iMax); + i = (i - iMin) / (iMax - iMin); + return clamp(i, 0.0f, 1.0f); +} \ No newline at end of file diff --git a/openspace-data b/openspace-data index 23fa37b495..4b77a977d7 160000 --- a/openspace-data +++ b/openspace-data @@ -1 +1 @@ -Subproject commit 23fa37b495961c714cc361c0f9e5f37326b6d64f +Subproject commit 4b77a977d78069d4472a7d1a1b3bcb2012057595 diff --git a/src/rendering/renderablevolumeexpert.cpp b/src/rendering/renderablevolumeexpert.cpp index de0fa16c61..7c2556365d 100644 --- a/src/rendering/renderablevolumeexpert.cpp +++ b/src/rendering/renderablevolumeexpert.cpp @@ -64,7 +64,8 @@ RenderableVolumeExpert::RenderableVolumeExpert(const ghoul::Dictionary& dictiona RenderableVolume(dictionary), _output(nullptr), _clBackTexture(0), _clFrontTexture(0), _clOutput(0), - _kernelSourceFile(nullptr), _programUpdateOnSave(false), _colorBoxRenderer(nullptr) { + _kernelSourceFile(nullptr), _programUpdateOnSave(false), _colorBoxRenderer(nullptr), + _boxScaling(1.0,1.0,1.0) { _kernelLock = new std::mutex; _textureLock = new std::mutex; @@ -161,6 +162,30 @@ RenderableVolumeExpert::RenderableVolumeExpert(const ghoul::Dictionary& dictiona } } } + + ghoul::Dictionary includeDictionary; + if (kernelDictionary.hasKey("Includes") && kernelDictionary.getValue("Includes", includeDictionary)) { + auto keys = includeDictionary.keys(); + for(auto key: keys) { + std::string includePath; + if(includeDictionary.getValue(key, includePath)) { + if(FileSys.directoryExists(includePath)) { + _kernelIncludes.push_back(absPath(includePath)); + } + } + } + } + + ghoul::Dictionary defineDictionary; + if (kernelDictionary.hasKey("Definitions") && kernelDictionary.getValue("Definitions", defineDictionary)) { + auto keys = defineDictionary.keys(); + for(auto key: keys) { + std::string defintion; + if(defineDictionary.getValue(key, defintion)) { + _kernelDefinitions.push_back(std::make_pair(key, defintion)); + } + } + } } } if (kernelPath != "") { @@ -168,6 +193,22 @@ RenderableVolumeExpert::RenderableVolumeExpert(const ghoul::Dictionary& dictiona } _colorBoxRenderer = new VolumeRaycasterBox(); + double tempValue; + if(dictionary.hasKey("BoxScaling.1") && dictionary.getValue("BoxScaling.1", tempValue)) { + if(tempValue > 0.0) { + _boxScaling[0] = tempValue; + } + } + if(dictionary.hasKey("BoxScaling.2") && dictionary.getValue("BoxScaling.2", tempValue)) { + if(tempValue > 0.0) { + _boxScaling[1] = tempValue; + } + } + if(dictionary.hasKey("BoxScaling.3") && dictionary.getValue("BoxScaling.3", tempValue)) { + if(tempValue > 0.0) { + _boxScaling[2] = tempValue; + } + } } RenderableVolumeExpert::~RenderableVolumeExpert() { @@ -207,6 +248,8 @@ bool RenderableVolumeExpert::initialize() { ghoul::opengl::Texture* volume = loadVolume(_volumePaths.at(i), _volumeHints.at(i)); if(volume) { volume->uploadTexture(); + + LDEBUG("Creating CL texture from GL texture with path '" << _volumePaths.at(i) << "'"); cl_mem volumeTexture = _context.createTextureFromGLTexture(CL_MEM_READ_ONLY, *volume); _volumes.push_back(volume); @@ -304,7 +347,7 @@ void RenderableVolumeExpert::render(const Camera *camera, const psc &thisPositio double factor = pow(10.0,thisPosition[3]); transform = glm::translate(transform, glm::vec3(thisPosition[0]*factor, thisPosition[1]*factor, thisPosition[2]*factor)); transform = glm::rotate(transform, time*speed, glm::vec3(0.0f, 1.0f, 0.0f)); - + transform = glm::scale(transform, _boxScaling); _colorBoxRenderer->render(transform); _textureLock->lock(); @@ -354,6 +397,14 @@ void RenderableVolumeExpert::safeKernelCompilation() { for(auto option: _kernelOptions) { tmpProgram.setOption(option.first, option.second); } + + for(auto defintion: _kernelDefinitions) { + tmpProgram.addDefinition(defintion.first, defintion.second); + } + + // add the include directories + tmpProgram.addIncludeDirectory(_kernelIncludes); + if(tmpProgram.build()) { ghoul::opencl::CLKernel tmpKernel = tmpProgram.createKernel("volumeraycaster"); if(tmpKernel.isValidKernel()) { diff --git a/src/util/kameleonwrapper.cpp b/src/util/kameleonwrapper.cpp index c7526be5e3..d44370676b 100644 --- a/src/util/kameleonwrapper.cpp +++ b/src/util/kameleonwrapper.cpp @@ -30,12 +30,14 @@ #include #include +#include + namespace openspace { std::string _loggerCat = "KameleonWrapper"; -KameleonWrapper::KameleonWrapper(const std::string& filename, Model model) { - switch (model) { +KameleonWrapper::KameleonWrapper(const std::string& filename, Model model): _type(model) { + switch (_type) { case Model::BATSRUS: _model = new ccmc::BATSRUS(); if(!_model) LERROR("BATSRUS:Failed to create model instance"); @@ -53,7 +55,7 @@ KameleonWrapper::KameleonWrapper(const std::string& filename, Model model) { if (!_interpolator) LERROR("Failed to create interpolator"); break; default: - LERROR("Only the BATSRUS model is supported for now. Sorry."); + LERROR("No valid model type provided!"); } } @@ -65,47 +67,177 @@ KameleonWrapper::~KameleonWrapper() { float* KameleonWrapper::getUniformSampledValues(const std::string& var, glm::size3_t outDimensions) { assert(_model && _interpolator); assert(outDimensions.x > 0 && outDimensions.y > 0 && outDimensions.z > 0); - LDEBUG("getUniformSampledValues"); + assert(_type == Model::ENLIL || _type == Model::BATSRUS); + LINFO("Loading CDF data"); int size = outDimensions.x*outDimensions.y*outDimensions.z; float* data = new float[size]; + + // get the grid system string + std::string gridSystem = _model->getGlobalAttribute("grid_system_1").getAttributeString(); + + // remove leading and trailing brackets + gridSystem = gridSystem.substr(1,gridSystem.length()-2); + + // remove all whitespaces + gridSystem.erase(remove_if(gridSystem.begin(), gridSystem.end(), isspace), gridSystem.end()); + + // replace all comma signs with whitespaces + std::replace( gridSystem.begin(), gridSystem.end(), ',', ' '); + + // tokenize + std::istringstream iss(gridSystem); + std::vector tokens{std::istream_iterator{iss},std::istream_iterator{}}; + + // validate + if (tokens.size() != 3) { + LERROR("Something went wrong"); + delete[] data; + return 0; + } - // TODO Check which coordinate system the model use. Currently assumes {x,y,z} - float xMin = _model->getVariableAttribute("x", "actual_min").getAttributeFloat(); - float xMax = _model->getVariableAttribute("x", "actual_max").getAttributeFloat(); - float yMin = _model->getVariableAttribute("y", "actual_min").getAttributeFloat(); - float yMax = _model->getVariableAttribute("y", "actual_max").getAttributeFloat(); - float zMin = _model->getVariableAttribute("z", "actual_min").getAttributeFloat(); - float zMax = _model->getVariableAttribute("z", "actual_max").getAttributeFloat(); - float varMin = _model->getVariableAttribute(var, "actual_min").getAttributeFloat(); - float varMax = _model->getVariableAttribute(var, "actual_max").getAttributeFloat(); + std::string v_x = tokens.at(0), v_y = tokens.at(1), v_z = tokens.at(2); + /* + for(auto t: tokens) + LDEBUG("t: " << t); + */ + /* + LERROR("getVariableAttributeNames"); + std::vector attributeNames = _model->getVariableAttributeNames(); + for(auto name : attributeNames) + LDEBUG(name); + */ + //_model->getVa + + //auto fan = std::find(attributeNames.begin(), attributeNames.end(), ""); + + + //KameleonWrapper (Debug) grid_system_1 + //KameleonWrapper (Debug) grid_1_type + + LDEBUG("Using coordinate system: " << v_x << ", " << v_y << ", " << v_z); + + float xMin = _model->getVariableAttribute(v_x, "actual_min").getAttributeFloat(); + float xMax = _model->getVariableAttribute(v_x, "actual_max").getAttributeFloat(); + float yMin = _model->getVariableAttribute(v_y, "actual_min").getAttributeFloat(); + float yMax = _model->getVariableAttribute(v_y, "actual_max").getAttributeFloat(); + float zMin = _model->getVariableAttribute(v_z, "actual_min").getAttributeFloat(); + float zMax = _model->getVariableAttribute(v_z, "actual_max").getAttributeFloat(); + float varMin = _model->getVariableAttribute(var, "actual_min").getAttributeFloat(); + float varMax = _model->getVariableAttribute(var, "actual_max").getAttributeFloat(); - float stepX = (xMax-xMin)/((float)outDimensions.x-1.0); - float stepY = (yMax-yMin)/((float)outDimensions.y-1.0); - float stepZ = (zMax-zMin)/((float)outDimensions.z-1.0); - - // Temporary hack to keep data proportions correct, will give a lot of empty - // voxels in y and z. TODO Add spacing/voxel dimensions - float step = std::max(stepX, stepY); - step = std::max(step, stepZ); - - for (int x = 0; x < outDimensions.x; ++x) { - unsigned int progress = (unsigned int)(((float)x/(float)outDimensions.x)*100.f); - if (progress % 10 == 0) { - std::cout << "Getting data from kameleon: "<< progress << "% \r" << std::flush; + float stepX = (xMax-xMin)/(static_cast(outDimensions.x)); + float stepY = (yMax-yMin)/(static_cast(outDimensions.y)); + float stepZ = (zMax-zMin)/(static_cast(outDimensions.z)); + + + LDEBUG(v_x << "Min: " << xMin); + LDEBUG(v_x << "Max: " << xMax); + LDEBUG(v_y << "Min: " << yMin); + LDEBUG(v_y << "Max: " << yMax); + LDEBUG(v_z << "Min: " << zMin); + LDEBUG(v_z << "Max: " << zMax); + LDEBUG(var << "Min: " << varMin); + LDEBUG(var << "Max: " << varMax); + + int barWidth = 70; + int lastiProgress = -1; + for (int x = 0; x < outDimensions.x; ++x) { + float progress = static_cast(x) / static_cast(outDimensions.x-1); + int iprogress = static_cast(progress*100.0f); + if (iprogress != lastiProgress) { + + int pos = barWidth * progress; + int eqWidth = pos+1; + int spWidth = barWidth - pos + 2; + std::cout << "[" << std::setfill('=') << std::setw(eqWidth) + << ">" << std::setfill(' ') << std::setw(spWidth) + << "] " << iprogress << " % \r" << std::flush; } + lastiProgress = iprogress; + for (int y = 0; y < outDimensions.y; ++y) { for (int z = 0; z < outDimensions.z; ++z) { - - float xPos = xMin + step*x; - float yPos = yMin + step*y; - float zPos = zMin + step*z; - int index = x + y*outDimensions.x + z*outDimensions.x*outDimensions.y; - // get interpolated data value for (xPos, yPos, zPos) and scale to [0,1] - data[index] = (_interpolator->interpolate(var, xPos, yPos, zPos)-varMin)/(varMax-varMin); + + int index = x + y*outDimensions.x + z*outDimensions.x*outDimensions.y; + + if(_type == Model::BATSRUS) { + float xPos = xMin + stepX*x; + float yPos = yMin + stepY*y; + float zPos = zMin + stepZ*z; + + // get interpolated data value for (xPos, yPos, zPos) + float value = _interpolator->interpolate(var, xPos, yPos, zPos); + + // scale to [0,1] + data[index] = (value-varMin)/(varMax-varMin); + } else if (_type == Model::ENLIL) { + //LDEBUG("data: " << theval); + + // Calculate array index + //unsigned int index = r + theta*xDim_ + phi*xDim_*yDim_; + + // Put r in the [0..sqrt(3)] range + float rNorm = sqrt(3.0)*(float)x/(float)(outDimensions.x-1); + + // Put theta in the [0..PI] range + float thetaNorm = M_PI*(float)y/(float)(outDimensions.y-1); + + // Put phi in the [0..2PI] range + float phiNorm = 2.0*M_PI*(float)z/(float)(outDimensions.z-1); + + // Go to physical coordinates before sampling + float rPh = xMin + rNorm*(xMax-xMin); + float thetaPh = thetaNorm; + //phi range needs to be mapped to the slightly different + // model range to avoid gaps in the data + // Subtract a small term to avoid rounding errors when comparing + // to phiMax. + float phiPh = zMin + phiNorm/(2.0*M_PI)*(zMax-zMin-0.000001); + + // Hardcoded variables (rho or rho - rho_back) + // TODO Don't hardcode, make more flexible + float varValue = 0.f;//, rho_back = 0.f, diff = 0.f; + // See if sample point is inside domain + if (rPh < xMin || rPh > xMax || thetaPh < yMin || + thetaPh > yMax || phiPh < zMin || phiPh > zMax) { + if (phiPh > zMax) { + std::cout << "Warning: There might be a gap in the data\n"; + } + // Leave values at zero if outside domain + } else { // if inside + + // ENLIL CDF specific hacks! + // Convert from meters to AU for interpolator + rPh /= ccmc::constants::AU_in_meters; + // Convert from colatitude [0, pi] rad to latitude [-90, 90] degrees + thetaPh = -thetaPh*180.f/M_PI+90.f; + // Convert from [0, 2pi] rad to [0, 360] degrees + phiPh = phiPh*180.f/M_PI; + // Sample + varValue = _interpolator->interpolate(var, rPh, thetaPh, phiPh); + //rho_back = _interpolator->interpolate("rho-back",rPh,thetaPh,phiPh); + + // Calculate difference (or just rho) + //diff = rho; + //diff = rho - rho_back; + + // Clamp to 0 + //if (diff < 0.f) diff = 0.f; + } + //if(var < 0.0f) var = 0.0f; + //data[index] = var; + data[index] = (varValue-varMin)/(varMax-varMin); + //LDEBUG("varValue:" << varValue); + //LDEBUG("data[index]:" << data[index]); + //data[index] = var; + //data[index] = diff; + } } } } + std::cout << std::endl; + LINFO("Done!"); return data; }