bring in support for multiresolution volume rendering

data/scene/enlilnh/enlilnh.mod

@@ -0,0 +1,22 @@
+return {
+    -- Volume module
+    {   
+        Name = "Enlil New Horizons",
+        Parent = "Root",
+        Ephemeris = {
+            Type = "Static",
+            Position = { 0.0, 0.0, 0.0, 0}
+        },
+        Renderable = {
+            Type = "RenderableMultiresVolume",
+            Translation = {0, 0, 0},
+            Rotation = {2.1, 0, 0},
+            Scaling = {1.1, 1.1, 1.1},
+            ScalingExponent = 12,
+			Source = "tsp/enlil_nh_128_128_16.tsp",
+            TransferFunction = "transferfunctions/fire.txt",
+			BrickSelector = "tf",
+        },
+        GuiName = "/Volumes/ENLIL New Horizons"
+    }
+}

data/scene/enlilnh/transferfunctions/fire.txt

@@ -0,0 +1,10 @@
+width 1024
+lower 0.0
+upper 1.0
+
+mappingkey	0.05	254	0	0	0
+mappingkey	0.10	254	0	0	110
+mappingkey	0.12	254	220	0	254
+mappingkey	0.15	254	100	0	254 
+mappingkey	0.20	180	0	0	254
+mappingkey	0.25	180	0	80	254

include/openspace/rendering/transferfunction.h

@@ -0,0 +1,70 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2015                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __TRANSFERFUNCTION_H__
+#define __TRANSFERFUNCTION_H__
+
+#include <string>
+#include <glm/glm.hpp>
+#include <functional>
+#include <ghoul/opengl/texture.h>
+#include <ghoul/filesystem/file.h>
+#include <memory>
+
+namespace openspace {
+
+    class TransferFunction {
+    public:
+        typedef std::function<void (const TransferFunction&)> TfChangedCallback;
+
+        TransferFunction(const std::string& filepath, TfChangedCallback tfChangedCallback = TfChangedCallback());
+        void setPath(const std::string& filepath);
+        ghoul::opengl::Texture& getTexture();
+        void update();
+        glm::vec4 sample(size_t t);
+        size_t width();
+        void setCallback(TfChangedCallback callback);
+    private:
+        void setTextureFromTxt();
+        void setTextureFromImage();
+        void uploadTexture();
+
+        std::string _filepath;
+        std::unique_ptr<ghoul::filesystem::File> _file = nullptr;
+        std::unique_ptr<ghoul::opengl::Texture> _texture = nullptr;
+        bool _needsUpdate = false;
+        TfChangedCallback _tfChangedCallback;
+    };
+
+    struct MappingKey {
+        float position{0.0f};
+        glm::vec4 color{0.0f,0.0f,0.0f,0.0f};
+        MappingKey(float p, const glm::vec4& c): position(p), color(c) {};
+        MappingKey(float p): position(p), color(glm::vec4(0.0f)) {};
+        bool operator<(const MappingKey& rhs) {return position < rhs.position;};
+    };
+} // namespace openspace
+
+#endif
+

modules/multiresvolume/CMakeLists.txt

@@ -0,0 +1,68 @@
+#########################################################################################
+#                                                                                       #
+# OpenSpace                                                                             #
+#                                                                                       #
+# Copyright (c) 2014-2015                                                               #
+#                                                                                       #
+# Permission is hereby granted, free of charge, to any person obtaining a copy of this  #
+# software and associated documentation files (the "Software"), to deal in the Software #
+# without restriction, including without limitation the rights to use, copy, modify,    #
+# merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    #
+# permit persons to whom the Software is furnished to do so, subject to the following   #
+# conditions:                                                                           #
+#                                                                                       #
+# The above copyright notice and this permission notice shall be included in all copies #
+# or substantial portions of the Software.                                              #
+#                                                                                       #
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   #
+# INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         #
+# PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    #
+# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  #
+# CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  #
+# OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         #
+#########################################################################################
+
+include(${OPENSPACE_CMAKE_EXT_DIR}/module_definition.cmake)
+
+set(HEADER_FILES
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/atlasmanager.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/brickmanager.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/brickselector.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/brickcover.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/brickselection.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/multiresvolumeraycaster.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/shenbrickselector.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/tfbrickselector.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/localtfbrickselector.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/simpletfbrickselector.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/renderablemultiresvolume.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/tsp.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/histogram.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/histogrammanager.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/errorhistogrammanager.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/localerrorhistogrammanager.h
+)
+source_group("Header Files" FILES ${HEADER_FILES})
+
+set(SOURCE_FILES
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/atlasmanager.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/brickmanager.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/multiresvolumeraycaster.cpp    
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/shenbrickselector.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/tfbrickselector.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/localtfbrickselector.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/simpletfbrickselector.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/renderablemultiresvolume.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/tsp.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/histogram.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/histogrammanager.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/errorhistogrammanager.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/rendering/localerrorhistogrammanager.cpp
+)
+source_group("Source Files" FILES ${SOURCE_FILES})
+
+create_new_module(
+    "MultiresVolume"
+    multiresvolume_module
+    ${HEADER_FILES} ${SOURCE_FILES}
+)

modules/multiresvolume/include.cmake

@@ -0,0 +1,4 @@
+set (DEFAULT_MODULE ON)
+set (OPENSPACE_DEPENDENCIES
+    volume
+)

modules/multiresvolume/multiresvolumemodule.cpp

@@ -0,0 +1,45 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2015                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/multiresvolumemodule.h>
+
+#include <openspace/rendering/renderable.h>
+#include <openspace/util/factorymanager.h>
+
+#include <ghoul/misc/assert.h>
+
+#include <modules/multiresvolume/rendering/renderablemultiresvolume.h>
+
+namespace openspace {
+
+MultiresVolumeModule::MultiresVolumeModule() : OpenSpaceModule("MultiresVolume") {}
+
+void MultiresVolumeModule::internalInitialize() {
+    auto fRenderable = FactoryManager::ref().factory<Renderable>();
+    ghoul_assert(fRenderable, "No renderable factory existed");
+
+    fRenderable->registerClass<RenderableMultiresVolume>("RenderableMultiresVolume");
+}
+
+} // namespace openspace

modules/multiresvolume/multiresvolumemodule.h

@@ -0,0 +1,40 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2015                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __MULTIRESVOLUMEMODULE_H__
+#define __MULTIRESVOLUMEMODULE_H__
+
+#include <openspace/util/openspacemodule.h>
+
+namespace openspace {
+
+class MultiresVolumeModule : public OpenSpaceModule {
+public:
+    MultiresVolumeModule();
+    void internalInitialize() override;
+};
+
+} // namespace openspace
+
+#endif // __MULTIRESVOLUMEMODULE_H__

modules/multiresvolume/rendering/atlasmanager.cpp

@@ -0,0 +1,263 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/atlasmanager.h>
+
+#include <ghoul/filesystem/filesystem.h>
+#include <ghoul/logging/logmanager.h>
+#include <ghoul/opengl/texture.h>
+
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <cstring>
+
+namespace {
+    const std::string _loggerCat = "AtlasManager";
+}
+
+namespace openspace {
+
+AtlasManager::AtlasManager(TSP* tsp) : _tsp(tsp) {}
+
+AtlasManager::~AtlasManager() {}
+
+bool AtlasManager::initialize() {
+    TSP::Header header = _tsp->header();
+
+    _nBricksPerDim = header.xNumBricks_;
+    _nOtLeaves = _nBricksPerDim * _nBricksPerDim * _nBricksPerDim;
+    _nOtNodes = _tsp->numOTNodes();
+    _nOtLevels = log(_nOtLeaves)/log(8) + 1;
+    _paddedBrickDim = _tsp->paddedBrickDim();
+    _nBricksInMap = _nBricksPerDim * _nBricksPerDim * _nBricksPerDim;
+    _atlasDim = _nBricksPerDim * _paddedBrickDim;
+    _nBrickVals = _paddedBrickDim*_paddedBrickDim*_paddedBrickDim;
+    _brickSize = _nBrickVals * sizeof(float);
+    _volumeSize = _brickSize * _nOtLeaves;
+    _atlasMap = std::vector<unsigned int>(_nOtLeaves, NOT_USED);
+    _nBricksInAtlas = _nBricksInMap;
+
+    _freeAtlasCoords = std::vector<unsigned int>(_nBricksInAtlas, 0);
+
+    for (unsigned int i = 0; i < _nBricksInAtlas; i++) {
+        _freeAtlasCoords[i] = i;
+    }
+
+    _textureAtlas = new ghoul::opengl::Texture(
+        glm::size3_t(_atlasDim, _atlasDim, _atlasDim), 
+        ghoul::opengl::Texture::Format::RGBA, 
+        GL_RGBA, 
+        GL_FLOAT);
+    _textureAtlas->uploadTexture();
+
+    glGenBuffers(2, _pboHandle);
+
+    glGenBuffers(1, &_atlasMapBuffer);
+    glBindBuffer(GL_SHADER_STORAGE_BUFFER, _atlasMapBuffer);
+    glBufferData(GL_SHADER_STORAGE_BUFFER, sizeof(GLint)*_nBricksInMap, NULL, GL_DYNAMIC_READ);
+    glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
+
+    return true;
+}
+
+std::vector<unsigned int> AtlasManager::atlasMap() {
+    return _atlasMap;
+}
+
+unsigned int AtlasManager::atlasMapBuffer() {
+    return _atlasMapBuffer;
+}
+
+void AtlasManager::updateAtlas(BUFFER_INDEX bufferIndex, std::vector<int>& brickIndices) {
+    int nBrickIndices = brickIndices.size();
+
+    _requiredBricks.clear();
+    for (int i = 0; i < nBrickIndices; i++) {
+        _requiredBricks.insert(brickIndices[i]);
+    }
+
+    for (unsigned int it : _prevRequiredBricks) {
+        if (!_requiredBricks.count(it)) {
+            removeFromAtlas(it);
+        }
+    }
+
+    // Stats
+    _nUsedBricks = _requiredBricks.size();
+    _nStreamedBricks = 0;
+    _nDiskReads = 0;
+
+    glBindBuffer(GL_PIXEL_UNPACK_BUFFER, _pboHandle[bufferIndex]);
+    glBufferData(GL_PIXEL_UNPACK_BUFFER, _volumeSize, 0, GL_STREAM_DRAW);
+    float* mappedBuffer = reinterpret_cast<float*>(glMapBuffer(GL_PIXEL_UNPACK_BUFFER, GL_WRITE_ONLY));
+
+    if (!mappedBuffer) {
+        LERROR("Failed to map PBO");
+        std::cout << glGetError() << std::endl;
+        return;
+    }
+
+    for (auto itStart = _requiredBricks.begin(); itStart != _requiredBricks.end();) {
+        int firstBrick = *itStart;
+        int lastBrick = firstBrick;
+
+        auto itEnd = itStart;
+        for (itEnd++; itEnd != _requiredBricks.end() && *itEnd == lastBrick + 1; itEnd++) {
+            lastBrick = *itEnd;
+        }
+
+        addToAtlas(firstBrick, lastBrick, mappedBuffer);
+
+        itStart = itEnd;
+    }
+
+    glUnmapBuffer(GL_PIXEL_UNPACK_BUFFER);
+    glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0);
+
+    for (int i = 0; i < nBrickIndices; i++) {
+        _atlasMap[i] = _brickMap[brickIndices[i]];
+    }
+
+    std::swap(_prevRequiredBricks, _requiredBricks);
+
+    pboToAtlas(bufferIndex);
+
+    glBindBuffer(GL_SHADER_STORAGE_BUFFER, _atlasMapBuffer);
+    GLint *to = (GLint*)glMapBuffer(GL_SHADER_STORAGE_BUFFER, GL_WRITE_ONLY);
+    memcpy(to, _atlasMap.data(), sizeof(GLint)*_atlasMap.size());
+    glUnmapBuffer(GL_SHADER_STORAGE_BUFFER);
+    glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
+
+}
+
+void AtlasManager::addToAtlas(int firstBrickIndex, int lastBrickIndex, float* mappedBuffer) {
+    while (_brickMap.count(firstBrickIndex) && firstBrickIndex <= lastBrickIndex) firstBrickIndex++;
+    while (_brickMap.count(lastBrickIndex) && lastBrickIndex >= firstBrickIndex) lastBrickIndex--;
+    if (lastBrickIndex < firstBrickIndex) return;
+    
+    int sequenceLength = lastBrickIndex - firstBrickIndex + 1;
+    float* sequenceBuffer = new float[sequenceLength*_nBrickVals];
+    size_t bufferSize = sequenceLength * _brickSize;
+
+    long long offset = TSP::dataPosition() + static_cast<long long>(firstBrickIndex) * static_cast<long long>(_brickSize);
+    _tsp->file().seekg(offset);
+    _tsp->file().read(reinterpret_cast<char*>(sequenceBuffer), bufferSize);
+    _nDiskReads++;
+
+    for (int brickIndex = firstBrickIndex; brickIndex <= lastBrickIndex; brickIndex++) {
+        if (!_brickMap.count(brickIndex)) {
+            unsigned int atlasCoords = _freeAtlasCoords.back();
+            _freeAtlasCoords.pop_back();
+            int level = _nOtLevels - floor(log((7.0 * (float(brickIndex % _nOtNodes)) + 1.0))/log(8)) - 1;
+            assert(atlasCoords <= 0x0FFFFFFF);
+            unsigned int atlasData = (level << 28) + atlasCoords;
+            _brickMap.insert(std::pair<unsigned int, unsigned int>(brickIndex, atlasData));
+            _nStreamedBricks++;
+            fillVolume(&sequenceBuffer[_nBrickVals*(brickIndex - firstBrickIndex)], mappedBuffer, atlasCoords);
+        }
+    }
+
+    delete[] sequenceBuffer;
+}
+
+void AtlasManager::removeFromAtlas(int brickIndex) {
+    unsigned int atlasData = _brickMap[brickIndex];
+    unsigned int atlasCoords = atlasData & 0x0FFFFFFF;
+    _brickMap.erase(brickIndex);
+    _freeAtlasCoords.push_back(atlasCoords);
+}
+
+void AtlasManager::fillVolume(float* in, float* out, unsigned int linearAtlasCoords) {
+    int x = linearAtlasCoords % _nBricksPerDim;
+    int y = (linearAtlasCoords / _nBricksPerDim) % _nBricksPerDim;
+    int z = linearAtlasCoords / _nBricksPerDim / _nBricksPerDim;
+
+
+    unsigned int xMin = x*_paddedBrickDim;
+    unsigned int yMin = y*_paddedBrickDim;
+    unsigned int zMin = z*_paddedBrickDim;
+    unsigned int xMax = xMin + _paddedBrickDim;
+    unsigned int yMax = yMin + _paddedBrickDim;
+    unsigned int zMax = zMin + _paddedBrickDim;
+
+    unsigned int from = 0;
+    for (unsigned int zValCoord = zMin; zValCoord<zMax; ++zValCoord) {
+        for (unsigned int yValCoord = yMin; yValCoord<yMax; ++yValCoord) {
+            for (unsigned int xValCoord = xMin; xValCoord<xMax; ++xValCoord) {
+
+                unsigned int idx =
+                    xValCoord +
+                    yValCoord*_atlasDim +
+                    zValCoord*_atlasDim*_atlasDim;
+
+                out[idx] = in[from];
+                from++;
+            }
+        }
+    }
+}
+
+void AtlasManager::pboToAtlas(BUFFER_INDEX bufferIndex) {
+    glBindBuffer(GL_PIXEL_UNPACK_BUFFER, _pboHandle[bufferIndex]);
+    glm::size3_t dim = _textureAtlas->dimensions();
+    glBindTexture(GL_TEXTURE_3D, *_textureAtlas);
+    glTexSubImage3D(GL_TEXTURE_3D,  // target
+        0,                          // level
+        0,                          // xoffset
+        0,                          // yoffset
+        0,                          // zoffset
+        dim[0],                     // width
+        dim[1],                     // height
+        dim[2],                     // depth
+        GL_RED,                     // format
+        GL_FLOAT,                   // type
+        NULL);                      // *pixels
+    glBindTexture(GL_TEXTURE_3D, 0);
+    glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0);
+}
+
+ghoul::opengl::Texture& AtlasManager::textureAtlas() {
+    ghoul_assert(_textureAtlas != nullptr, "Texture atlas is nullptr");
+    return *_textureAtlas;
+}
+
+unsigned int AtlasManager::getNumDiskReads() {
+    return _nDiskReads;
+}
+
+unsigned int AtlasManager::getNumUsedBricks() {
+    return _nUsedBricks;
+}
+
+unsigned int AtlasManager::getNumStreamedBricks() {
+    return _nStreamedBricks;
+}
+
+
+glm::size3_t AtlasManager::textureSize() {
+    return _textureAtlas->dimensions();
+}
+
+}

modules/multiresvolume/rendering/atlasmanager.h

@@ -0,0 +1,103 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __ATLASMANAGER_H__
+#define __ATLASMANAGER_H__
+
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <ghoul/glm.h>
+#include <glm/gtx/std_based_type.hpp>
+
+#include <string>
+#include <vector>
+#include <climits>
+#include <map>
+#include <set>
+
+namespace ghoul {
+    namespace opengl {
+        class Texture;
+    }
+}
+
+namespace openspace {
+
+class AtlasManager {
+public:
+    enum BUFFER_INDEX { EVEN = 0, ODD = 1 };
+
+    AtlasManager(TSP* tsp);
+    ~AtlasManager();
+
+    void updateAtlas(BUFFER_INDEX bufferIndex, std::vector<int>& brickIndices);
+    void addToAtlas(int firstBrickIndex, int lastBrickIndex, float* mappedBuffer);
+    void removeFromAtlas(int brickIndex);
+    bool initialize();
+    std::vector<unsigned int> atlasMap();
+    unsigned int atlasMapBuffer();
+
+    void pboToAtlas(BUFFER_INDEX bufferIndex);
+    ghoul::opengl::Texture& textureAtlas();
+    glm::size3_t textureSize();
+
+    unsigned int getNumDiskReads();
+    unsigned int getNumUsedBricks();
+    unsigned int getNumStreamedBricks();
+private:
+    const unsigned int NOT_USED = UINT_MAX;
+    TSP* _tsp;
+    unsigned int _pboHandle[2];
+    unsigned int _atlasMapBuffer;
+
+    std::vector<unsigned int> _atlasMap;
+    std::map<unsigned int, unsigned int> _brickMap;
+    std::vector<unsigned int> _freeAtlasCoords;
+    std::set<unsigned int> _requiredBricks;
+    std::set<unsigned int> _prevRequiredBricks;
+
+    ghoul::opengl::Texture* _textureAtlas;
+
+    // Stats
+    unsigned int _nUsedBricks;
+    unsigned int _nStreamedBricks;
+    unsigned int _nDiskReads;
+
+    unsigned int _nBricksPerDim,
+                 _nOtLeaves,
+                 _nOtNodes,
+                 _nOtLevels,
+                 _brickSize,
+                 _nBrickVals,
+                 _volumeSize,
+                 _paddedBrickDim,
+                 _nBricksInAtlas,
+                 _nBricksInMap,
+                 _atlasDim;
+
+    void fillVolume(float* in, float* out, unsigned int linearAtlasCoords);
+};
+
+} // namespace openspace
+
+#endif

modules/multiresvolume/rendering/brickcover.h

@@ -0,0 +1,68 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __BRICKCOVER_H__
+#define __BRICKCOVER_H__
+
+namespace openspace {
+
+struct BrickCover {
+    int lowX, highX, lowY, highY, lowZ, highZ;
+
+    BrickCover() {}
+
+    BrickCover(int numBricks) {
+        lowX = lowY = lowZ = 0;
+        highX = highY = highZ = numBricks;
+    }
+
+    BrickCover split(bool x, bool y, bool z) {
+        BrickCover child;
+        if (x) {
+            child.lowX = lowX + (highX - lowX) / 2;
+            child.highX = highX;
+        } else {
+            child.lowX = lowX;
+            child.highX = lowX + (highX - lowX) / 2;
+        }
+        if (y) {
+            child.lowY = lowY + (highY - lowY) / 2;
+            child.highY = highY;
+        } else {
+            child.lowY = lowY;
+            child.highY = lowY + (highY - lowY) / 2;
+        }
+        if (z) {
+            child.lowZ = lowZ + (highZ - lowZ) / 2;
+            child.highZ = highZ;
+        } else {
+            child.lowZ = lowZ;
+            child.highZ = lowZ + (highZ - lowZ) / 2;
+        }
+        return child;
+    }
+};
+}
+
+#endif // __BRICKCOVER_H__

modules/multiresvolume/rendering/brickmanager.cpp

@@ -0,0 +1,476 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/brickmanager.h>
+
+#include <ghoul/filesystem/filesystem.h>
+#include <ghoul/logging/logmanager.h>
+#include <ghoul/opengl/texture.h>
+
+#include <iostream>
+#include <fstream>
+
+namespace {
+	const std::string _loggerCat = "BrickManager";
+}
+
+namespace openspace {
+
+BrickManager::BrickManager(TSP* tsp) 
+	: _tsp(tsp)
+	, numBricks_(0)
+	, brickDim_(0)
+	, paddedBrickDim_(0)
+	, atlasDim_(0)
+	, numBrickVals_(0)
+	, numBricksFrame_(0)
+	, numBricksTree_(0)
+	, brickSize_(0)
+	, volumeSize_(0)
+	, numValsTot_(0)
+	, xCoord_(0)
+	, yCoord_(0)
+	, zCoord_(0)
+	, textureAtlas_(nullptr)
+	, hasReadHeader_(false)
+	, atlasInitialized_(false)
+{}
+
+BrickManager::~BrickManager() {
+
+}
+
+bool BrickManager::readHeader() {
+
+	if (!_tsp->file().is_open())
+		return false;
+
+	_header = _tsp->header();
+
+	LDEBUG("Grid type: " << _header.gridType_);
+	LDEBUG("Original num timesteps: " << _header.numOrigTimesteps_);
+	LDEBUG("Num timesteps: " << _header.numTimesteps_);
+	LDEBUG("Brick dims: " << _header.xBrickDim_ << " " << _header.yBrickDim_ << " " << _header.zBrickDim_);
+	LDEBUG("Num bricks: " << _header.xNumBricks_ << " " << _header.yNumBricks_ << " " << _header.zNumBricks_);
+	LDEBUG("");
+
+	brickDim_ = _header.xBrickDim_;
+	numBricks_ = _header.xNumBricks_;
+	paddedBrickDim_ = brickDim_ + paddingWidth_ * 2;
+	atlasDim_ = paddedBrickDim_*numBricks_;
+
+	LDEBUG("Padded brick dim: " << paddedBrickDim_);
+	LDEBUG("Atlas dim: " << atlasDim_);
+
+	numBrickVals_ = paddedBrickDim_*paddedBrickDim_*paddedBrickDim_;
+	// Number of bricks per frame
+	numBricksFrame_ = numBricks_*numBricks_*numBricks_;
+
+	// Calculate number of bricks in tree
+	unsigned int numOTLevels = static_cast<unsigned int>(log((int)numBricks_) / log(2) + 1);
+	unsigned int numOTNodes = static_cast<unsigned int>((pow(8, numOTLevels) - 1) / 7);
+	unsigned int numBSTNodes = static_cast<unsigned int>(_header.numTimesteps_ * 2 - 1);
+	numBricksTree_ = numOTNodes * numBSTNodes;
+	LDEBUG("Num OT levels: " << numOTLevels);
+	LDEBUG("Num OT nodes: " << numOTNodes);
+	LDEBUG("Num BST nodes: " << numBSTNodes);
+	LDEBUG("Num bricks in tree: " << numBricksTree_);
+	LDEBUG("Num values per brick: " << numBrickVals_);
+
+	brickSize_ = sizeof(float)*numBrickVals_;
+	volumeSize_ = brickSize_*numBricksFrame_;
+	numValsTot_ = numBrickVals_*numBricksFrame_;
+
+	_tsp->file().seekg(0, _tsp->file().end);
+	long long fileSize = _tsp->file().tellg();
+	long long calcFileSize = static_cast<long long>(numBricksTree_)*
+		static_cast<long long>(brickSize_) + TSP::dataPosition();
+
+
+	if (fileSize != calcFileSize) {
+		LERROR("Sizes don't match");
+		LERROR("calculated file size: " << calcFileSize);
+		LERROR("file size: " << fileSize);
+		return false;
+	}
+
+	hasReadHeader_ = true;
+
+	// Hold two brick lists
+	brickLists_.resize(2);
+	// Make sure the brick list can hold the maximum number of bricks
+	// Each entry holds tree coordinates
+	brickLists_[EVEN].resize(numBricksTree_ * 3, -1);
+	brickLists_[ODD].resize(numBricksTree_ * 3, -1);
+
+	// Allocate space for keeping tracks of bricks in PBO
+	bricksInPBO_.resize(2);
+	bricksInPBO_[EVEN].resize(numBricksTree_, -1);
+	bricksInPBO_[ODD].resize(numBricksTree_, -1);
+
+	// Allocate space for keeping track of the used coordinates in atlas
+	usedCoords_.resize(2);
+	usedCoords_[EVEN].resize(numBricksFrame_, false);
+	usedCoords_[ODD].resize(numBricksFrame_, false);
+
+	return true;
+}
+
+bool BrickManager::initialize() {
+	if (atlasInitialized_) {
+		LWARNING("InitAtlas() - already initialized");
+	}
+
+	if (!hasReadHeader_) {
+		LWARNING("InitAtlas() - Has not read header, trying to read");
+		return readHeader();
+	}
+
+	// Prepare the 3D texture
+	std::vector<unsigned int> dims;
+	dims.push_back(atlasDim_);
+	dims.push_back(atlasDim_);
+	dims.push_back(atlasDim_);
+	textureAtlas_ = new ghoul::opengl::Texture(
+		glm::size3_t(atlasDim_, atlasDim_, atlasDim_), 
+		ghoul::opengl::Texture::Format::RGBA, 
+		GL_RGBA, 
+		GL_FLOAT);
+	textureAtlas_->uploadTexture();
+	//textureAtlas_ = Texture3D::New(dims);
+
+	//if (!textureAtlas_->Init()) return false;
+
+	atlasInitialized_ = true;
+
+
+	glGenBuffers(2, pboHandle_);
+
+	return true;
+}
+
+bool BrickManager::BuildBrickList(BUFFER_INDEX _bufIdx,
+	std::vector<int> &_brickRequest) {
+
+	// Keep track of number bricks used and number of bricks cached
+	// (for benchmarking)
+	int numBricks = 0;
+	int numCached = 0;
+
+	// For every non-zero entry in the request list, assign a texture atlas
+	// coordinate. For zero entries, signal "no brick" using -1.
+
+	for (unsigned int i = 0; i<_brickRequest.size(); ++i) {
+
+		if (_brickRequest[i] > 0) {
+
+			numBricks++;
+
+			//INFO("Checking brick " << i);
+
+			// If the brick is already in the atlas, keep the coordinate
+			if (bricksInPBO_[_bufIdx][i] != -1) {
+
+				numCached++;
+
+				// Get the corresponding coordinates from index
+				int x, y, z;
+				CoordsFromLin(bricksInPBO_[_bufIdx][i], x, y, z);
+				brickLists_[_bufIdx][3 * i + 0] = x;
+				brickLists_[_bufIdx][3 * i + 1] = y;
+				brickLists_[_bufIdx][3 * i + 2] = z;
+
+				// Mark coordinate as used
+				usedCoords_[_bufIdx][bricksInPBO_[_bufIdx][i]] = true;
+
+			}
+			else {
+
+				// If coord is already usedi by another brick, 
+				// skip it and try the next one
+				while (usedCoords_[_bufIdx][LinearCoord(xCoord_, yCoord_, zCoord_)]) {
+					IncCoord();
+				}
+
+				brickLists_[_bufIdx][3 * i + 0] = xCoord_;
+				brickLists_[_bufIdx][3 * i + 1] = yCoord_;
+				brickLists_[_bufIdx][3 * i + 2] = zCoord_;
+				usedCoords_[_bufIdx][LinearCoord(xCoord_, yCoord_, zCoord_)] = true;
+
+				IncCoord();
+			}
+
+
+		}
+		else {
+
+			// -1 is for "not used"
+			brickLists_[_bufIdx][3 * i + 0] = -1;
+			brickLists_[_bufIdx][3 * i + 1] = -1;
+			brickLists_[_bufIdx][3 * i + 2] = -1;
+
+		}
+
+		// Reset brick list during iteration
+		_brickRequest[i] = 0;
+
+	}
+
+
+	// Brick list is build, reset coordinate list
+	for (auto it = usedCoords_[_bufIdx].begin();
+		it != usedCoords_[_bufIdx].end(); ++it) {
+		*it = false;
+	}
+
+	//INFO("bricks NOT used: " << (float)(numBricksFrame_-numBricks) / (float)(numBricksFrame_));
+	//INFO("bricks cached: " << (float)numCached / (float)(numBricksFrame_));
+
+	return true;
+}
+
+bool BrickManager::FillVolume(float *_in, float *_out,
+	unsigned int _x,
+	unsigned int _y,
+	unsigned int _z) {
+
+	//timer_.start();
+	unsigned int xMin = _x*paddedBrickDim_;
+	unsigned int yMin = _y*paddedBrickDim_;
+	unsigned int zMin = _z*paddedBrickDim_;
+	unsigned int xMax = xMin + paddedBrickDim_;
+	unsigned int yMax = yMin + paddedBrickDim_;
+	unsigned int zMax = zMin + paddedBrickDim_;
+
+	// Loop over the brick using three loops
+	unsigned int from = 0;
+	for (unsigned int zValCoord = zMin; zValCoord<zMax; ++zValCoord) {
+		for (unsigned int yValCoord = yMin; yValCoord<yMax; ++yValCoord) {
+			for (unsigned int xValCoord = xMin; xValCoord<xMax; ++xValCoord) {
+
+				unsigned int idx =
+					xValCoord +
+					yValCoord*atlasDim_ +
+					zValCoord*atlasDim_*atlasDim_;
+
+				_out[idx] = _in[from];
+				from++;
+			}
+		}
+	}
+	return true;
+}
+
+void BrickManager::IncCoord() {
+	// Update atlas coordinate
+	xCoord_++;
+	if (xCoord_ == _header.xNumBricks_) {
+		xCoord_ = 0;
+		yCoord_++;
+		if (yCoord_ == _header.yNumBricks_) {
+			yCoord_ = 0;
+			zCoord_++;
+			if (zCoord_ == _header.zNumBricks_) {
+				zCoord_ = 0;
+			}
+		}
+	}
+}
+
+unsigned int BrickManager::LinearCoord(int _x, int _y, int _z) {
+	return _x + _y*_header.xNumBricks_ + _z*_header.xNumBricks_*_header.yNumBricks_;
+}
+
+void BrickManager::CoordsFromLin(int _idx, int &_x, int &_y, int &_z) {
+	_x = _idx % _header.xNumBricks_;
+	_idx /= _header.xNumBricks_;
+	_y = _idx % _header.yNumBricks_;
+	_idx /= _header.yNumBricks_;
+	_z = _idx;
+}
+
+
+bool BrickManager::DiskToPBO(BUFFER_INDEX _pboIndex) {
+
+	// Map PBO
+	glBindBuffer(GL_PIXEL_UNPACK_BUFFER, pboHandle_[_pboIndex]);
+	glBufferData(GL_PIXEL_UNPACK_BUFFER, volumeSize_, 0, GL_STREAM_DRAW);
+	float *mappedBuffer = reinterpret_cast<float*>(
+		glMapBuffer(GL_PIXEL_UNPACK_BUFFER, GL_WRITE_ONLY));
+
+	if (!mappedBuffer) {
+		LERROR("Failed to map PBO");
+		return false;
+	}
+
+	// Loop over brick request list
+	unsigned int brickIndex = 0;
+	while (brickIndex < brickLists_[_pboIndex].size() / 3) {
+
+		// Find first active brick index in list
+		while (brickIndex<brickLists_[_pboIndex].size() / 3 &&
+			brickLists_[_pboIndex][3 * brickIndex] == -1) {
+			// If not used, remove from PBO cache list
+			bricksInPBO_[_pboIndex][brickIndex] = -1;
+			brickIndex++;
+		}
+
+		// If we are at the end of the list, exit
+		if (brickIndex == brickLists_[_pboIndex].size() / 3) {
+			break;
+		}
+
+		// Find a sequence of consecutive bricks in list
+		unsigned int sequence = 0;
+		// Count number of bricks already in PBO
+		unsigned int inPBO = 0;
+		unsigned int brickIndexProbe = brickIndex;
+		while (brickIndexProbe < brickLists_[_pboIndex].size() / 3 &&
+			brickLists_[_pboIndex][3 * brickIndexProbe] != -1) {
+			sequence++;
+			if (bricksInPBO_[_pboIndex][brickIndexProbe] != -1) {
+				inPBO++;
+			}
+			brickIndexProbe++;
+		}
+		//INFO("Reading " << sequence << " bricks");
+
+		// Read the sequence into a buffer
+		float *seqBuffer = new float[sequence*numBrickVals_];
+		size_t bufSize = sequence*numBrickVals_*sizeof(float);
+		/*
+		std::ios::pos_type offset = dataPos_ +
+		static_cast<std::ios::pos_type>(brickIndex) *
+		static_cast<std::ios::pos_type>(brickSize_);
+		*/
+
+		long offset = TSP::dataPosition() +
+			static_cast<long>(brickIndex)*
+			static_cast<long>(brickSize_);
+
+		// Skip reading if all bricks in sequence is already in PBO
+		if (inPBO != sequence) {
+
+			//timer_.start();
+
+			/*
+			std::streamoff off = static_cast<std::streamoff>(offset);
+			in_.seekg(off);
+			if (in_.tellg() == -1) {
+			ERROR("Failed to get input stream position");
+			INFO("offset: " << offset);
+			INFO("streamoff max: " << std::numeric_limits<std::streamoff>::max());
+			INFO("size_t max: " << std::numeric_limits<size_t>::max());
+			return false;
+			}
+			INFO("in.tellg(): " << in_.tellg());
+			in_.read(reinterpret_cast<char*>(seqBuffer), brickSize_*sequence);
+			*/
+
+			_tsp->file().seekg(offset);
+			_tsp->file().read(reinterpret_cast<char*>(seqBuffer), bufSize);
+
+
+			//timer_.stop();
+			//double time = timer_.elapsed().wall / 1.0e9;
+			//double mb = (brickSize_*sequence) / 1048576.0;
+			//INFO("Disk read "<<mb<<" MB in "<<time<<" s, "<< mb/time<<" MB/s");
+
+			// For each brick in the buffer, put it the correct buffer spot
+			for (unsigned int i = 0; i<sequence; ++i) {
+
+
+				// Only upload if needed
+				// Pointless if implementation only skips reading when ALL bricks in
+				// sequence are in PBO, but could be useful if other solutions that
+				// considers part of the buffer are implemented
+				if (bricksInPBO_[_pboIndex][brickIndex + i] == -1) {
+
+					unsigned int x = static_cast<unsigned int>(
+						brickLists_[_pboIndex][3 * (brickIndex + i) + 0]);
+					unsigned int y = static_cast<unsigned int>(
+						brickLists_[_pboIndex][3 * (brickIndex + i) + 1]);
+					unsigned int z = static_cast<unsigned int>(
+						brickLists_[_pboIndex][3 * (brickIndex + i) + 2]);
+
+					// Put each brick in the correct buffer place.
+					// This needs to be done because the values are in brick order, and
+					// the volume needs to be filled with one big float array.
+					FillVolume(&seqBuffer[numBrickVals_*i], mappedBuffer, x, y, z);
+					// Update the atlas list since the brick will be uploaded
+					//INFO(brickIndex+i);
+					bricksInPBO_[_pboIndex][brickIndex + i] = LinearCoord(x, y, z);
+
+				}
+			}
+
+			delete[] seqBuffer;
+
+		} // if in pbo
+
+		// Update the brick index
+		brickIndex += sequence;
+
+	}
+
+	glUnmapBuffer(GL_PIXEL_UNPACK_BUFFER);
+	glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0);
+
+	return true;
+}
+
+bool BrickManager::PBOToAtlas(BUFFER_INDEX _pboIndex) {
+	glBindBuffer(GL_PIXEL_UNPACK_BUFFER, pboHandle_[_pboIndex]);
+	glm::size3_t dim = textureAtlas_->dimensions();
+	glBindTexture(GL_TEXTURE_3D, *textureAtlas_);
+	glTexSubImage3D(GL_TEXTURE_3D,	// target
+		0,							// level
+		0,							// xoffset
+		0,							// yoffset
+		0,							// zoffset
+		dim[0],						// width
+		dim[1],						// height
+		dim[2],						// depth
+		GL_RED,						// format
+		GL_FLOAT,					// type
+		NULL);						// *pixels
+	glBindTexture(GL_TEXTURE_3D, 0);
+	glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0);
+
+	return true;
+}
+
+ghoul::opengl::Texture* BrickManager::textureAtlas() {
+	return textureAtlas_;
+}
+
+unsigned int BrickManager::pbo(BUFFER_INDEX _pboIndex) {
+	return pboHandle_[_pboIndex];
+}
+
+const std::vector<int>& BrickManager::brickList(BUFFER_INDEX index) const {
+	return brickLists_.at(index);
+}
+
+}

modules/multiresvolume/rendering/brickmanager.h

@@ -0,0 +1,111 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __BRICKMANAGER_H__
+#define __BRICKMANAGER_H__
+
+#include <modules/multiresvolume/rendering/tsp.h>
+
+#include <string>
+#include <vector>
+
+namespace ghoul {
+	namespace opengl {
+		class Texture;
+	}
+}
+
+namespace openspace {
+
+class BrickManager {
+public:
+	enum BUFFER_INDEX { EVEN = 0, ODD = 1 };
+
+	BrickManager(TSP* tsp);
+	~BrickManager();
+
+	bool readHeader();
+
+	bool initialize(); 
+
+	bool BuildBrickList(BUFFER_INDEX _bufIdx, std::vector<int> &_brickRequest);
+
+	bool FillVolume(float *_in, float *_out,
+		unsigned int _x,
+		unsigned int _y,
+		unsigned int _z);
+	bool DiskToPBO(BUFFER_INDEX _pboIndex);
+	bool PBOToAtlas(BUFFER_INDEX _pboIndex);
+
+	ghoul::opengl::Texture* textureAtlas();
+	unsigned int pbo(BUFFER_INDEX _pboIndex);
+	const std::vector<int>& brickList(BUFFER_INDEX index) const;
+
+private:
+
+	void IncCoord();
+	unsigned int LinearCoord(int _x, int _y, int _z);
+	void CoordsFromLin(int _idx, int &_x, int &_y, int &_z);
+
+	TSP* _tsp;
+	TSP::Header _header;
+
+	unsigned int numBricks_;
+	unsigned int brickDim_;
+	unsigned int paddedBrickDim_;
+	unsigned int atlasDim_;
+
+	const unsigned int paddingWidth_ = 1;
+
+	unsigned int numBrickVals_;
+	unsigned int numBricksFrame_;
+	unsigned int numBricksTree_;
+	unsigned int brickSize_;
+	unsigned int volumeSize_;
+	unsigned int numValsTot_;
+
+	// Texture coordinates to be assigned
+	int xCoord_;
+	int yCoord_;
+	int zCoord_;
+
+	// Texture where the actual atlas is kept
+	ghoul::opengl::Texture* textureAtlas_;
+
+	std::vector<std::vector<int> > brickLists_;
+
+	bool hasReadHeader_;
+	bool atlasInitialized_;
+
+	// PBOs
+	unsigned int pboHandle_[2];
+
+	// Caching, one for each PBO
+	std::vector<std::vector<int> > bricksInPBO_;
+	std::vector<std::vector<bool> > usedCoords_;
+};
+
+} // namespace openspace
+
+#endif

modules/multiresvolume/rendering/brickselection.h

@@ -0,0 +1,108 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __BRICKSELECTION_H__
+#define __BRICKSELECTION_H__
+
+#include <modules/multiresvolume/rendering/brickcover.h>
+
+
+namespace openspace {
+
+struct BrickSelection {
+    enum SplitType {
+        None = 0,
+        Temporal = 1,
+        Spatial = 2
+    };
+
+    unsigned int brickIndex;
+    float splitPoints;
+    BrickSelection::SplitType splitType;
+    BrickCover cover;
+    int lowT;
+    int highT;
+    int nSpatialSplits;
+    int nTemporalSplits;
+
+    BrickSelection() {}
+
+    BrickSelection(int numBricks, int numTimeSteps, SplitType splitType, float splitPoints) {
+        this->cover = BrickCover(numBricks);
+        this->lowT = 0;
+        this->highT = numTimeSteps;
+        this->brickIndex = 0;
+        this->splitType = splitType;
+        this->splitPoints = splitPoints;
+        this->nSpatialSplits = 0;
+        this->nTemporalSplits = 0;
+    }
+
+    BrickSelection splitSpatially(bool x, bool y, bool z, unsigned int childBrickIndex, SplitType childSplitType, float childSplitPoints) {
+        BrickSelection child;
+        child.cover = cover.split(x, y, z);
+        child.brickIndex = childBrickIndex;
+        child.splitPoints = childSplitPoints;
+        child.splitType = childSplitType;
+        child.nSpatialSplits = nSpatialSplits + 1;
+        child.nTemporalSplits = nTemporalSplits;
+        child.lowT = lowT;
+        child.highT = highT;
+        return child;
+    }
+
+    BrickSelection splitTemporally(bool t, unsigned int childBrickIndex, SplitType childSplitType, float childSplitPoints) {
+        BrickSelection child;
+        child.cover = cover;
+        child.brickIndex = childBrickIndex;
+        child.splitPoints = childSplitPoints;
+        child.splitType = childSplitType;
+        if (t) {
+            child.lowT = centerT();
+            child.highT = highT;
+        } else {
+            child.lowT = lowT;
+            child.highT = centerT();
+        }
+        child.nSpatialSplits = nSpatialSplits;
+        child.nTemporalSplits = nTemporalSplits + 1;
+
+        return child;
+    }
+
+    int centerT() {
+        return lowT + (highT - lowT) / 2;
+    }
+
+    bool timestepInRightChild(int timestep) {
+        return timestep >= centerT();
+    }
+
+    static bool compareSplitPoints(const BrickSelection& a, const BrickSelection& b) {
+        return a.splitPoints < b.splitPoints;
+    }
+};
+}
+
+#endif // __BRICKSELECTION_H__

modules/multiresvolume/rendering/brickselector.h

@@ -0,0 +1,42 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __BRICKSELECTOR_H__
+#define __BRICKSELECTOR_H__
+
+#include <vector>
+
+namespace openspace {
+
+class BrickSelector {
+public:
+    virtual ~BrickSelector() {};
+	virtual bool initialize() { return true; };
+    virtual void selectBricks(int timestep,
+                      std::vector<int>& bricks) = 0;
+};
+
+} // namespace openspace
+
+#endif // __BRICKSELECTOR_H__

modules/multiresvolume/rendering/errorhistogrammanager.cpp

@@ -0,0 +1,361 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <float.h>
+#include <map>
+#include <string.h>
+#include <cmath>
+
+#include <modules/multiresvolume/rendering/errorhistogrammanager.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+
+#include <openspace/util/progressbar.h>
+
+#include <ghoul/logging/logmanager.h>
+
+namespace {
+    const std::string _loggerCat = "ErrorHistogramManager";
+}
+
+namespace openspace {
+
+ErrorHistogramManager::ErrorHistogramManager(TSP* tsp) : _tsp(tsp) {}
+
+ErrorHistogramManager::~ErrorHistogramManager() {}
+
+bool ErrorHistogramManager::buildHistograms(int numBins) {
+    _numBins = numBins;
+
+    _file = &(_tsp->file());
+    if (!_file->is_open()) {
+        return false;
+    }
+    _minBin = 0.0; // Should be calculated from tsp file
+    _maxBin = 1.0; // Should be calculated from tsp file as (maxValue - minValue)
+
+    unsigned int numOtLevels = _tsp->numOTLevels();
+    unsigned int numOtLeaves = pow(8, numOtLevels - 1);
+    unsigned int numBstLeaves = pow(2, _tsp->numBSTLevels() - 1);
+
+    _numInnerNodes = _tsp->numTotalNodes() - numOtLeaves * numBstLeaves;
+    _histograms = std::vector<Histogram>(_numInnerNodes);
+    LINFO("Build " << _numInnerNodes << " histograms with " << numBins << " bins each");
+
+    // All TSP Leaves
+    int numOtNodes = _tsp->numOTNodes();
+    int otOffset = (pow(8, numOtLevels - 1) - 1) / 7;
+
+    int numBstNodes = _tsp->numBSTNodes();
+    int bstOffset = numBstNodes / 2;
+
+    int numberOfLeaves = (numBstNodes - bstOffset) * (numOtNodes - otOffset);
+    ProgressBar pb(numberOfLeaves);
+    int processedLeaves = 0;
+    bool success = true;
+    for (int bst = bstOffset; bst < numBstNodes; bst++) {
+        for (int ot = otOffset; ot < numOtNodes; ot++) {
+            success &= buildFromLeaf(bst, ot);
+            if (!success) return false;
+            pb.print(++processedLeaves);
+        }
+    }
+    
+    return success;
+}
+
+bool ErrorHistogramManager::buildFromLeaf(unsigned int bstOffset, unsigned int octreeOffset) {
+    // Traverse all ancestors of leaf and add errors to their histograms
+
+    unsigned int brickDim = _tsp->brickDim();
+    unsigned int paddedBrickDim = _tsp->paddedBrickDim();
+    unsigned int padding = (paddedBrickDim - brickDim) / 2;
+
+    int numOtNodes = _tsp->numOTNodes();
+    unsigned int leafIndex = bstOffset * numOtNodes + octreeOffset;
+    std::vector<float> leafValues = readValues(leafIndex);
+    int numVoxels = leafValues.size();
+
+    int bstNode = bstOffset;
+    bool bstRightOnly = true;
+    unsigned int bstLevel = 0;
+
+    do {
+
+        glm::vec3 leafOffset(0.0); // Leaf offset in leaf sized voxels
+        unsigned int octreeLevel = 0;
+        unsigned int octreeNode = octreeOffset;
+        bool octreeLastOnly = true;
+        do {
+            // Visit ancestor
+            if (bstNode != bstOffset || octreeNode != octreeOffset) {
+                // Is actually an ancestor
+
+                std::vector<float> ancestorVoxels;
+                unsigned int ancestorBrickIndex = bstNode * numOtNodes + octreeNode;
+                unsigned int innerNodeIndex = brickToInnerNodeIndex(ancestorBrickIndex);
+                auto it = _voxelCache.find(innerNodeIndex);
+                if (it == _voxelCache.end()) {
+                    // First visit
+                    _histograms[innerNodeIndex] = Histogram(_minBin, _maxBin, _numBins);
+                    ancestorVoxels = readValues(ancestorBrickIndex);
+                    _voxelCache[innerNodeIndex] = ancestorVoxels;
+                } else {
+                    ancestorVoxels = it->second;
+                }
+
+                float voxelScale = pow(2, octreeLevel);
+                float invVoxelScale = 1.0 / voxelScale;
+
+                // Calculate leaf offset in ancestor sized voxels
+                glm::vec3 ancestorOffset = (leafOffset * invVoxelScale) + glm::vec3(padding - 0.5);
+
+                for (int z = 0; z < brickDim; z++) {
+                    for (int y = 0; y < brickDim; y++) {
+                        for (int x = 0; x < brickDim; x++) {
+                            glm::vec3 leafSamplePoint = glm::vec3(x, y, z) + glm::vec3(padding);
+                            glm::vec3 ancestorSamplePoint = ancestorOffset + (glm::vec3(x, y, z) + glm::vec3(0.5)) * invVoxelScale;
+                            float leafValue = leafValues[linearCoords(leafSamplePoint)];
+                            float ancestorValue = interpolate(ancestorSamplePoint, ancestorVoxels);
+
+                            _histograms[innerNodeIndex].addRectangle(leafValue, ancestorValue, std::abs(leafValue - ancestorValue));
+                        }
+                    }
+                }
+
+                if (bstRightOnly && octreeLastOnly) {
+                    _voxelCache.erase(innerNodeIndex);
+                }
+            }
+
+            // Traverse to next octree ancestor
+            int octreeChild = (octreeNode - 1) % 8;
+            octreeLastOnly &= octreeChild == 7;
+            octreeNode = parentOffset(octreeNode, 8);
+
+            int childSize = pow(2, octreeLevel) * brickDim;
+            leafOffset.x += (octreeChild % 2) * childSize;
+            leafOffset.y += ((octreeChild / 2) % 2) * childSize;
+            leafOffset.z += (octreeChild / 4) * childSize;
+
+            octreeLevel++;
+        } while (octreeNode != -1);
+
+        bstRightOnly &= (bstNode % 2 == 0);
+        bstNode = parentOffset(bstNode, 2);
+
+        bstLevel++;
+    } while (bstNode != -1);
+
+    return true;
+}
+
+bool ErrorHistogramManager::loadFromFile(const std::string& filename) {
+    std::ifstream file(filename, std::ios::in | std::ios::binary);
+    if (!file.is_open()) {
+        return false;
+    }
+
+    file.read(reinterpret_cast<char*>(&_numInnerNodes), sizeof(int));
+    file.read(reinterpret_cast<char*>(&_numBins), sizeof(int));
+    file.read(reinterpret_cast<char*>(&_minBin), sizeof(float));
+    file.read(reinterpret_cast<char*>(&_maxBin), sizeof(float));
+
+    int nFloats = _numInnerNodes * _numBins;
+    float* histogramData = new float[nFloats];
+    file.read(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+
+    _histograms = std::vector<Histogram>(_numInnerNodes);
+
+    for (int i = 0; i < _numInnerNodes; ++i) {
+        int offset = i*_numBins;
+        float* data = new float[_numBins];
+        memcpy(data, &histogramData[offset], sizeof(float) * _numBins);
+        _histograms[i] = Histogram(_minBin, _maxBin, _numBins, data);
+    }
+
+    delete[] histogramData;
+    // No need to deallocate histogram data, since histograms take ownership.
+    file.close();
+    return true;
+}
+
+
+bool ErrorHistogramManager::saveToFile(const std::string& filename) {
+    std::ofstream file(filename, std::ios::out | std::ios::binary);
+    if (!file.is_open()) {
+        return false;
+    }
+
+    file.write(reinterpret_cast<char*>(&_numInnerNodes), sizeof(int));
+    file.write(reinterpret_cast<char*>(&_numBins), sizeof(int));
+    file.write(reinterpret_cast<char*>(&_minBin), sizeof(float));
+    file.write(reinterpret_cast<char*>(&_maxBin), sizeof(float));
+
+    int nFloats = _numInnerNodes * _numBins;
+    float* histogramData = new float[nFloats];
+
+    for (int i = 0; i < _numInnerNodes; ++i) {
+        int offset = i*_numBins;
+        memcpy(&histogramData[offset], _histograms[i].data(), sizeof(float) * _numBins);
+    }
+
+    file.write(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+    delete[] histogramData;
+
+    file.close();
+    return true;
+}
+
+unsigned int ErrorHistogramManager::linearCoords(glm::vec3 coords) const {
+    return linearCoords(glm::ivec3(coords));
+}
+
+unsigned int ErrorHistogramManager::linearCoords(int x, int y, int z) const {
+    return linearCoords(glm::ivec3(x, y, z));
+}
+
+unsigned int ErrorHistogramManager::linearCoords(glm::ivec3 coords) const {
+    unsigned int paddedBrickDim = _tsp->paddedBrickDim();
+    return coords.z * paddedBrickDim * paddedBrickDim + coords.y * paddedBrickDim + coords.x;
+}
+
+float ErrorHistogramManager::interpolate(glm::vec3 samplePoint, const std::vector<float>& voxels) const {
+    int lowX = samplePoint.x;
+    int lowY = samplePoint.y;
+    int lowZ = samplePoint.z;
+
+    int highX = ceil(samplePoint.x);
+    int highY = ceil(samplePoint.y);
+    int highZ = ceil(samplePoint.z);
+
+    float interpolatorX = 1.0 - (samplePoint.x - lowX);
+    float interpolatorY = 1.0 - (samplePoint.y - lowY);
+    float interpolatorZ = 1.0 - (samplePoint.z - lowZ);
+
+    float v000 = voxels[linearCoords(lowX, lowY, lowZ)];
+    float v001 = voxels[linearCoords(lowX, lowY, highZ)];
+    float v010 = voxels[linearCoords(lowX, highY, lowZ)];
+    float v011 = voxels[linearCoords(lowX, highY, highZ)];
+    float v100 = voxels[linearCoords(highX, lowY, lowZ)];
+    float v101 = voxels[linearCoords(highX, lowY, highZ)];
+    float v110 = voxels[linearCoords(highX, highY, lowZ)];
+    float v111 = voxels[linearCoords(highX, highY, highZ)];
+
+    float v00 = interpolatorZ * v000 + (1.0 - interpolatorZ) * v001;
+    float v01 = interpolatorZ * v010 + (1.0 - interpolatorZ) * v011;
+    float v10 = interpolatorZ * v100 + (1.0 - interpolatorZ) * v101;
+    float v11 = interpolatorZ * v110 + (1.0 - interpolatorZ) * v111;
+
+    float v0 = interpolatorY * v00 + (1.0 - interpolatorY) * v01;
+    float v1 = interpolatorY * v10 + (1.0 - interpolatorY) * v11;
+
+    return interpolatorX * v0 + (1.0 - interpolatorX) * v1;
+
+}
+
+const Histogram* ErrorHistogramManager::getHistogram(unsigned int brickIndex) const {
+    unsigned int innerNodeIndex = brickToInnerNodeIndex(brickIndex);
+    if (innerNodeIndex < _numInnerNodes) {
+        return &(_histograms[innerNodeIndex]);
+    } else {
+        return nullptr;
+    }
+}
+
+int ErrorHistogramManager::parentOffset(int offset, int base) const {
+    if (offset == 0) {
+        return -1;
+    }
+    int depth = floor(log(((base - 1) * offset + 1.0)) / log(base));
+    int firstInLevel = (pow(base, depth) - 1) / (base - 1);
+    int inLevelOffset = offset - firstInLevel;
+
+    int parentDepth = depth - 1;
+    int firstInParentLevel = (pow(base, parentDepth) - 1) / (base - 1);
+    int parentInLevelOffset = inLevelOffset / base;
+
+    int parentOffset = firstInParentLevel + parentInLevelOffset;
+    return parentOffset;
+}
+
+std::vector<float> ErrorHistogramManager::readValues(unsigned int brickIndex) const {
+    unsigned int paddedBrickDim = _tsp->paddedBrickDim();
+    unsigned int numBrickVals = paddedBrickDim * paddedBrickDim * paddedBrickDim;
+    std::vector<float> voxelValues(numBrickVals);
+
+    std::streampos offset = _tsp->dataPosition() + static_cast<long long>(brickIndex*numBrickVals*sizeof(float));
+    _file->seekg(offset);
+
+    _file->read(reinterpret_cast<char*>(&voxelValues[0]),
+        static_cast<size_t>(numBrickVals)*sizeof(float));
+
+    return voxelValues;
+}
+
+unsigned int ErrorHistogramManager::brickToInnerNodeIndex(unsigned int brickIndex) const {
+    unsigned int numOtNodes = _tsp->numOTNodes();
+    unsigned int numBstLevels = _tsp->numBSTLevels();
+
+    unsigned int numInnerBstNodes = (pow(2, numBstLevels - 1) - 1) * numOtNodes;
+    if (brickIndex < numInnerBstNodes) return brickIndex;
+
+    unsigned int numOtLeaves = pow(8, _tsp->numOTLevels() - 1);
+    unsigned int numOtInnerNodes = (numOtNodes - numOtLeaves);
+
+    unsigned int innerBstOffset = brickIndex - numInnerBstNodes;
+    unsigned int rowIndex = innerBstOffset / numOtNodes;
+    unsigned int indexInRow = innerBstOffset % numOtNodes;
+
+    if (indexInRow >= numOtInnerNodes) return -1;
+
+    unsigned int offset = rowIndex * numOtInnerNodes;
+    unsigned int leavesOffset = offset + indexInRow;
+
+    return numInnerBstNodes + leavesOffset;
+}
+
+unsigned int ErrorHistogramManager::innerNodeToBrickIndex(unsigned int innerNodeIndex) const {
+    if (innerNodeIndex < 0 || innerNodeIndex >= _numInnerNodes) return -1; // Not an inner node
+
+    unsigned int numOtNodes = _tsp->numOTNodes();
+    unsigned int numBstLevels = _tsp->numBSTLevels();
+
+    unsigned int numInnerBstNodes = (pow(2, numBstLevels - 1) - 1) * numOtNodes;
+    if (innerNodeIndex < numInnerBstNodes) return innerNodeIndex;
+
+    unsigned int numOtLeaves = pow(8, _tsp->numOTLevels() - 1);
+    unsigned int numOtInnerNodes = (numOtNodes - numOtLeaves);
+
+    unsigned int innerBstOffset = innerNodeIndex - numInnerBstNodes;
+    unsigned int rowIndex = innerBstOffset / numOtInnerNodes;
+    unsigned int indexInRow = innerBstOffset % numOtInnerNodes;
+
+    unsigned int offset = rowIndex * numOtNodes;
+    unsigned int leavesOffset = offset + indexInRow;
+
+    return numInnerBstNodes + leavesOffset;
+}
+
+} // namespace openspace
+

modules/multiresvolume/rendering/errorhistogrammanager.h

@@ -0,0 +1,76 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __ERRORHISTOGRAMMANAGER_H__
+#define __ERRORHISTOGRAMMANAGER_H__
+
+#include <fstream>
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+#include <map>
+
+#include <ghoul/glm.h>
+
+namespace openspace {
+
+class ErrorHistogramManager {
+public:
+    ErrorHistogramManager(TSP* tsp);
+    ~ErrorHistogramManager();
+
+    bool buildHistograms(int numBins);
+    const Histogram* getHistogram(unsigned int brickIndex) const;
+
+    bool loadFromFile(const std::string& filename);
+    bool saveToFile(const std::string& filename);
+
+private:
+    TSP* _tsp;
+    std::ifstream* _file;
+
+    std::vector<Histogram> _histograms;
+    unsigned int _numInnerNodes;
+    float _minBin;
+    float _maxBin;
+    int _numBins;
+
+    std::map<unsigned int, std::vector<float>> _voxelCache;
+
+    bool buildFromLeaf(unsigned int bstOffset, unsigned int octreeOffset);
+    std::vector<float> readValues(unsigned int brickIndex) const;
+
+    int parentOffset(int offset, int base) const;
+
+    unsigned int brickToInnerNodeIndex(unsigned int brickIndex) const;
+    unsigned int innerNodeToBrickIndex(unsigned int innerNodeIndex) const;
+    unsigned int linearCoords(glm::vec3 coords) const;
+    unsigned int linearCoords(int x, int y, int z) const;
+    unsigned int linearCoords(glm::ivec3 coords) const;
+
+    float interpolate(glm::vec3 samplePoint, const std::vector<float>& voxels) const;
+};
+
+} // namespace openspace
+
+#endif // __ERRORHISTOGRAMMANAGER_H__

modules/multiresvolume/rendering/histogram.cpp

@@ -0,0 +1,221 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/histogram.h>
+
+#include <ghoul/logging/logmanager.h>
+
+#include <cmath>
+#include <cassert>
+namespace {
+    const std::string _loggerCat = "Histogram";
+}
+
+namespace openspace {
+
+Histogram::Histogram()
+    : _minBin(0)
+    , _maxBin(0)
+    , _numBins(-1)
+    , _data(nullptr) {}
+
+Histogram::Histogram(float minBin, float maxBin, int numBins)
+    : _minBin(minBin)
+    , _maxBin(maxBin)
+    , _numBins(numBins)
+    , _data(nullptr) {
+
+    _data = new float[numBins];
+    for (int i = 0; i < numBins; ++i) {
+        _data[i] = 0.0;
+    }
+}
+
+Histogram::Histogram(float minBin, float maxBin, int numBins, float *data)
+    : _minBin(minBin)
+    , _maxBin(maxBin)
+    , _numBins(numBins)
+    , _data(data) {}
+
+Histogram::Histogram(Histogram&& other) {
+    _minBin = other._minBin;
+    _maxBin = other._maxBin;
+    _numBins = other._numBins;
+    _data = other._data;
+    other._data = nullptr;
+}
+
+Histogram& Histogram::operator=(Histogram&& other) {
+    _minBin = other._minBin;
+    _maxBin = other._maxBin;
+    _numBins = other._numBins;
+    _data = other._data;
+    other._data = nullptr;
+	return *this;
+}
+
+
+Histogram::~Histogram() {
+    if (_data) {
+        delete[] _data;
+    }
+}
+
+
+int Histogram::numBins() const {
+    return _numBins;
+}
+
+float Histogram::minBin() const {
+    return _minBin;
+}
+
+float Histogram::maxBin() const {
+    return _maxBin;
+}
+
+bool Histogram::isValid() const {
+    return _numBins != -1;
+}
+
+
+bool Histogram::add(float bin, float value) {
+    if (bin < _minBin || bin > _maxBin) {
+        // Out of range
+        return false;
+    }
+
+    float normalizedBin = (bin - _minBin) / (_maxBin - _minBin);    // [0.0, 1.0]
+    int binIndex = floor(normalizedBin * _numBins);                 // [0, _numBins]
+    if (binIndex == _numBins) binIndex--;                           // [0, _numBins[
+
+    _data[binIndex] += value;
+    return true;
+}
+
+bool Histogram::add(const Histogram& histogram) {
+    if (_minBin == histogram.minBin() && _maxBin == histogram.maxBin() && _numBins == histogram.numBins()) {
+
+        const float* data = histogram.data();
+        for (int i = 0; i < _numBins; i++) {
+
+            _data[i] += data[i];
+
+        }
+        return true;
+    } else {
+        LERROR("Dimension mismatch");
+        return false;
+    }
+}
+
+bool Histogram::addRectangle(float lowBin, float highBin, float value) {
+    if (lowBin == highBin) return true;
+    if (lowBin > highBin) {
+        std::swap(lowBin, highBin);
+    }
+    if (lowBin < _minBin || highBin > _maxBin) {
+        // Out of range
+        return false;
+    }
+
+    float normalizedLowBin = (lowBin - _minBin) / (_maxBin - _minBin);
+    float normalizedHighBin = (highBin - _minBin) / (_maxBin - _minBin);
+
+    float lowBinIndex = normalizedLowBin * _numBins;
+    float highBinIndex = normalizedHighBin * _numBins;
+
+    int fillLow = floor(lowBinIndex);
+    int fillHigh = ceil(highBinIndex);
+
+    for (int i = fillLow; i < fillHigh; i++) {
+        _data[i] += value;
+    }
+
+    if (lowBinIndex > fillLow) {
+        float diff = lowBinIndex - fillLow;
+        _data[fillLow] -= diff * value;
+    }
+    if (highBinIndex < fillHigh) {
+        float diff = -highBinIndex + fillHigh;
+        _data[fillHigh - 1] -= diff * value;
+    }
+
+    return true;
+}
+
+
+float Histogram::interpolate(float bin) const {
+    float normalizedBin = (bin - _minBin) / (_maxBin - _minBin);
+    float binIndex = normalizedBin * _numBins - 0.5; // Center
+
+    float interpolator = binIndex - floor(binIndex);
+    int binLow = floor(binIndex);
+    int binHigh = ceil(binIndex);
+
+    // Clamp bins
+    if (binLow < 0) binLow = 0;
+    if (binHigh >= _numBins) binHigh = _numBins - 1;
+
+    return (1.0 - interpolator) * _data[binLow] + interpolator * _data[binHigh];
+}
+
+float Histogram::sample(int binIndex) const {
+    assert(binIndex >= 0 && binIndex < _numBins);
+    return _data[binIndex];
+}
+
+
+const float* Histogram::data() const {
+    return _data;
+}
+
+std::vector<std::pair<float,float>> Histogram::getDecimated(int numBins) const {
+    // Return a copy of _data decimated as in Ljung 2004
+    return std::vector<std::pair<float,float>>();
+}
+
+
+void Histogram::normalize() {
+    float sum = 0.0;
+    for (int i = 0; i < _numBins; i++) {
+        sum += _data[i];
+    }
+    for (int i = 0; i < _numBins; i++) {
+        _data[i] /= sum;
+    }
+}
+
+void Histogram::print() const {
+    std::cout << "number of bins: " << _numBins << std::endl
+              << "range: " << _minBin << " - " << _maxBin << std::endl << std::endl;
+    for (int i = 0; i < _numBins; i++) {
+        float low = _minBin + float(i) / _numBins * (_maxBin - _minBin);
+        float high = low + (_maxBin - _minBin) / float(_numBins);
+        std::cout << "[" << low << ", " << high << "[" << std::endl
+                  << "   " << _data[i] << std::endl;
+    }
+}
+
+}

modules/multiresvolume/rendering/histogram.h

@@ -0,0 +1,70 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __HISTOGRAM_H__
+#define __HISTOGRAM_H__
+
+#include <vector>
+#include <iostream>
+
+namespace openspace {
+class Histogram {
+
+public:
+    Histogram();
+    Histogram(float minBin, float maxBin, int numBins);
+    Histogram(float minBin, float maxBin, int numBins, float *data);
+    Histogram(Histogram&& other);
+    ~Histogram();
+
+    Histogram& operator=(Histogram&& other);
+
+    int numBins() const;
+    float minBin() const;
+    float maxBin() const;
+    bool isValid() const;
+
+    bool add(float bin, float value);
+    bool add(const Histogram& histogram);
+    bool addRectangle(float lowBin, float highBin, float value);
+
+    float interpolate(float bin) const;
+    float sample(int binIndex) const;
+    const float* data() const;
+    std::vector<std::pair<float,float>> getDecimated(int numBins) const;
+
+    void normalize();
+    void print() const;
+
+private:
+    int _numBins;
+    float _minBin;
+    float _maxBin;
+
+    float* _data;
+
+}; // class Histogram
+}  // namespace openspace
+
+#endif //__HISTOGRAM_H__

modules/multiresvolume/rendering/histogrammanager.cpp

@@ -0,0 +1,192 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <float.h>
+#include <math.h>
+#include <string.h>
+#include <cassert>
+
+#include <modules/multiresvolume/rendering/histogrammanager.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+
+namespace {
+    const std::string _loggerCat = "HistogramManager";
+}
+
+namespace openspace {
+
+HistogramManager::HistogramManager() {}
+
+HistogramManager::~HistogramManager() {}
+
+bool HistogramManager::buildHistograms(TSP* tsp, int numBins) {
+    std::cout << "Build histograms with " << numBins << " bins each" << std::endl;
+    _numBins = numBins;
+
+    std::ifstream& file = tsp->file();
+    if (!file.is_open()) {
+        return false;
+    }
+    _minBin = 0.0; // Should be calculated from tsp file
+    _maxBin = 1.0; // Should be calculated from tsp file
+
+    int numTotalNodes = tsp->numTotalNodes();
+    _histograms = std::vector<Histogram>(numTotalNodes);
+
+    bool success = buildHistogram(tsp, 0);
+
+    return success;
+}
+
+Histogram* HistogramManager::getHistogram(unsigned int brickIndex) {
+    return &_histograms[brickIndex];
+}
+
+bool HistogramManager::buildHistogram(TSP* tsp, unsigned int brickIndex) {
+    Histogram histogram(_minBin, _maxBin, _numBins);
+    bool isBstLeaf = tsp->isBstLeaf(brickIndex);
+    bool isOctreeLeaf = tsp->isOctreeLeaf(brickIndex);
+
+    if (isBstLeaf && isOctreeLeaf) {
+        // TSP leaf, read from file and build histogram
+        std::vector<float> voxelValues = readValues(tsp, brickIndex);
+        unsigned int numVoxels = voxelValues.size();
+
+        for (unsigned int v = 0; v < numVoxels; ++v) {
+            histogram.add(voxelValues[v], 1.0);
+        }
+    } else {
+        // Has children
+        auto children = std::vector<unsigned int>();
+
+        if (!isBstLeaf) {
+            // Push BST children
+            children.push_back(tsp->getBstLeft(brickIndex));
+            children.push_back(tsp->getBstRight(brickIndex));
+        }
+        if (!isOctreeLeaf) {
+            // Push Octree children
+            unsigned int firstChild = tsp->getFirstOctreeChild(brickIndex);
+            for (int c = 0; c < 8; c++) {
+                children.push_back(firstChild + c);
+            }
+        }
+        int numChildren = children.size();
+        for (int c = 0; c < numChildren; c++) {
+            // Visit child
+            unsigned int childIndex = children[c];
+            if (_histograms[childIndex].isValid() || buildHistogram(tsp, childIndex)) {
+                if (numChildren <= 8 || c < 2) {
+                    // If node has both BST and Octree children, only add BST ones
+                    histogram.add(_histograms[childIndex]);
+                }
+            } else {
+                return false;
+            }
+        }
+    }
+
+    //histogram.normalize();
+    _histograms[brickIndex] = std::move(histogram);
+
+
+    return true;
+}
+
+std::vector<float> HistogramManager::readValues(TSP* tsp, unsigned int brickIndex) {
+    unsigned int paddedBrickDim = tsp->paddedBrickDim();
+    unsigned int numBrickVals = paddedBrickDim * paddedBrickDim * paddedBrickDim;
+    std::vector<float> voxelValues(numBrickVals);
+
+    std::streampos offset = tsp->dataPosition() + static_cast<long long>(brickIndex*numBrickVals*sizeof(float));
+    std::ifstream& file = tsp->file();
+    file.seekg(offset);
+
+    file.read(reinterpret_cast<char*>(&voxelValues[0]),
+        static_cast<size_t>(numBrickVals)*sizeof(float));
+
+    return voxelValues;
+}
+
+bool HistogramManager::loadFromFile(const std::string& filename) {
+    std::ifstream file(filename, std::ios::in | std::ios::binary);
+    if (!file.is_open()) {
+	return false;
+    }
+
+    int numHistograms;
+    file.read(reinterpret_cast<char*>(&numHistograms), sizeof(int));
+    file.read(reinterpret_cast<char*>(&_numBins), sizeof(int));
+    file.read(reinterpret_cast<char*>(&_minBin), sizeof(float));
+    file.read(reinterpret_cast<char*>(&_maxBin), sizeof(float));
+
+    int nFloats = numHistograms * _numBins;
+    float* histogramData = new float[nFloats];
+    file.read(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+
+    _histograms = std::vector<Histogram>(numHistograms);
+
+    for (int i = 0; i < numHistograms; ++i) {
+	int offset = i*_numBins;
+	float* data = new float[_numBins];
+	memcpy(data, &histogramData[offset], sizeof(float) * _numBins);
+	_histograms[i] = Histogram(_minBin, _maxBin, _numBins, data);
+    }
+
+    delete[] histogramData;
+    // No need to deallocate histogram data, since histograms take ownership.
+    file.close();
+    return true;
+}
+
+
+bool HistogramManager::saveToFile(const std::string& filename) {
+    std::ofstream file(filename, std::ios::out | std::ios::binary);
+    if (!file.is_open()) {
+	return false;
+    }
+
+    int numHistograms = _histograms.size();
+    file.write(reinterpret_cast<char*>(&numHistograms), sizeof(int));
+    file.write(reinterpret_cast<char*>(&_numBins), sizeof(int));
+    file.write(reinterpret_cast<char*>(&_minBin), sizeof(float));
+    file.write(reinterpret_cast<char*>(&_maxBin), sizeof(float));
+
+    int nFloats = numHistograms * _numBins;
+    float* histogramData = new float[nFloats];
+
+    for (int i = 0; i < numHistograms; ++i) {
+	int offset = i*_numBins;
+	memcpy(&histogramData[offset], _histograms[i].data(), sizeof(float) * _numBins);
+    }
+
+    file.write(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+    delete[] histogramData;
+
+    file.close();
+    return true;
+}
+
+} // namespace openspace
+

modules/multiresvolume/rendering/histogrammanager.h

@@ -0,0 +1,54 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __HISTOGRAMMANAGER_H__
+#define __HISTOGRAMMANAGER_H__
+
+#include <fstream>
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+
+namespace openspace {
+
+class HistogramManager {
+public:
+    HistogramManager();
+    ~HistogramManager();
+    bool buildHistograms(TSP* tsp, int numBins);
+    Histogram* getHistogram(unsigned int brickIndex);
+    bool loadFromFile(const std::string& filename);
+    bool saveToFile(const std::string& filename);
+private:
+    std::vector<Histogram> _histograms;
+    float _minBin;
+    float _maxBin;
+    int _numBins;
+
+    bool buildHistogram(TSP* tsp, unsigned int brickIndex);
+    std::vector<float> readValues(TSP* tsp, unsigned int brickIndex);
+};
+
+} // namespace openspace
+
+#endif // __HISTOGRAMMANAGER_H__

modules/multiresvolume/rendering/localerrorhistogrammanager.cpp

@@ -0,0 +1,490 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <float.h>
+#include <map>
+#include <string.h>
+#include <cmath>
+
+#include <modules/multiresvolume/rendering/localerrorhistogrammanager.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+
+#include <openspace/util/progressbar.h>
+
+#include <ghoul/logging/logmanager.h>
+
+namespace {
+    const std::string _loggerCat = "LocalErrorHistogramManager";
+}
+
+namespace openspace {
+
+LocalErrorHistogramManager::LocalErrorHistogramManager(TSP* tsp) : _tsp(tsp) {}
+
+LocalErrorHistogramManager::~LocalErrorHistogramManager() {}
+
+bool LocalErrorHistogramManager::buildHistograms(int numBins) {
+    LINFO("Build histograms with " << numBins << " bins each");
+    _numBins = numBins;
+
+    _file = &(_tsp->file());
+    if (!_file->is_open()) {
+        return false;
+    }
+    _minBin = 0.0; // Should be calculated from tsp file
+    _maxBin = 1.0; // Should be calculated from tsp file as (maxValue - minValue)
+
+    unsigned int numOtLevels = _tsp->numOTLevels();
+    unsigned int numOtLeaves = pow(8, numOtLevels - 1);
+    unsigned int numBstLeaves = pow(2, _tsp->numBSTLevels() - 1);
+
+    _numInnerNodes = _tsp->numTotalNodes() - numOtLeaves * numBstLeaves;
+
+    _spatialHistograms = std::vector<Histogram>(_numInnerNodes);
+    _temporalHistograms = std::vector<Histogram>(_numInnerNodes);
+    for (unsigned int i = 0; i < _numInnerNodes; i++) {
+        _spatialHistograms[i] = Histogram(_minBin, _maxBin, numBins);
+        _temporalHistograms[i] = Histogram(_minBin, _maxBin, numBins);
+    }
+
+    // All TSP Leaves
+    int numOtNodes = _tsp->numOTNodes();
+    int otOffset = (pow(8, numOtLevels - 1) - 1) / 7;
+
+    int numBstNodes = _tsp->numBSTNodes();
+    int bstOffset = numBstNodes / 2;
+
+    int numberOfLeaves = numOtLeaves * numBstLeaves;
+
+    LINFO("Building spatial histograms");
+    ProgressBar pb1(numberOfLeaves);
+    int processedLeaves = 0;
+    pb1.print(processedLeaves);
+    bool success = true;
+    for (int bst = bstOffset; bst < numBstNodes; bst++) {
+        for (int ot = otOffset; ot < numOtNodes; ot++) {
+            success &= buildFromOctreeChild(bst, ot);
+            if (!success) LERROR("Failed in buildFromOctreeChild");
+            if (!success) return false;
+            pb1.print(processedLeaves++);
+        }
+    }
+    //pb1.stop();
+
+
+    LINFO("Building temporal histograms");
+    ProgressBar pb2(numberOfLeaves);
+    processedLeaves = 0;
+    pb2.print(processedLeaves);
+    for (int ot = otOffset; ot < numOtNodes; ot++) {
+        for (int bst = bstOffset; bst < numBstNodes; bst++) {
+            success &= buildFromBstChild(bst, ot);
+            if (!success) LERROR("Failed in buildFromBstChild");
+            if (!success) return false;
+            pb2.print(processedLeaves++);
+        }
+    }
+    //pb2.stop();
+
+    return success;
+}
+
+bool LocalErrorHistogramManager::buildFromOctreeChild(unsigned int bstOffset, unsigned int octreeOffset) {
+    // Add errors to octree parent histogram
+    int numOtNodes = _tsp->numOTNodes();
+    unsigned int childIndex = bstOffset * numOtNodes + octreeOffset;
+    bool isOctreeLeaf = _tsp->isOctreeLeaf(childIndex);
+
+    if (octreeOffset > 0) {
+        // Not octree root
+        std::vector<float> childValues;
+        std::vector<float> parentValues;
+
+        int octreeParent = parentOffset(octreeOffset, 8);
+        unsigned int parentIndex = bstOffset * numOtNodes + octreeParent;
+        unsigned int parentInnerNodeIndex = brickToInnerNodeIndex(parentIndex);
+
+        if (isOctreeLeaf) {
+            childValues = readValues(childIndex);
+        } else {
+            unsigned int childInnerNodeIndex = brickToInnerNodeIndex(childIndex);
+            auto it = _voxelCache.find(childInnerNodeIndex);
+            if (it != _voxelCache.end()) {
+                childValues = it->second;
+            } else {
+                LERROR("Child " << childIndex << " visited without cache, " << bstOffset << ", " << octreeOffset);
+                return false;
+            }
+        }
+
+        int octreeChildIndex = (octreeOffset - 1) % 8;
+        if (octreeChildIndex == 0) {
+            parentValues = readValues(parentIndex);
+            _voxelCache[parentInnerNodeIndex] = parentValues;
+        } else {
+            auto it = _voxelCache.find(parentInnerNodeIndex);
+            if (it != _voxelCache.end()) {
+                parentValues = it->second;
+            } else {
+                LERROR("Parent " << parentIndex << " visited without cache");
+                return false;
+            }
+        }
+
+        // Compare values and add errors to parent histogram
+        unsigned int paddedBrickDim = _tsp->paddedBrickDim();
+        unsigned int brickDim = _tsp->brickDim();
+        unsigned int padding = (paddedBrickDim - brickDim) / 2;
+
+        glm::vec3 parentOffset = glm::vec3(octreeChildIndex % 2, (octreeChildIndex / 2) % 2, octreeChildIndex / 4) * float(brickDim) / 2.f;
+
+        for (int z = 0; z < brickDim; z++) {
+            for (int y = 0; y < brickDim; y++) {
+                for (int x = 0; x < brickDim; x++) {
+                    glm::vec3 childSamplePoint = glm::vec3(x, y, z) + glm::vec3(padding);
+                    glm::vec3 parentSamplePoint = parentOffset + (glm::vec3(x, y, z) + glm::vec3(0.5)) * 0.5f;
+                    float childValue = childValues[linearCoords(childSamplePoint)];
+                    float parentValue = interpolate(parentSamplePoint, parentValues);
+
+                    // Divide by number of child voxels that will be taken into account
+                    float rectangleHeight = std::abs(childValue - parentValue) / 8.0;
+                    _spatialHistograms[parentInnerNodeIndex].addRectangle(childValue, parentValue, rectangleHeight);
+                }
+            }
+        }
+
+        bool isLastOctreeChild = octreeOffset > 0 && octreeChildIndex == 7;
+        if (isLastOctreeChild) {
+            buildFromOctreeChild(bstOffset, octreeParent);
+        }
+    }
+
+    if (!isOctreeLeaf) {
+        unsigned int childInnerNodeIndex = brickToInnerNodeIndex(childIndex);
+        _voxelCache.erase(childInnerNodeIndex);
+    }
+
+    int bstChildIndex = bstOffset % 2;
+    bool isLastBstChild = bstOffset > 0 && bstChildIndex == 0;
+    if (isOctreeLeaf && isLastBstChild) {
+        int bstParent = parentOffset(bstOffset, 2);
+        buildFromOctreeChild(bstParent, octreeOffset);
+    }
+
+    return true;
+}
+
+bool LocalErrorHistogramManager::buildFromBstChild(unsigned int bstOffset, unsigned int octreeOffset) {
+    // Add errors to bst parent histogram
+    int numOtNodes = _tsp->numOTNodes();
+    unsigned int childIndex = bstOffset * numOtNodes + octreeOffset;
+    bool isBstLeaf = _tsp->isBstLeaf(childIndex);
+
+    if (bstOffset > 0) {
+        // Not BST root
+        std::vector<float> childValues;
+        std::vector<float> parentValues;
+
+        int bstParent = parentOffset(bstOffset, 2);
+        unsigned int parentIndex = bstParent * numOtNodes + octreeOffset;
+        unsigned int parentInnerNodeIndex = brickToInnerNodeIndex(parentIndex);
+
+        if (isBstLeaf) {
+            childValues = readValues(childIndex);
+        } else {
+            unsigned int childInnerNodeIndex = brickToInnerNodeIndex(childIndex);
+            auto it = _voxelCache.find(childInnerNodeIndex);
+            if (it != _voxelCache.end()) {
+                childValues = it->second;
+            } else {
+                LERROR("Child " << childIndex << " visited without cache");
+                return false;
+            }
+        }
+
+        int bstChildIndex = bstOffset % 2;
+        if (bstChildIndex == 1) {
+            parentValues = readValues(parentIndex);
+            _voxelCache[parentInnerNodeIndex] = parentValues;
+        } else {
+            auto it = _voxelCache.find(parentInnerNodeIndex);
+            if (it != _voxelCache.end()) {
+                parentValues = it->second;
+            } else {
+                LERROR("Parent " << parentIndex << " visited without cache");
+                return false;
+            }
+        }
+
+        // Compare values and add errors to parent histogram
+        unsigned int paddedBrickDim = _tsp->paddedBrickDim();
+        unsigned int brickDim = _tsp->brickDim();
+        unsigned int padding = (paddedBrickDim - brickDim) / 2;
+
+        for (int z = 0; z < brickDim; z++) {
+            for (int y = 0; y < brickDim; y++) {
+                for (int x = 0; x < brickDim; x++) {
+                    glm::vec3 samplePoint = glm::vec3(x, y, z) + glm::vec3(padding);
+                    unsigned int linearSamplePoint = linearCoords(samplePoint);
+                    float childValue = childValues[linearSamplePoint];
+                    float parentValue = parentValues[linearSamplePoint];
+
+                    // Divide by number of child voxels that will be taken into account
+                    float rectangleHeight = std::abs(childValue - parentValue) / 2.0;
+                    _temporalHistograms[parentInnerNodeIndex].addRectangle(childValue, parentValue, rectangleHeight);
+                }
+            }
+        }
+
+        bool isLastBstChild = bstOffset > 0 && bstChildIndex == 0;
+        if (isLastBstChild) {
+            buildFromBstChild(bstParent, octreeOffset);
+        }
+    }
+
+    if (!isBstLeaf) {
+        unsigned int childInnerNodeIndex = brickToInnerNodeIndex(childIndex);
+        _voxelCache.erase(childInnerNodeIndex);
+    }
+
+    int octreeChildIndex = (octreeOffset - 1) % 8;
+    bool isLastOctreeChild = octreeOffset > 0 && octreeChildIndex == 7;
+    if (isBstLeaf && isLastOctreeChild) {
+        int octreeParent = parentOffset(octreeOffset, 8);
+        buildFromBstChild(bstOffset, octreeParent);
+    }
+
+    return true;
+}
+
+bool LocalErrorHistogramManager::loadFromFile(const std::string& filename) {
+    std::ifstream file(filename, std::ios::in | std::ios::binary);
+    if (!file.is_open()) {
+        return false;
+    }
+
+    file.read(reinterpret_cast<char*>(&_numInnerNodes), sizeof(int));
+    file.read(reinterpret_cast<char*>(&_numBins), sizeof(int));
+    file.read(reinterpret_cast<char*>(&_minBin), sizeof(float));
+    file.read(reinterpret_cast<char*>(&_maxBin), sizeof(float));
+
+    int nFloats = _numInnerNodes * _numBins;
+    float* histogramData = new float[nFloats];
+
+    file.read(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+    _spatialHistograms = std::vector<Histogram>(_numInnerNodes);
+    for (int i = 0; i < _numInnerNodes; ++i) {
+        int offset = i*_numBins;
+        float* data = new float[_numBins];
+        memcpy(data, &histogramData[offset], sizeof(float) * _numBins);
+        _spatialHistograms[i] = Histogram(_minBin, _maxBin, _numBins, data);
+    }
+
+    file.read(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+    _temporalHistograms = std::vector<Histogram>(_numInnerNodes);
+    for (int i = 0; i < _numInnerNodes; ++i) {
+        int offset = i*_numBins;
+        float* data = new float[_numBins];
+        memcpy(data, &histogramData[offset], sizeof(float) * _numBins);
+        _temporalHistograms[i] = Histogram(_minBin, _maxBin, _numBins, data);
+    }
+
+    delete[] histogramData;
+    // No need to deallocate histogram data, since histograms take ownership.
+    file.close();
+    return true;
+}
+
+
+bool LocalErrorHistogramManager::saveToFile(const std::string& filename) {
+    std::ofstream file(filename, std::ios::out | std::ios::binary);
+    if (!file.is_open()) {
+        return false;
+    }
+
+    file.write(reinterpret_cast<char*>(&_numInnerNodes), sizeof(int));
+    file.write(reinterpret_cast<char*>(&_numBins), sizeof(int));
+    file.write(reinterpret_cast<char*>(&_minBin), sizeof(float));
+    file.write(reinterpret_cast<char*>(&_maxBin), sizeof(float));
+
+    int nFloats = _numInnerNodes * _numBins;
+    float* histogramData = new float[nFloats];
+
+    for (int i = 0; i < _numInnerNodes; ++i) {
+        int offset = i*_numBins;
+        memcpy(&histogramData[offset], _spatialHistograms[i].data(), sizeof(float) * _numBins);
+    }
+    file.write(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+
+    for (int i = 0; i < _numInnerNodes; ++i) {
+        int offset = i*_numBins;
+        memcpy(&histogramData[offset], _temporalHistograms[i].data(), sizeof(float) * _numBins);
+    }
+    file.write(reinterpret_cast<char*>(histogramData), sizeof(float) * nFloats);
+
+    delete[] histogramData;
+
+    file.close();
+    return true;
+}
+
+unsigned int LocalErrorHistogramManager::linearCoords(glm::vec3 coords) const {
+    return linearCoords(glm::ivec3(coords));
+}
+
+unsigned int LocalErrorHistogramManager::linearCoords(int x, int y, int z) const {
+    return linearCoords(glm::ivec3(x, y, z));
+}
+
+unsigned int LocalErrorHistogramManager::linearCoords(glm::ivec3 coords) const {
+    unsigned int paddedBrickDim = _tsp->paddedBrickDim();
+    return coords.z * paddedBrickDim * paddedBrickDim + coords.y * paddedBrickDim + coords.x;
+}
+
+float LocalErrorHistogramManager::interpolate(glm::vec3 samplePoint, const std::vector<float>& voxels) const {
+    int lowX = samplePoint.x;
+    int lowY = samplePoint.y;
+    int lowZ = samplePoint.z;
+
+    int highX = ceil(samplePoint.x);
+    int highY = ceil(samplePoint.y);
+    int highZ = ceil(samplePoint.z);
+
+    float interpolatorX = 1.0 - (samplePoint.x - lowX);
+    float interpolatorY = 1.0 - (samplePoint.y - lowY);
+    float interpolatorZ = 1.0 - (samplePoint.z - lowZ);
+
+    float v000 = voxels[linearCoords(lowX, lowY, lowZ)];
+    float v001 = voxels[linearCoords(lowX, lowY, highZ)];
+    float v010 = voxels[linearCoords(lowX, highY, lowZ)];
+    float v011 = voxels[linearCoords(lowX, highY, highZ)];
+    float v100 = voxels[linearCoords(highX, lowY, lowZ)];
+    float v101 = voxels[linearCoords(highX, lowY, highZ)];
+    float v110 = voxels[linearCoords(highX, highY, lowZ)];
+    float v111 = voxels[linearCoords(highX, highY, highZ)];
+
+    float v00 = interpolatorZ * v000 + (1.0 - interpolatorZ) * v001;
+    float v01 = interpolatorZ * v010 + (1.0 - interpolatorZ) * v011;
+    float v10 = interpolatorZ * v100 + (1.0 - interpolatorZ) * v101;
+    float v11 = interpolatorZ * v110 + (1.0 - interpolatorZ) * v111;
+
+    float v0 = interpolatorY * v00 + (1.0 - interpolatorY) * v01;
+    float v1 = interpolatorY * v10 + (1.0 - interpolatorY) * v11;
+
+    return interpolatorX * v0 + (1.0 - interpolatorX) * v1;
+
+}
+
+const Histogram* LocalErrorHistogramManager::getSpatialHistogram(unsigned int brickIndex) const {
+    unsigned int innerNodeIndex = brickToInnerNodeIndex(brickIndex);
+    if (innerNodeIndex < _numInnerNodes) {
+        return &(_spatialHistograms[innerNodeIndex]);
+    } else {
+        return nullptr;
+    }
+}
+
+const Histogram* LocalErrorHistogramManager::getTemporalHistogram(unsigned int brickIndex) const {
+    unsigned int innerNodeIndex = brickToInnerNodeIndex(brickIndex);
+    if (innerNodeIndex < _numInnerNodes) {
+        return &(_temporalHistograms[innerNodeIndex]);
+    } else {
+        return nullptr;
+    }
+}
+
+int LocalErrorHistogramManager::parentOffset(int offset, int base) const {
+    if (offset == 0) {
+        return -1;
+    }
+    int depth = floor(log(((base - 1) * offset + 1.0)) / log(base));
+    int firstInLevel = (pow(base, depth) - 1) / (base - 1);
+    int inLevelOffset = offset - firstInLevel;
+
+    int parentDepth = depth - 1;
+    int firstInParentLevel = (pow(base, parentDepth) - 1) / (base - 1);
+    int parentInLevelOffset = inLevelOffset / base;
+
+    int parentOffset = firstInParentLevel + parentInLevelOffset;
+    return parentOffset;
+}
+
+std::vector<float> LocalErrorHistogramManager::readValues(unsigned int brickIndex) const {
+    unsigned int paddedBrickDim = _tsp->paddedBrickDim();
+    unsigned int numBrickVals = paddedBrickDim * paddedBrickDim * paddedBrickDim;
+    std::vector<float> voxelValues(numBrickVals);
+
+    std::streampos offset = _tsp->dataPosition() + static_cast<long long>(brickIndex*numBrickVals*sizeof(float));
+    _file->seekg(offset);
+
+    _file->read(reinterpret_cast<char*>(&voxelValues[0]),
+        static_cast<size_t>(numBrickVals)*sizeof(float));
+
+    return voxelValues;
+}
+
+unsigned int LocalErrorHistogramManager::brickToInnerNodeIndex(unsigned int brickIndex) const {
+    unsigned int numOtNodes = _tsp->numOTNodes();
+    unsigned int numBstLevels = _tsp->numBSTLevels();
+
+    unsigned int numInnerBstNodes = (pow(2, numBstLevels - 1) - 1) * numOtNodes;
+    if (brickIndex < numInnerBstNodes) return brickIndex;
+
+    unsigned int numOtLeaves = pow(8, _tsp->numOTLevels() - 1);
+    unsigned int numOtInnerNodes = (numOtNodes - numOtLeaves);
+
+    unsigned int innerBstOffset = brickIndex - numInnerBstNodes;
+    unsigned int rowIndex = innerBstOffset / numOtNodes;
+    unsigned int indexInRow = innerBstOffset % numOtNodes;
+
+    if (indexInRow >= numOtInnerNodes) return _numInnerNodes;
+
+    unsigned int offset = rowIndex * numOtInnerNodes;
+    unsigned int leavesOffset = offset + indexInRow;
+
+    return numInnerBstNodes + leavesOffset;
+}
+
+unsigned int LocalErrorHistogramManager::innerNodeToBrickIndex(unsigned int innerNodeIndex) const {
+    unsigned int numOtNodes = _tsp->numOTNodes();
+    unsigned int numBstLevels = _tsp->numBSTLevels();
+
+    unsigned int numInnerBstNodes = (pow(2, numBstLevels - 1) - 1) * numOtNodes;
+    if (innerNodeIndex < numInnerBstNodes) return innerNodeIndex;
+
+    unsigned int numOtLeaves = pow(8, _tsp->numOTLevels() - 1);
+    unsigned int numOtInnerNodes = (numOtNodes - numOtLeaves);
+
+    unsigned int innerBstOffset = innerNodeIndex - numInnerBstNodes;
+    unsigned int rowIndex = innerBstOffset / numOtInnerNodes;
+    unsigned int indexInRow = innerBstOffset % numOtInnerNodes;
+
+    unsigned int offset = rowIndex * numOtNodes;
+    unsigned int leavesOffset = offset + indexInRow;
+
+    return numInnerBstNodes + leavesOffset;
+}
+
+} // namespace openspace
+

modules/multiresvolume/rendering/localerrorhistogrammanager.h

@@ -0,0 +1,80 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __LOCALERRORHISTOGRAMMANAGER_H__
+#define __LOCALERRORHISTOGRAMMANAGER_H__
+
+#include <fstream>
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+#include <map>
+
+#include <ghoul/glm.h>
+
+namespace openspace {
+
+class LocalErrorHistogramManager {
+public:
+    LocalErrorHistogramManager(TSP* tsp);
+    ~LocalErrorHistogramManager();
+
+    bool buildHistograms(int numBins);
+    const Histogram* getSpatialHistogram(unsigned int brickIndex) const;
+    const Histogram* getTemporalHistogram(unsigned int brickIndex) const;
+
+    bool loadFromFile(const std::string& filename);
+    bool saveToFile(const std::string& filename);
+
+private:
+    TSP* _tsp;
+    std::ifstream* _file;
+
+    std::vector<Histogram> _spatialHistograms;
+    std::vector<Histogram> _temporalHistograms;
+    unsigned int _numInnerNodes;
+    float _minBin;
+    float _maxBin;
+    int _numBins;
+
+    std::map<unsigned int, std::vector<float>> _voxelCache;
+
+    bool buildFromOctreeChild(unsigned int bstOffset, unsigned int octreeOffset);
+    bool buildFromBstChild(unsigned int bstOffset, unsigned int octreeOffset);
+
+    std::vector<float> readValues(unsigned int brickIndex) const;
+
+    int parentOffset(int offset, int base) const;
+
+    unsigned int brickToInnerNodeIndex(unsigned int brickIndex) const;
+    unsigned int innerNodeToBrickIndex(unsigned int innerNodeIndex) const;
+    unsigned int linearCoords(glm::vec3 coords) const;
+    unsigned int linearCoords(int x, int y, int z) const;
+    unsigned int linearCoords(glm::ivec3 coords) const;
+
+    float interpolate(glm::vec3 samplePoint, const std::vector<float>& voxels) const;
+};
+
+} // namespace openspace
+
+#endif // __LOCALERRORHISTOGRAMMANAGER_H__

modules/multiresvolume/rendering/localtfbrickselector.cpp

@@ -0,0 +1,348 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/localtfbrickselector.h>
+#include <modules/multiresvolume/rendering/localerrorhistogrammanager.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+#include <openspace/rendering/transferfunction.h>
+#include <algorithm>
+#include <cassert>
+
+namespace {
+    const std::string _loggerCat = "LocalTfBrickSelector";
+}
+
+namespace openspace {
+
+LocalTfBrickSelector::LocalTfBrickSelector(TSP* tsp, LocalErrorHistogramManager* hm, TransferFunction* tf, int memoryBudget, int streamingBudget)
+    : _tsp(tsp)
+    , _histogramManager(hm)
+    , _transferFunction(tf)
+    , _memoryBudget(memoryBudget)
+    , _streamingBudget(streamingBudget) {}
+
+LocalTfBrickSelector::~LocalTfBrickSelector() {}
+
+bool LocalTfBrickSelector::initialize() {
+    return calculateBrickErrors();
+}
+
+void LocalTfBrickSelector::setMemoryBudget(int memoryBudget) {
+    _memoryBudget = memoryBudget;
+}
+
+void LocalTfBrickSelector::setStreamingBudget(int streamingBudget) {
+    _streamingBudget = streamingBudget;
+}
+
+void LocalTfBrickSelector::selectBricks(int timestep, std::vector<int>& bricks) {
+    int numTimeSteps = _tsp->header().numTimesteps_;
+    int numBricksPerDim = _tsp->header().xNumBricks_;
+
+    unsigned int rootNode = 0;
+    BrickSelection::SplitType splitType;
+    float rootSplitPoints = splitPoints(rootNode, splitType);
+    BrickSelection brickSelection = BrickSelection(numBricksPerDim, numTimeSteps, splitType, rootSplitPoints);
+
+    std::vector<BrickSelection> priorityQueue;
+    std::vector<BrickSelection> leafSelections;
+    std::vector<BrickSelection> temporalSplitQueue;
+    std::vector<BrickSelection> deadEnds;
+
+    if (splitType != BrickSelection::SplitType::None) {
+        priorityQueue.push_back(brickSelection);
+    } else {
+        leafSelections.push_back(brickSelection);
+    }
+
+    int memoryBudget = _memoryBudget;
+    int totalStreamingBudget = _streamingBudget * numTimeSteps;
+    int nBricksInMemory = 1;
+    int nStreamedBricks = 1;
+
+    while (nBricksInMemory <= memoryBudget - 7 && priorityQueue.size() > 0) {
+        std::pop_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+        BrickSelection bs = priorityQueue.back();
+
+        // TODO: handle edge case when we can only afford temporal splits or no split (only 1 spot left)
+
+        unsigned int brickIndex = bs.brickIndex;
+        priorityQueue.pop_back();
+        if (bs.splitType == BrickSelection::SplitType::Temporal) {
+            int timeSpanCenter = bs.centerT();
+            unsigned int childBrickIndex;
+            bool pickRightTimeChild = bs.timestepInRightChild(timestep);
+
+            // On average on the whole time period, splitting this spatial brick in two time steps
+            // would generate twice as much streaming. Current number of streams of this spatial brick
+            // is 2^nTemporalSplits over the whole time period.
+            int newStreams = std::pow(2, bs.nTemporalSplits);
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // Reached dead end (streaming budget would be exceeded)
+                deadEnds.push_back(bs);
+                break;
+            }
+            nStreamedBricks += newStreams;
+
+            if (pickRightTimeChild) {
+                childBrickIndex = _tsp->getBstRight(brickIndex);
+            } else {
+                childBrickIndex = _tsp->getBstLeft(brickIndex);
+            }
+
+            BrickSelection::SplitType childSplitType;
+            float childSplitPoints = splitPoints(childBrickIndex, childSplitType);
+            BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, childSplitType, childSplitPoints);
+
+            if (childSplitType != BrickSelection::SplitType::None) {
+                priorityQueue.push_back(childSelection);
+                std::push_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                leafSelections.push_back(childSelection);
+            }
+        } else if (bs.splitType == BrickSelection::SplitType::Spatial) {
+            nBricksInMemory += 7; // Remove one and add eight.
+            unsigned int firstChild = _tsp->getFirstOctreeChild(brickIndex);
+
+            // On average on the whole time period, splitting this spatial brick into eight spatial bricks
+            // would generate eight times as much streaming. Current number of streams of this spatial brick
+            // is 2^nTemporalStreams over the whole time period.
+            int newStreams = 7*std::pow(2, bs.nTemporalSplits);
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // Reached dead end (streaming budget would be exceeded)
+                // However, temporal split might be possible
+                if (bs.splitType != BrickSelection::SplitType::Temporal) {
+                    bs.splitType = BrickSelection::SplitType::Temporal;
+                    bs.splitPoints = temporalSplitPoints(bs.brickIndex);
+                }
+                if (bs.splitPoints > -1) {
+                    temporalSplitQueue.push_back(bs);
+                } else {
+                    deadEnds.push_back(bs);
+                }
+                break;
+            }
+            nStreamedBricks += newStreams;
+
+            for (unsigned int i = 0; i < 8; i++) {
+                unsigned int childBrickIndex = firstChild + i;
+
+                BrickSelection::SplitType childSplitType;
+                float childSplitPoints = splitPoints(childBrickIndex, childSplitType);
+                BrickSelection childSelection = bs.splitSpatially(i % 2, (i/2) % 2, i/4, childBrickIndex, childSplitType, childSplitPoints);
+
+                if (childSplitType != BrickSelection::SplitType::None) {
+                    priorityQueue.push_back(childSelection);
+                    std::push_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+                } else {
+                    leafSelections.push_back(childSelection);
+                }
+            }
+        }
+    }
+
+    // Is it possible that we may stream more bricks?
+    if (nStreamedBricks < totalStreamingBudget) {
+        while (priorityQueue.size() > 0) {
+            BrickSelection bs = priorityQueue.back();
+            if (bs.splitType != BrickSelection::SplitType::Temporal) {
+                bs.splitType = BrickSelection::SplitType::Temporal;
+                bs.splitPoints = temporalSplitPoints(bs.brickIndex);
+            }
+            priorityQueue.pop_back();
+            if (bs.splitPoints > -1) {
+                temporalSplitQueue.push_back(bs);
+                std::push_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                deadEnds.push_back(bs);
+            }
+        }
+
+        while (nStreamedBricks < totalStreamingBudget - 1 && temporalSplitQueue.size() > 0) {
+            std::pop_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            BrickSelection bs = temporalSplitQueue.back();
+            temporalSplitQueue.pop_back();
+
+            unsigned int brickIndex = bs.brickIndex;
+            int newStreams = std::pow(2, bs.nTemporalSplits);
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // The current best choice would make us exceed the streaming budget, try next instead.
+                deadEnds.push_back(bs);
+                continue;
+            }
+
+            nStreamedBricks += newStreams;
+            unsigned int childBrickIndex;
+            bool pickRightTimeChild = bs.timestepInRightChild(timestep);
+
+            if (pickRightTimeChild) {
+                childBrickIndex = _tsp->getBstRight(brickIndex);
+            } else {
+                childBrickIndex = _tsp->getBstLeft(brickIndex);
+            }
+
+            float childSplitPoints = temporalSplitPoints(childBrickIndex);
+
+            if (childSplitPoints > -1) {
+                BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, BrickSelection::SplitType::Temporal, childSplitPoints);
+                temporalSplitQueue.push_back(childSelection);
+                std::push_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, BrickSelection::SplitType::None, -1);
+                deadEnds.push_back(childSelection);
+            }
+        }
+    } else {
+        // Write selected inner nodes to brickSelection vector
+        for (const BrickSelection& bs : priorityQueue) {
+            writeSelection(bs, bricks);
+        } 
+    }
+
+    // Write selected inner nodes to brickSelection vector
+    for (const BrickSelection& bs : temporalSplitQueue) {
+        writeSelection(bs, bricks);
+    }
+    // Write dead end nodes to brickSelection vector
+    for (const BrickSelection& bs : deadEnds) {
+        writeSelection(bs, bricks);
+    }
+    // Write selected leaf nodes to brickSelection vector
+    for (const BrickSelection& bs : leafSelections) {
+        writeSelection(bs, bricks);
+    }
+
+    //std::cout << "Bricks in memory: " << nBricksInMemory << "/" << _memoryBudget << "___\t\t"
+    //          << "Streamed bricks:  " << nStreamedBricks << "/" << totalStreamingBudget << std::flush << "___\r";
+}
+
+float LocalTfBrickSelector::temporalSplitPoints(unsigned int brickIndex) {
+    if (_tsp->isBstLeaf(brickIndex)) {
+        return -1;
+    }
+    return _brickErrors[brickIndex].temporal * 0.5;
+}
+
+float LocalTfBrickSelector::spatialSplitPoints(unsigned int brickIndex) {
+    if (_tsp->isOctreeLeaf(brickIndex)) {
+        return -1;
+    }
+    return _brickErrors[brickIndex].spatial * 0.125;
+}
+
+float LocalTfBrickSelector::splitPoints(unsigned int brickIndex, BrickSelection::SplitType& splitType) {
+    float temporalPoints = temporalSplitPoints(brickIndex);
+    float spatialPoints = spatialSplitPoints(brickIndex);
+
+    float splitPoints;
+
+    if (spatialPoints > 0 && spatialPoints > temporalPoints) {
+        splitPoints = spatialPoints;
+        splitType = BrickSelection::SplitType::Spatial;
+    } else if (temporalPoints > 0) {
+        splitPoints = temporalPoints;
+        splitType = BrickSelection::SplitType::Temporal;
+    } else {
+        splitPoints = -1;
+        splitType = BrickSelection::SplitType::None;
+    }
+    return splitPoints;
+}
+
+
+bool LocalTfBrickSelector::calculateBrickErrors() {
+    TransferFunction *tf = _transferFunction;
+    if (!tf) return false;
+
+    size_t tfWidth = tf->width();
+    if (tfWidth <= 0) return false;
+
+    std::vector<float> gradients(tfWidth - 1);
+    for (size_t offset = 0; offset < tfWidth - 1; offset++) {
+        glm::vec4 prevRgba = tf->sample(offset);
+        glm::vec4 nextRgba = tf->sample(offset + 1);
+
+        float colorDifference = glm::distance(prevRgba, nextRgba);
+        float alpha = (prevRgba.w + nextRgba.w) * 0.5;
+
+        gradients[offset] = colorDifference*alpha;
+    }
+
+    unsigned int nHistograms = _tsp->numTotalNodes();
+    _brickErrors = std::vector<Error>(nHistograms);
+
+    for (unsigned int brickIndex = 0; brickIndex < nHistograms; brickIndex++) {
+        if (_tsp->isOctreeLeaf(brickIndex)) {
+            _brickErrors[brickIndex].spatial = 0.0;
+        } else {
+            const Histogram* histogram = _histogramManager->getSpatialHistogram(brickIndex);
+            float error = 0;
+            for (int i = 0; i < gradients.size(); i++) {
+                float x = (i + 0.5) / tfWidth;
+                float sample = histogram->interpolate(x);
+                assert(sample >= 0);
+                assert(gradients[i] >= 0);
+                error += sample * gradients[i];
+            }
+            _brickErrors[brickIndex].spatial = error;
+        }
+
+        if (_tsp->isBstLeaf(brickIndex)) {
+            _brickErrors[brickIndex].temporal = 0.0;
+        } else {
+            const Histogram* histogram = _histogramManager->getTemporalHistogram(brickIndex);
+            float error = 0;
+            for (int i = 0; i < gradients.size(); i++) {
+                float x = (i + 0.5) / tfWidth;
+                float sample = histogram->interpolate(x);
+                assert(sample >= 0);
+                assert(gradients[i] >= 0);
+                error += sample * gradients[i];
+            }
+            _brickErrors[brickIndex].temporal = error;
+        }
+    }
+
+    return true;
+}
+
+int LocalTfBrickSelector::linearCoords(int x, int y, int z) {
+    const TSP::Header &header = _tsp->header();
+    return x + (header.xNumBricks_ * y) + (header.xNumBricks_ * header.yNumBricks_ * z);
+}
+
+void LocalTfBrickSelector::writeSelection(BrickSelection brickSelection, std::vector<int>& bricks) {
+    BrickCover coveredBricks = brickSelection.cover;
+    for (int z = coveredBricks.lowZ; z < coveredBricks.highZ; z++) {
+        for (int y = coveredBricks.lowY; y < coveredBricks.highY; y++) {
+            for (int x = coveredBricks.lowX; x < coveredBricks.highX; x++) {
+                bricks[linearCoords(x, y, z)] = brickSelection.brickIndex;
+            }
+        }
+    }
+}
+
+
+} // namespace openspace

modules/multiresvolume/rendering/localtfbrickselector.h

@@ -0,0 +1,78 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __LOCALTFBRICKSELECTOR_H__
+#define __LOCALTFBRICKSELECTOR_H__
+
+#include <vector>
+#include <modules/multiresvolume/rendering/brickselection.h>
+#include <modules/multiresvolume/rendering/brickselector.h>
+#include <modules/multiresvolume/rendering/brickcover.h>
+
+
+namespace openspace {
+
+class TSP;
+class LocalErrorHistogramManager;
+class TransferFunction;
+
+class LocalTfBrickSelector : public BrickSelector {
+public:
+    struct Error {
+        float spatial;
+        float temporal;
+    };
+
+    LocalTfBrickSelector(TSP* tsp, LocalErrorHistogramManager* hm, TransferFunction* tf, int memoryBudget, int streamingBudget);
+    ~LocalTfBrickSelector();
+
+    virtual bool initialize();
+
+    void selectBricks(int timestep, std::vector<int>& bricks);
+    void setMemoryBudget(int memoryBudget);
+    void setStreamingBudget(int streamingBudget);
+    bool calculateBrickErrors();
+ private:
+
+    TSP* _tsp;
+    LocalErrorHistogramManager* _histogramManager;
+    TransferFunction* _transferFunction;
+    std::vector<Error> _brickErrors;
+
+    float spatialSplitPoints(unsigned int brickIndex);
+    float temporalSplitPoints(unsigned int brickIndex);
+    float splitPoints(unsigned int brickIndex, BrickSelection::SplitType& splitType);
+
+    int linearCoords(int x, int y, int z);
+    void writeSelection(BrickSelection coveredBricks, std::vector<int>& bricks);
+
+    int _memoryBudget;
+    int _streamingBudget;
+    
+};
+
+} // namespace openspace
+
+#endif // __LOCALTFBRICKSELECTOR_H__
+

modules/multiresvolume/rendering/multiresvolumeraycaster.cpp

@@ -0,0 +1,158 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2016                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <sstream>
+#include <memory>
+
+#include <modules/multiresvolume/rendering/multiresvolumeraycaster.h>
+#include <glm/glm.hpp>
+#include <ghoul/opengl/ghoul_gl.h>
+#include <ghoul/opengl/bufferbinding.h>
+
+
+#include <ghoul/opengl/programobject.h>
+#include <openspace/util/powerscaledcoordinate.h>
+#include <openspace/util/updatestructures.h>
+#include <openspace/rendering/renderable.h>
+#include <ghoul/opengl/texture.h>
+
+
+namespace {
+    const std::string GlslRaycastPath = "${MODULES}/multiresvolume/shaders/raycast.glsl";
+    const std::string GlslHelperPath = "${MODULES}/multiresvolume/shaders/helper.glsl";
+    const std::string GlslBoundsVsPath = "${MODULES}/multiresvolume/shaders/boundsVs.glsl";
+    const std::string GlslBoundsFsPath = "${MODULES}/multiresvolume/shaders/boundsFs.glsl";
+}
+
+namespace openspace {
+
+MultiresVolumeRaycaster::MultiresVolumeRaycaster(std::shared_ptr<TSP> tsp,
+    std::shared_ptr<AtlasManager> atlasManager,
+    std::shared_ptr<TransferFunction> transferFunction)
+    : _tsp(tsp)
+    , _atlasManager(atlasManager)
+    , _transferFunction(transferFunction)
+    , _boundingBox(glm::vec3(1.0)) {}
+    
+MultiresVolumeRaycaster::~MultiresVolumeRaycaster() {}
+
+void MultiresVolumeRaycaster::initialize() {
+    _boundingBox.initialize();
+}
+    
+void MultiresVolumeRaycaster::deinitialize() {
+}
+    
+void MultiresVolumeRaycaster::renderEntryPoints(const RenderData& data, ghoul::opengl::ProgramObject& program) {
+    program.setUniform("modelTransform", _modelTransform);
+    program.setUniform("viewProjection", data.camera.viewProjectionMatrix());
+    Renderable::setPscUniforms(program, data.camera, data.position);
+    
+    // Cull back face
+    glEnable(GL_CULL_FACE);
+    glCullFace(GL_BACK);
+    
+    // Render bounding geometry
+    _boundingBox.render();
+}
+    
+void MultiresVolumeRaycaster::renderExitPoints(const RenderData& data, ghoul::opengl::ProgramObject& program) {   
+    // Uniforms
+    program.setUniform("modelTransform", _modelTransform);
+    program.setUniform("viewProjection", data.camera.viewProjectionMatrix());
+    Renderable::setPscUniforms(program, data.camera, data.position);
+
+    // Cull front face
+    glEnable(GL_CULL_FACE);
+    glCullFace(GL_FRONT);
+    
+    // Render bounding geometry
+    _boundingBox.render();
+    
+    // Restore defaults
+    glCullFace(GL_BACK);
+}
+    
+void MultiresVolumeRaycaster::preRaycast(const RaycastData& data, ghoul::opengl::ProgramObject& program) {
+    std::string timestepUniformName = "timestep" + std::to_string(data.id);
+    std::string stepSizeUniformName = "maxStepSize" + std::to_string(data.id);
+
+    std::string id = std::to_string(data.id);
+    //program.setUniform("opacity_" + std::to_string(id), visible ? 1.0f : 0.0f);
+    program.setUniform("stepSizeCoefficient_" + id, _stepSizeCoefficient);
+
+    _tfUnit = std::make_unique<ghoul::opengl::TextureUnit>();
+    _tfUnit->activate();
+    _transferFunction->getTexture().bind();
+    program.setUniform("transferFunction_" + id, _tfUnit->unitNumber());
+    
+    _atlasUnit = std::make_unique<ghoul::opengl::TextureUnit>();
+    _atlasUnit->activate();
+    _atlasManager->textureAtlas().bind();
+    program.setUniform("textureAtlas_" + id, _atlasUnit->unitNumber());
+
+    
+    _atlasMapBinding = std::make_unique<ghoul::opengl::BufferBinding<ghoul::opengl::bufferbinding::Buffer::ShaderStorage>>();
+    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, _atlasMapBinding->bindingNumber(), _atlasManager->atlasMapBuffer());
+    program.setSsboBinding("atlasMapBlock_" + id, _atlasMapBinding->bindingNumber());
+
+    
+    program.setUniform("gridType_" + id, static_cast<int>(_tsp->header().gridType_));
+    program.setUniform("maxNumBricksPerAxis_" + id, static_cast<unsigned int>(_tsp->header().xNumBricks_));
+    program.setUniform("paddedBrickDim_" + id, static_cast<unsigned int>(_tsp->paddedBrickDim()));
+    
+    glm::size3_t size = _atlasManager->textureSize();
+    glm::ivec3 atlasSize(size.x, size.y, size.z);
+    program.setUniform("atlasSize_" + id, atlasSize);
+}
+    
+void MultiresVolumeRaycaster::postRaycast(const RaycastData& data, ghoul::opengl::ProgramObject& program) {
+    // For example: release texture units
+}
+    
+std::string MultiresVolumeRaycaster::getBoundsVsPath() const {
+    return GlslBoundsVsPath;
+}
+    
+std::string MultiresVolumeRaycaster::getBoundsFsPath() const {
+    return GlslBoundsFsPath;
+}
+
+std::string MultiresVolumeRaycaster::getRaycastPath() const {
+    return GlslRaycastPath;
+}
+
+std::string MultiresVolumeRaycaster::getHelperPath() const {
+    return GlslHelperPath; // no helper file
+}
+
+void MultiresVolumeRaycaster::setModelTransform(glm::mat4 transform) {
+    _modelTransform = transform;
+}
+    
+void MultiresVolumeRaycaster::setStepSizeCoefficient(float stepSizeCoefficient) {
+    _stepSizeCoefficient = stepSizeCoefficient;
+}
+    
+}

modules/multiresvolume/rendering/multiresvolumeraycaster.h

@@ -0,0 +1,97 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2016                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __MULTIRESVOLUMERAYCASTER_H__
+#define __MULTIRESVOLUMERAYCASTER_H__
+
+#include <ghoul/glm.h>
+#include <string>
+#include <vector>
+#include <openspace/rendering/volumeraycaster.h>
+#include <openspace/util/boxgeometry.h>
+#include <openspace/util/blockplaneintersectiongeometry.h>
+#include <openspace/rendering/transferfunction.h>
+#include <ghoul/opengl/textureunit.h>
+#include <ghoul/opengl/bufferbinding.h>
+
+#include <modules/multiresvolume/rendering/atlasmanager.h>
+#include <modules/multiresvolume/rendering/tsp.h>
+
+
+namespace ghoul {
+    namespace opengl {
+        class Texture;
+        class ProgramObject;
+    }
+}
+
+namespace openspace {
+
+class RenderData;
+class RaycastData;
+
+class MultiresVolumeRaycaster : public VolumeRaycaster {
+public:
+
+    MultiresVolumeRaycaster(std::shared_ptr<TSP> tsp,
+        std::shared_ptr<AtlasManager> atlasManager,
+        std::shared_ptr<TransferFunction> transferFunction);
+
+    virtual ~MultiresVolumeRaycaster();
+
+    void initialize();
+    void deinitialize();
+    void renderEntryPoints(const RenderData& data, ghoul::opengl::ProgramObject& program) override;
+    void renderExitPoints(const RenderData& data, ghoul::opengl::ProgramObject& program) override;
+    void preRaycast(const RaycastData& data, ghoul::opengl::ProgramObject& program) override;
+    void postRaycast(const RaycastData& data, ghoul::opengl::ProgramObject& program) override;
+
+    std::string getBoundsVsPath() const override;
+    std::string getBoundsFsPath() const override;
+    std::string getRaycastPath() const override;
+    std::string getHelperPath() const override;
+
+    void setModelTransform(glm::mat4 transform);
+    //void setTime(double time);
+    void setStepSizeCoefficient(float coefficient);
+private:
+    BoxGeometry _boundingBox;
+    glm::mat4 _modelTransform;
+    float _stepSizeCoefficient;
+    double _time;
+
+    std::shared_ptr<TSP> _tsp;
+    std::shared_ptr<AtlasManager> _atlasManager;
+    std::shared_ptr<TransferFunction> _transferFunction;
+
+    std::unique_ptr<ghoul::opengl::TextureUnit> _tfUnit;
+    std::unique_ptr<ghoul::opengl::TextureUnit> _atlasUnit;
+    std::unique_ptr < ghoul::opengl::BufferBinding<ghoul::opengl::bufferbinding::Buffer::ShaderStorage>> _atlasMapBinding;
+
+
+}; // MultiresVolumeRaycaster
+
+} // openspace
+
+#endif  // __MULTIRESVOLUMERAYCASTER_H__ 

modules/multiresvolume/rendering/renderablemultiresvolume.cpp

@@ -0,0 +1,638 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2015                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/renderablemultiresvolume.h>
+
+#include <openspace/engine/configurationmanager.h>
+#include <openspace/engine/openspaceengine.h>
+#include <modules/kameleon/include/kameleonwrapper.h>
+#include <openspace/rendering/renderengine.h>
+#include <openspace/rendering/raycastermanager.h>
+
+#include <glm/glm.hpp>
+
+#include <ghoul/filesystem/filesystem.h>
+#include <ghoul/filesystem/file.h>
+#include <ghoul/filesystem/cachemanager.h>
+
+#include <ghoul/opengl/framebufferobject.h>
+#include <ghoul/opengl/programobject.h>
+#include <ghoul/io/texture/texturereader.h>
+#include <ghoul/opengl/texture.h>
+
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/atlasmanager.h>
+#include <modules/multiresvolume/rendering/shenbrickselector.h>
+#include <modules/multiresvolume/rendering/tfbrickselector.h>
+#include <modules/multiresvolume/rendering/simpletfbrickselector.h>
+#include <modules/multiresvolume/rendering/localtfbrickselector.h>
+
+#include <modules/multiresvolume/rendering/histogrammanager.h>
+#include <modules/multiresvolume/rendering/errorhistogrammanager.h>
+#include <modules/multiresvolume/rendering/localerrorhistogrammanager.h>
+#include <openspace/properties/scalarproperty.h>
+#include <openspace/properties/vectorproperty.h>
+
+#include <openspace/util/time.h>
+#include <openspace/util/spicemanager.h>
+
+#include <algorithm>
+#include <iterator>
+#include <fstream>
+#include <algorithm>
+#include <chrono>
+
+
+
+namespace {
+    const std::string _loggerCat = "RenderableMultiresVolume";
+    const std::string KeyDataSource = "Source";
+    const std::string KeyHints = "Hints";
+    const std::string KeyTransferFunction = "TransferFunction";
+
+    const std::string KeyVolumeName = "VolumeName";
+    const std::string KeyBrickSelector = "BrickSelector";
+    const std::string KeyStartTime = "StartTime";
+    const std::string KeyEndTime = "EndTime";
+    const std::string GlslHelpersPath = "${MODULES}/multiresvolume/shaders/helpers_fs.glsl";
+    const std::string GlslHelperPath = "${MODULES}/multiresvolume/shaders/helper.glsl";
+    const std::string GlslHeaderPath = "${MODULES}/multiresvolume/shaders/header.glsl";
+    bool registeredGlslHelpers = false;
+}
+
+namespace openspace {
+
+RenderableMultiresVolume::RenderableMultiresVolume (const ghoul::Dictionary& dictionary)
+    :  Renderable(dictionary)
+    , _transferFunction(nullptr)
+    , _timestep(0)
+    , _atlasMapSize(0)
+    , _tfBrickSelector(nullptr)
+    , _simpleTfBrickSelector(nullptr)
+    , _localTfBrickSelector(nullptr)
+    , _errorHistogramManager(nullptr)
+    , _histogramManager(nullptr)
+    , _localErrorHistogramManager(nullptr)
+    , _stepSizeCoefficient("stepSizeCoefficient", "Stepsize Coefficient", 1.f, 0.01f, 10.f)
+    , _currentTime("currentTime", "Current Time", 0, 0, 0)
+    , _memoryBudget("memoryBudget", "Memory Budget", 0, 0, 0)
+    , _streamingBudget("streamingBudget", "Streaming Budget", 0, 0, 0)
+    , _useGlobalTime("useGlobalTime", "Global Time", false)
+    , _loop("loop", "Loop", false)
+    , _selectorName("selector", "Brick Selector")
+    , _gatheringStats(false)
+    , _statsToFile("printStats", "Print Stats", false)
+    , _statsToFileName("printStatsFileName", "Stats Filename")
+    , _scalingExponent("scalingExponent", "Scaling Exponent", 1, -10, 20)
+    , _scaling("scaling", "Scaling", glm::vec3(1.0, 1.0, 1.0), glm::vec3(0.0), glm::vec3(10.0))
+    , _translation("translation", "Translation", glm::vec3(0.0, 0.0, 0.0), glm::vec3(0.0), glm::vec3(10.0))
+    , _rotation("rotation", "Euler rotation", glm::vec3(0.0, 0.0, 0.0), glm::vec3(0), glm::vec3(6.28))
+{
+    std::string name;
+    //bool success = dictionary.getValue(constants::scenegraphnode::keyName, name);
+    //assert(success);
+
+    _filename = "";
+    bool success = dictionary.getValue(KeyDataSource, _filename);
+    if (!success) {
+        LERROR("Node '" << name << "' did not contain a valid '" <<  KeyDataSource << "'");
+        return;
+    }
+    _filename = absPath(_filename);
+    if (_filename == "") {
+        return;
+    }
+
+    float scalingExponent, stepSizeCoefficient;
+    glm::vec3 scaling, translation, rotation;
+
+    if (dictionary.getValue("ScalingExponent", scalingExponent)) {
+        _scalingExponent = scalingExponent;
+    }
+    if (dictionary.getValue("Scaling", scaling)) {
+        _scaling = scaling;
+    }
+    if (dictionary.getValue("Translation", translation)) {
+        _translation = translation;
+    }
+    if (dictionary.getValue("Rotation", rotation)) {
+        _rotation = rotation;
+    }
+    if (dictionary.getValue("StepSizeCoefficient", stepSizeCoefficient)) {
+        _stepSizeCoefficient = stepSizeCoefficient;
+    }
+
+    std::string startTimeString, endTimeString;
+    bool hasTimeData = true;
+    hasTimeData &= dictionary.getValue(KeyStartTime, startTimeString);
+    hasTimeData &= dictionary.getValue(KeyEndTime, endTimeString);
+    if (hasTimeData) {
+        _startTime = SpiceManager::ref().ephemerisTimeFromDate(startTimeString);
+        _endTime = SpiceManager::ref().ephemerisTimeFromDate(endTimeString);
+    }
+    if (hasTimeData) {
+        _loop = false;
+    } else {
+	_loop = true;
+        LWARNING("Node " << name << " does not provide valid time information. Viewing one image per frame.");
+    }
+
+
+    _transferFunction = nullptr;
+    _transferFunctionPath = "";
+    success = dictionary.getValue(KeyTransferFunction, _transferFunctionPath);
+    if (!success) {
+        LERROR("Node '" << name << "' did not contain a valid '" <<
+            KeyTransferFunction << "'");
+        return;
+    }
+    _transferFunctionPath = absPath(_transferFunctionPath);
+    _transferFunction = std::make_shared<TransferFunction>(_transferFunctionPath);
+
+    //_pscOffset = psc(glm::vec4(0.0));
+    //_boxScaling = glm::vec3(1.0);
+
+
+    /*if (dictionary.hasKey(KeyBoxScaling)) {
+        glm::vec4 scalingVec4(_boxScaling, _w);
+        success = dictionary.getValue(KeyBoxScaling, scalingVec4);
+        if (success) {
+            _boxScaling = scalingVec4.xyz;
+            _w = scalingVec4.w;
+        }
+        else {
+            success = dictionary.getValue(KeyBoxScaling, _boxScaling);
+            if (!success) {
+                LERROR("Node '" << name << "' did not contain a valid '" <<
+                    KeyBoxScaling << "'");
+                return;
+            }
+        }
+    }*/
+
+    //setBoundingSphere(PowerScaledScalar::CreatePSS(glm::length(_boxScaling)*pow(10,_w)));
+
+    _tsp = std::make_shared<TSP>(_filename);
+    _atlasManager = std::make_shared<AtlasManager>(_tsp.get());
+
+    _selectorName = "tf";
+    std::string brickSelectorType;
+    if (dictionary.hasKey(KeyBrickSelector)) {
+        success = dictionary.getValue(KeyBrickSelector, brickSelectorType);
+        if (success) {
+            _selectorName = brickSelectorType;
+        }
+    }
+
+    std::string selectorName = _selectorName;
+    if (selectorName == "simple") {
+        _selector = Selector::SIMPLE;
+    } else if (selectorName == "local") {
+        _selector = Selector::LOCAL;
+    } else {
+        _selector = Selector::TF;
+    }
+
+    addProperty(_selectorName);
+    _selectorName.onChange([&] {
+        Selector s;
+        std::string newSelectorName = _selectorName;
+        if (newSelectorName == "simple") {
+            s = Selector::SIMPLE;
+        } else if (newSelectorName == "local") {
+            s = Selector::LOCAL;
+        } else if (newSelectorName == "tf") {
+            s = Selector::TF;
+        } else {
+            return;
+        }
+        setSelectorType(s);
+    });
+
+    addProperty(_stepSizeCoefficient);
+    addProperty(_useGlobalTime);
+    addProperty(_loop);
+    addProperty(_statsToFile);
+    addProperty(_statsToFileName);
+    addProperty(_scaling);
+    addProperty(_scalingExponent);
+    addProperty(_translation);
+    addProperty(_rotation);
+
+
+    //_brickSelector = new ShenBrickSelector(_tsp, -1, -1);
+}
+
+RenderableMultiresVolume::~RenderableMultiresVolume() {
+    //OsEng.renderEngine()->aBuffer()->removeVolume(this);
+
+    if (_tfBrickSelector)
+        delete _tfBrickSelector;
+    if (_simpleTfBrickSelector)
+        delete _simpleTfBrickSelector;
+    if (_localTfBrickSelector)
+        delete _localTfBrickSelector;
+
+    if (_errorHistogramManager)
+        delete _errorHistogramManager;
+    if (_histogramManager)
+        delete _histogramManager;
+    if (_localErrorHistogramManager)
+        delete _localErrorHistogramManager;
+}
+
+bool RenderableMultiresVolume::setSelectorType(Selector selector) {
+    _selector = selector;
+    switch (_selector) {
+        case Selector::TF:
+            if (!_tfBrickSelector) {
+                TfBrickSelector* tbs;
+                _errorHistogramManager = new ErrorHistogramManager(_tsp.get());
+                _tfBrickSelector = tbs = new TfBrickSelector(_tsp.get(), _errorHistogramManager, _transferFunction.get(), _memoryBudget, _streamingBudget);
+                _transferFunction->setCallback([tbs](const TransferFunction &tf) {
+                    tbs->calculateBrickErrors();
+                });
+                if (initializeSelector()) {
+                    tbs->calculateBrickErrors();
+                    return true;
+                }
+            }
+            break;
+
+        case Selector::SIMPLE:
+            if (!_simpleTfBrickSelector) {
+                SimpleTfBrickSelector *stbs;
+                _histogramManager = new HistogramManager();
+                _simpleTfBrickSelector = stbs = new SimpleTfBrickSelector(_tsp.get(), _histogramManager, _transferFunction.get(), _memoryBudget, _streamingBudget);
+                _transferFunction->setCallback([stbs](const TransferFunction &tf) {
+                    stbs->calculateBrickImportances();
+                });
+                if (initializeSelector()) {
+                    stbs->calculateBrickImportances();
+                    return true;
+                }
+            }
+            break;
+
+        case Selector::LOCAL:
+            if (!_localTfBrickSelector) {
+                LocalTfBrickSelector* ltbs;
+                _localErrorHistogramManager = new LocalErrorHistogramManager(_tsp.get());
+                _localTfBrickSelector = ltbs = new LocalTfBrickSelector(_tsp.get(), _localErrorHistogramManager, _transferFunction.get(), _memoryBudget, _streamingBudget);
+                _transferFunction->setCallback([ltbs](const TransferFunction &tf) {
+                    ltbs->calculateBrickErrors();
+                });
+                if (initializeSelector()) {
+                    ltbs->calculateBrickErrors();
+                    return true;
+                }
+            }
+            break;
+    }
+    return false;
+}
+
+bool RenderableMultiresVolume::initialize() {
+
+    bool success = _tsp && _tsp->load();
+
+    unsigned int maxNumBricks = _tsp->header().xNumBricks_ * _tsp->header().yNumBricks_ * _tsp->header().zNumBricks_;
+
+    unsigned int maxInitialBudget = 2048;
+    int initialBudget = std::min(maxInitialBudget, maxNumBricks);
+
+    _currentTime = properties::IntProperty("currentTime", "Current Time", 0, 0, _tsp->header().numTimesteps_ - 1);
+    _memoryBudget = properties::IntProperty("memoryBudget", "Memory Budget", initialBudget, 0, maxNumBricks);
+    _streamingBudget = properties::IntProperty("streamingBudget", "Streaming Budget", initialBudget, 0, maxNumBricks);
+    addProperty(_currentTime);
+    addProperty(_memoryBudget);
+    addProperty(_streamingBudget);
+
+    if (success) {
+        _brickIndices.resize(maxNumBricks, 0);
+        success &= setSelectorType(_selector);
+    }
+
+    success &= _atlasManager && _atlasManager->initialize();
+
+    _transferFunction->update();
+
+    success &= isReady();
+
+    _raycaster = std::make_unique<MultiresVolumeRaycaster>(_tsp, _atlasManager, _transferFunction);
+    _raycaster->initialize();
+
+    OsEng.renderEngine().raycasterManager().attachRaycaster(*_raycaster.get());
+
+    std::function<void(bool)> onChange = [&](bool enabled) {
+        if (enabled) {
+            OsEng.renderEngine().raycasterManager().attachRaycaster(*_raycaster.get());
+        }
+        else {
+            OsEng.renderEngine().raycasterManager().detachRaycaster(*_raycaster.get());
+        }
+    };
+
+
+
+    return success;
+}
+
+bool RenderableMultiresVolume::deinitialize() {
+    _tsp = nullptr;
+    _transferFunction = nullptr;
+    return true;
+}
+
+bool RenderableMultiresVolume::isReady() const {
+    return true;
+}
+
+
+bool RenderableMultiresVolume::initializeSelector() {
+    int nHistograms = 50;
+    bool success = true;
+
+    switch (_selector) {
+        case Selector::TF:
+            if (_errorHistogramManager) {
+                std::stringstream cacheName;
+                ghoul::filesystem::File f = _filename;
+                cacheName << f.baseName() << "_" << nHistograms << "_errorHistograms";
+                std::string cacheFilename;
+                cacheFilename = FileSys.cacheManager()->cachedFilename(
+                    cacheName.str(), "", ghoul::filesystem::CacheManager::Persistent::Yes);
+                std::ifstream cacheFile(cacheFilename, std::ios::in | std::ios::binary);
+                if (cacheFile.is_open()) {
+                    // Read histograms from cache.
+                    cacheFile.close();
+                    LINFO("Loading histograms from " << cacheFilename);
+                    success &= _errorHistogramManager->loadFromFile(cacheFilename);
+                } else {
+                    // Build histograms from tsp file.
+                    LWARNING("Failed to open " << cacheFilename);
+                    if (success &= _errorHistogramManager->buildHistograms(nHistograms)) {
+                        LINFO("Writing cache to " << cacheFilename);
+                        _errorHistogramManager->saveToFile(cacheFilename);
+                    }
+                }
+                success &= _tfBrickSelector && _tfBrickSelector->initialize();
+            }
+            break;
+
+        case Selector::SIMPLE:
+            if (_histogramManager) {
+                std::stringstream cacheName;
+                ghoul::filesystem::File f = _filename;
+                cacheName << f.baseName() << "_" << nHistograms << "_histograms";
+                std::string cacheFilename;
+                cacheFilename = FileSys.cacheManager()->cachedFilename(
+                    cacheName.str(), "", ghoul::filesystem::CacheManager::Persistent::Yes);
+                std::ifstream cacheFile(cacheFilename, std::ios::in | std::ios::binary);
+                if (cacheFile.is_open()) {
+                    // Read histograms from cache.
+                    cacheFile.close();
+                    LINFO("Loading histograms from " << cacheFilename);
+                    success &= _histogramManager->loadFromFile(cacheFilename);
+                } else {
+                    // Build histograms from tsp file.
+                    LWARNING("Failed to open " << cacheFilename);
+                    if (success &= _histogramManager->buildHistograms(_tsp.get(), nHistograms)) {
+                        LINFO("Writing cache to " << cacheFilename);
+                        _histogramManager->saveToFile(cacheFilename);
+                    }
+                }
+                success &= _simpleTfBrickSelector && _simpleTfBrickSelector->initialize();
+            }
+            break;
+
+        case Selector::LOCAL:
+            if (_localErrorHistogramManager) {
+                std::stringstream cacheName;
+                ghoul::filesystem::File f = _filename;
+                cacheName << f.baseName() << "_" << nHistograms << "_localErrorHistograms";
+                std::string cacheFilename;
+                cacheFilename = FileSys.cacheManager()->cachedFilename(
+                    cacheName.str(), "", ghoul::filesystem::CacheManager::Persistent::Yes);
+                std::ifstream cacheFile(cacheFilename, std::ios::in | std::ios::binary);
+                if (cacheFile.is_open()) {
+                    // Read histograms from cache.
+                    cacheFile.close();
+                    LINFO("Loading histograms from " << cacheFilename);
+                    success &= _localErrorHistogramManager->loadFromFile(cacheFilename);
+                } else {
+                    // Build histograms from tsp file.
+                    LWARNING("Failed to open " << cacheFilename);
+                    if (success &= _localErrorHistogramManager->buildHistograms(nHistograms)) {
+                        LINFO("Writing cache to " << cacheFilename);
+                        _localErrorHistogramManager->saveToFile(cacheFilename);
+                    }
+                }
+                success &= _localTfBrickSelector && _localTfBrickSelector->initialize();
+            }
+            break;
+    }
+
+    return success;
+}
+/*
+void RenderableMultiresVolume::preResolve(ghoul::opengl::ProgramObject* program) {
+    RenderableVolume::preResolve(program);
+
+
+    std::stringstream ss;
+    ss << "opacity_" << getId();
+    program->setUniform(ss.str(), visible ? 1.0f : 0.0f);
+
+    ss.str(std::string());
+    ss << "stepSizeCoefficient_" << getId();
+    program->setUniform(ss.str(), _stepSizeCoefficient);
+
+    ss.str(std::string());
+    ss << "transferFunction_" << getId();
+    program->setUniform(ss.str(), getTextureUnit(_transferFunction->getTexture()));
+
+    ss.str(std::string());
+    ss << "textureAtlas_" << getId();
+    program->setUniform(ss.str(), getTextureUnit(_atlasManager->textureAtlas()));
+
+    ss.str(std::string());
+    ss << "atlasMapBlock_" << getId();
+    program->setSsboBinding(ss.str(), getSsboBinding(_atlasManager->atlasMapBuffer()));
+
+    ss.str(std::string());
+    ss << "gridType_" << getId();
+    program->setUniform(ss.str(), static_cast<int>(_tsp->header().gridType_));
+
+    ss.str(std::string());
+    ss << "maxNumBricksPerAxis_" << getId();
+    program->setUniform(ss.str(), static_cast<unsigned int>(_tsp->header().xNumBricks_));
+
+    ss.str(std::string());
+    ss << "paddedBrickDim_" << getId();
+    program->setUniform(ss.str(), static_cast<unsigned int>(_tsp->paddedBrickDim()));
+
+    ss.str(std::string());
+    ss << "atlasSize_" << getId();
+    glm::size3_t size = _atlasManager->textureSize();
+    glm::ivec3 atlasSize(size.x, size.y, size.z);
+    program->setUniform(ss.str(), atlasSize);
+
+    _timestep++;
+}
+*/
+/*
+std::vector<ghoul::opengl::Texture*> RenderableMultiresVolume::getTextures() {
+    std::vector<ghoul::opengl::Texture*> textures{_transferFunction->getTexture(), _atlasManager->textureAtlas()};
+    return textures;
+}
+
+std::vector<unsigned int> RenderableMultiresVolume::getBuffers() {
+    std::vector<unsigned int> buffers{_atlasManager->atlasMapBuffer()};
+    return buffers;
+}*/
+
+void RenderableMultiresVolume::update(const UpdateData& data) {
+        _timestep++;
+    _time = data.time;
+
+    if (_gatheringStats) {
+        std::chrono::system_clock::time_point frameEnd = std::chrono::system_clock::now();
+        std::chrono::duration<double> frameDuration = frameEnd - _frameStart;
+
+        // Make sure that the directory exists
+        ghoul::filesystem::File file(_statsFileName);
+        ghoul::filesystem::Directory directory(file.directoryName());
+        FileSys.createDirectory(directory, ghoul::filesystem::FileSystem::Recursive::Yes);
+
+        std::ofstream ofs(_statsFileName, std::ofstream::out);
+
+        ofs << frameDuration.count() << " "
+            << _selectionDuration.count() << " "
+            << _uploadDuration.count() << " "
+            << _nUsedBricks << " "
+            << _nStreamedBricks << " "
+            << _nDiskReads;
+
+        ofs.close();
+
+        _gatheringStats = false;
+    }
+    if (_statsToFile) {
+        // Start frame timer
+        _frameStart = std::chrono::system_clock::now();
+        _statsFileName = _statsToFileName;
+
+        _gatheringStats = true;
+        _statsToFile = false;
+    }
+
+
+    int numTimesteps = _tsp->header().numTimesteps_;
+    int currentTimestep;
+    bool visible = true;
+    if (_loop) {
+        currentTimestep = _timestep % numTimesteps;
+    }
+    else if (_useGlobalTime) {
+        double t = (_time - _startTime) / (_endTime - _startTime);
+        currentTimestep = t * numTimesteps;
+        visible = currentTimestep >= 0 && currentTimestep < numTimesteps;
+    }
+    else {
+        currentTimestep = _currentTime;
+    }
+
+    if (visible) {
+
+        std::chrono::system_clock::time_point selectionStart;
+        if (_gatheringStats) {
+            selectionStart = std::chrono::system_clock::now();
+        }
+
+        switch (_selector) {
+        case Selector::TF:
+            if (_tfBrickSelector) {
+                _tfBrickSelector->setMemoryBudget(_memoryBudget);
+                _tfBrickSelector->setStreamingBudget(_streamingBudget);
+                _tfBrickSelector->selectBricks(currentTimestep, _brickIndices);
+            }
+            break;
+        case Selector::SIMPLE:
+            if (_simpleTfBrickSelector) {
+                _simpleTfBrickSelector->setMemoryBudget(_memoryBudget);
+                _simpleTfBrickSelector->setStreamingBudget(_streamingBudget);
+                _simpleTfBrickSelector->selectBricks(currentTimestep, _brickIndices);
+            }
+            break;
+        case Selector::LOCAL:
+            if (_localTfBrickSelector) {
+                _localTfBrickSelector->setMemoryBudget(_memoryBudget);
+                _localTfBrickSelector->setStreamingBudget(_streamingBudget);
+                _localTfBrickSelector->selectBricks(currentTimestep, _brickIndices);
+            }
+            break;
+        }
+
+        std::chrono::system_clock::time_point uploadStart;
+        if (_gatheringStats) {
+            std::chrono::system_clock::time_point selectionEnd = std::chrono::system_clock::now();
+            _selectionDuration = selectionEnd - selectionStart;
+            uploadStart = selectionEnd;
+        }
+
+        _atlasManager->updateAtlas(AtlasManager::EVEN, _brickIndices);
+
+        if (_gatheringStats) {
+            std::chrono::system_clock::time_point uploadEnd = std::chrono::system_clock::now();
+            _uploadDuration = uploadEnd - uploadStart;
+            _nDiskReads = _atlasManager->getNumDiskReads();
+            _nUsedBricks = _atlasManager->getNumUsedBricks();
+            _nStreamedBricks = _atlasManager->getNumStreamedBricks();
+        }
+    }
+
+    if (_raycaster) {
+
+        glm::mat4 transform = glm::translate(glm::mat4(1.0), static_cast<glm::vec3>(_translation) * std::pow(10.0f, static_cast<float>(_scalingExponent)));
+        glm::vec3 eulerRotation = static_cast<glm::vec3>(_rotation);
+        transform = glm::rotate(transform, eulerRotation.x, glm::vec3(1, 0, 0));
+        transform = glm::rotate(transform, eulerRotation.y, glm::vec3(0, 1, 0));
+        transform = glm::rotate(transform, eulerRotation.z, glm::vec3(0, 0, 1));
+        transform = glm::scale(transform, static_cast<glm::vec3>(_scaling) * std::pow(10.0f, static_cast<float>(_scalingExponent)));
+
+    
+        _raycaster->setStepSizeCoefficient(_stepSizeCoefficient);
+        _raycaster->setModelTransform(transform);
+        //_raycaster->setTime(data.time);
+    }
+
+}
+
+void RenderableMultiresVolume::render(const RenderData& data, RendererTasks& tasks) {
+    RaycasterTask task{ _raycaster.get(), data };
+    tasks.raycasterTasks.push_back(task);
+}
+
+} // namespace openspace

modules/multiresvolume/rendering/renderablemultiresvolume.h

@@ -0,0 +1,159 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2015                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __RENDERABLEMULTIRESVOLUME_H__
+#define __RENDERABLEMULTIRESVOLUME_H__
+
+#include <vector>
+#include <chrono>
+#include <memory>
+#include <openspace/rendering/renderable.h>
+#include <openspace/rendering/transferfunction.h>
+#include <openspace/util/powerscaledcoordinate.h>
+#include <ghoul/misc/dictionary.h>
+#include <openspace/properties/scalarproperty.h>
+#include <openspace/properties/vectorproperty.h>
+#include <openspace/properties/stringproperty.h>
+#include <modules/multiresvolume/rendering/multiresvolumeraycaster.h>
+
+// Forward declare to minimize dependencies
+namespace ghoul {
+    namespace filesystem {
+        class File;
+    }
+    namespace opengl {
+        class ProgramObject;
+        class Texture;
+    }
+}
+
+namespace openspace {
+// Forward declare
+class TSP;
+class AtlasManager;
+class BrickSelector;
+class TfBrickSelector;
+class SimpleTfBrickSelector;
+class LocalTfBrickSelector;
+class HistogramManager;
+class ErrorHistogramManager;
+class LocalErrorHistogramManager;
+
+class RenderableMultiresVolume : public Renderable {
+
+
+
+
+public:
+    RenderableMultiresVolume(const ghoul::Dictionary& dictionary);
+    ~RenderableMultiresVolume();
+
+    enum Selector {TF, SIMPLE, LOCAL};
+
+    bool setSelectorType(Selector selector);
+    bool initializeSelector();
+
+    bool initialize() override;
+    bool deinitialize() override;
+
+    bool isReady() const override;
+
+    virtual void update(const UpdateData& data) override;
+    virtual void RenderableMultiresVolume::render(const RenderData& data, RendererTasks& tasks);
+    //virtual void preResolve(ghoul::opengl::ProgramObject* program) override;
+    //virtual std::string getHeaderPath() override;
+    //virtual std::string getHelperPath() override;
+    //virtual std::vector<ghoul::opengl::Texture*> getTextures() override;
+    //virtual std::vector<unsigned int> getBuffers() override;
+
+private:
+    double _time;
+    double _startTime;
+    double _endTime;
+
+    properties::BoolProperty _useGlobalTime;
+    properties::BoolProperty _loop;
+    properties::IntProperty _currentTime; // used to vary time, if not using global time nor looping
+    properties::IntProperty _memoryBudget;
+    properties::IntProperty _streamingBudget;
+    properties::FloatProperty _stepSizeCoefficient;
+    properties::StringProperty _selectorName;
+    properties::BoolProperty _statsToFile;
+    properties::StringProperty _statsToFileName;
+
+    // Stats timers
+    std::string _statsFileName;
+    bool _gatheringStats;
+    std::chrono::system_clock::time_point _frameStart;
+    std::chrono::duration<double> _selectionDuration;
+    std::chrono::duration<double> _uploadDuration;
+    unsigned int _nDiskReads;
+    unsigned int _nUsedBricks;
+    unsigned int _nStreamedBricks;
+
+    int _timestep;
+
+    std::string _filename;
+
+    std::string _transferFunctionName;
+    std::string _volumeName;
+
+    std::string _transferFunctionPath;
+
+    std::shared_ptr<TransferFunction> _transferFunction;
+
+    float _spatialTolerance;
+    float _temporalTolerance;
+
+    std::shared_ptr<TSP> _tsp;
+    std::vector<int> _brickIndices;
+    int _atlasMapSize;
+
+    std::shared_ptr<AtlasManager> _atlasManager;
+
+    std::unique_ptr<MultiresVolumeRaycaster> _raycaster;
+
+    TfBrickSelector* _tfBrickSelector;
+    SimpleTfBrickSelector* _simpleTfBrickSelector;
+    LocalTfBrickSelector* _localTfBrickSelector;
+
+    Selector _selector;
+
+    HistogramManager* _histogramManager;
+    ErrorHistogramManager* _errorHistogramManager;
+    LocalErrorHistogramManager* _localErrorHistogramManager;
+
+    float _w;
+    PowerScaledCoordinate _pscOffset;
+
+    properties::IntProperty _scalingExponent;
+    properties::Vec3Property _translation;
+    properties::Vec3Property _rotation;
+    properties::Vec3Property _scaling;
+
+};
+
+} // namespace openspace
+
+#endif // __RENDERABLEMULTIRESVOLUME_H__

modules/multiresvolume/rendering/shenbrickselector.cpp

@@ -0,0 +1,134 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/shenbrickselector.h>
+
+namespace {
+    const std::string _loggerCat = "ShenBrickSelector";
+}
+
+namespace openspace {
+
+ShenBrickSelector::ShenBrickSelector(TSP* tsp, float spatialTolerance, float temporalTolerance)
+    : _tsp(tsp)
+    , _spatialTolerance(spatialTolerance)
+    , _temporalTolerance(temporalTolerance) {}
+
+ShenBrickSelector::~ShenBrickSelector() {
+
+}
+
+void ShenBrickSelector::setSpatialTolerance(float spatialTolerance) {
+    _spatialTolerance = spatialTolerance;
+}
+
+void ShenBrickSelector::setTemporalTolerance(float temporalTolerance) {
+    _temporalTolerance = temporalTolerance;
+}
+
+void ShenBrickSelector::selectBricks(int timestep, std::vector<int>& bricks) {
+    int numTimeSteps = _tsp->header().numTimesteps_;
+    BrickCover coveredBricks(_tsp->header().xNumBricks_);
+    selectBricks(timestep, 0, 0, 0, numTimeSteps, coveredBricks, bricks);
+}
+
+/**
+ * Traverse the Octree in the BST root
+ */
+void ShenBrickSelector::traverseOT(int timestep, unsigned int brickIndex, BrickCover coveredBricks, std::vector<int>& bricks) {
+    unsigned int firstChild = _tsp->getFirstOctreeChild(brickIndex);
+    int numTimeSteps = _tsp->header().numTimesteps_;
+    for (unsigned int i = 0; i < 8; i++) {
+        unsigned int child = firstChild + i;
+        BrickCover cover = coveredBricks.split(i % 2, (i/2) % 2, (i/4));
+        selectBricks(timestep, child, child, 0, numTimeSteps, cover, bricks);
+    }
+}
+
+void ShenBrickSelector::traverseBST(int timestep,
+                                unsigned int brickIndex,
+                                unsigned int bstRootBrickIndex,
+                                int timeSpanStart,
+                                int timeSpanEnd,
+                                BrickCover coveredBricks,
+                                std::vector<int>& bricks) {
+
+
+    int timeSpanCenter = timeSpanStart + (timeSpanEnd - timeSpanStart) / 2;
+    unsigned int bstChild;
+    if (timestep <= timeSpanCenter) {
+        bstChild = _tsp->getBstLeft(brickIndex);
+        timeSpanEnd = timeSpanCenter;
+    } else {
+        bstChild = _tsp->getBstRight(brickIndex);
+        timeSpanStart = timeSpanCenter;
+    }
+    selectBricks(timestep, bstChild, bstRootBrickIndex, timeSpanStart, timeSpanEnd, coveredBricks, bricks);
+}
+
+void ShenBrickSelector::selectBricks(int timestep,
+                                 unsigned int brickIndex,
+                                 unsigned int bstRootBrickIndex,
+                                 int timeSpanStart,
+                                 int timeSpanEnd,
+                                 BrickCover coveredBricks,
+                                 std::vector<int>& bricks) {
+
+    if (_tsp->getTemporalError(brickIndex) <= _temporalTolerance) {
+        if (_tsp->isOctreeLeaf(bstRootBrickIndex)) {
+            selectCover(coveredBricks, brickIndex, bricks);
+        } else if (_tsp->getSpatialError(brickIndex) <= _spatialTolerance) {
+            selectCover(coveredBricks, brickIndex, bricks);
+        } else if (_tsp->isBstLeaf(brickIndex)) {
+            traverseOT(timestep, bstRootBrickIndex, coveredBricks, bricks);
+        } else {
+            traverseBST(timestep, brickIndex, bstRootBrickIndex, timeSpanStart, timeSpanEnd, coveredBricks, bricks);
+        }
+    } else if (_tsp->isBstLeaf(brickIndex)) {
+        if (_tsp->isOctreeLeaf(bstRootBrickIndex)) {
+            selectCover(coveredBricks, brickIndex, bricks);
+        } else {
+            traverseOT(timestep, bstRootBrickIndex, coveredBricks, bricks);
+        }
+    } else {
+        traverseBST(timestep, brickIndex, bstRootBrickIndex, timeSpanStart, timeSpanEnd, coveredBricks, bricks);
+    }
+}
+
+int ShenBrickSelector::linearCoords(int x, int y, int z) {
+    const TSP::Header &header = _tsp->header();
+    return x + (header.xNumBricks_ * y) + (header.xNumBricks_ * header.yNumBricks_ * z);
+}
+
+void ShenBrickSelector::selectCover(BrickCover coveredBricks, unsigned int brickIndex, std::vector<int>& bricks) {
+    for (int z = coveredBricks.lowZ; z < coveredBricks.highZ; z++) {
+        for (int y = coveredBricks.lowY; y < coveredBricks.highY; y++) {
+            for (int x = coveredBricks.lowX; x < coveredBricks.highX; x++) {
+                bricks[linearCoords(x, y, z)] = brickIndex;
+            }
+        }
+    }
+}
+
+} // namespace openspace

modules/multiresvolume/rendering/shenbrickselector.h

@@ -0,0 +1,77 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __SHENBRICKSELECTOR_H__
+#define __SHENBRICKSELECTOR_H__
+
+#include <vector>
+#include <modules/multiresvolume/rendering/brickselector.h>
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/brickcover.h>
+
+namespace openspace {
+
+
+class ShenBrickSelector : public BrickSelector {
+public:
+    ShenBrickSelector(TSP* tsp, float spatialTolerance, float temporalTolerance);
+    ~ShenBrickSelector();
+    void setSpatialTolerance(float spatialTolerance);
+    void setTemporalTolerance(float temporalTolerance);
+    void selectBricks(int timestep,
+                      std::vector<int>& bricks);
+private:
+    TSP* _tsp;
+    float _spatialTolerance;
+    float _temporalTolerance;
+
+    void traverseOT(int timestep,
+                    unsigned int brickIndex,
+                    BrickCover coveredBricks,
+                    std::vector<int>& bricks);
+
+    void traverseBST(int timestep,
+                     unsigned int brickIndex,
+                     unsigned int bstRootBrickIndex,
+                     int timeSpanStart,
+                     int timeSpanEnd,
+                     BrickCover coveredBricks,
+                     std::vector<int>& bricks);
+
+    void selectBricks(int timestep,
+                      unsigned int brickIndex,
+                      unsigned int bstRootBrickIndex,
+                      int timeSpanStart,
+                      int timeSpanEnd,
+                      BrickCover coveredBricks,
+                      std::vector<int>& bricks);
+
+    int linearCoords(int x, int y, int z);
+
+    void selectCover(BrickCover coveredBricks, unsigned int brickIndex, std::vector<int>& bricks);
+};
+
+} // namespace openspace
+
+#endif // __SHENBRICKSELECTOR_H__

modules/multiresvolume/rendering/simpletfbrickselector.cpp

@@ -0,0 +1,337 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/simpletfbrickselector.h>
+#include <modules/multiresvolume/rendering/histogrammanager.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+#include <openspace/rendering/transferfunction.h>
+#include <algorithm>
+#include <cassert>
+#include <ghoul/logging/logmanager.h>
+
+namespace {
+    const std::string _loggerCat = "SimpleTfBrickSelector";
+}
+
+namespace openspace {
+
+    SimpleTfBrickSelector::SimpleTfBrickSelector(TSP* tsp, HistogramManager* hm, TransferFunction* tf, int memoryBudget, int streamingBudget)
+    : _tsp(tsp)
+    , _histogramManager(hm)
+    , _transferFunction(tf)
+    , _memoryBudget(memoryBudget)
+    , _streamingBudget(streamingBudget) {}
+
+SimpleTfBrickSelector::~SimpleTfBrickSelector() {}
+
+bool SimpleTfBrickSelector::initialize() {
+    return true;
+}
+
+void SimpleTfBrickSelector::setMemoryBudget(int memoryBudget) {
+    _memoryBudget = memoryBudget;
+}
+
+void SimpleTfBrickSelector::setStreamingBudget(int streamingBudget) {
+    _streamingBudget = streamingBudget;
+}
+
+void SimpleTfBrickSelector::selectBricks(int timestep, std::vector<int>& bricks) {
+    int numTimeSteps = _tsp->header().numTimesteps_;
+    int numBricksPerDim = _tsp->header().xNumBricks_;
+
+    unsigned int rootNode = 0;
+    BrickSelection::SplitType splitType;
+    float rootSplitPoints = splitPoints(rootNode, splitType);
+
+    BrickSelection brickSelection = BrickSelection(numBricksPerDim, numTimeSteps, splitType, rootSplitPoints);
+
+    std::vector<BrickSelection> priorityQueue;
+    std::vector<BrickSelection> leafSelections;
+    std::vector<BrickSelection> temporalSplitQueue;
+    std::vector<BrickSelection> deadEnds;
+
+    if (splitType != BrickSelection::SplitType::None) {
+        priorityQueue.push_back(brickSelection);
+    } else {
+        leafSelections.push_back(brickSelection);
+    }
+
+    int memoryBudget = _memoryBudget;
+    int totalStreamingBudget = _streamingBudget * numTimeSteps;
+    int nBricksInMemory = 1;
+    int nStreamedBricks = 1;
+
+    while (nBricksInMemory <= memoryBudget - 7 && priorityQueue.size() > 0) {
+        std::pop_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+        BrickSelection bs = priorityQueue.back();
+
+        // TODO: handle edge case when we can only afford temporal splits or no split (only 1 spot left)
+
+        unsigned int brickIndex = bs.brickIndex;
+        priorityQueue.pop_back();
+        if (bs.splitType == BrickSelection::SplitType::Temporal) {
+            unsigned int childBrickIndex;
+            bool pickRightTimeChild = bs.timestepInRightChild(timestep);
+
+            // On average on the whole time period, splitting this spatial brick in two time steps
+            // would generate twice as much streaming. Current number of streams of this spatial brick
+            // is 2^nTemporalSplits over the whole time period.
+            int newStreams = std::pow(2, bs.nTemporalSplits);
+
+            // Refining this one more step would require the double amount of streams
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // Reached dead end (streaming budget would be exceeded)
+                deadEnds.push_back(bs);
+                break;
+            }
+            nStreamedBricks += newStreams;
+
+            if (pickRightTimeChild) {
+                childBrickIndex = _tsp->getBstRight(brickIndex);
+            } else {
+                childBrickIndex = _tsp->getBstLeft(brickIndex);
+            }
+
+            BrickSelection::SplitType childSplitType;
+            float childSplitPoints = splitPoints(childBrickIndex, childSplitType);
+            BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, childSplitType, childSplitPoints);
+
+            if (childSplitType != BrickSelection::SplitType::None) {
+                priorityQueue.push_back(childSelection);
+                std::push_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                leafSelections.push_back(childSelection);
+            }
+        } else if (bs.splitType == BrickSelection::SplitType::Spatial) {
+            nBricksInMemory += 7; // Remove one and add eight.
+            unsigned int firstChild = _tsp->getFirstOctreeChild(brickIndex);
+
+            // On average on the whole time period, splitting this spatial brick into eight spatial bricks
+            // would generate eight times as much streaming. Current number of streams of this spatial brick
+            // is 2^nTemporalStreams over the whole time period.
+            int newStreams = 7*std::pow(2, bs.nTemporalSplits);
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // Reached dead end (streaming budget would be exceeded)
+                // However, temporal split might be possible
+                if (bs.splitType != BrickSelection::SplitType::Temporal) {
+                    bs.splitType = BrickSelection::SplitType::Temporal;
+                    bs.splitPoints = temporalSplitPoints(bs.brickIndex);
+                }
+                if (bs.splitPoints > -1) {
+                    temporalSplitQueue.push_back(bs);
+                } else {
+                    deadEnds.push_back(bs);
+                }
+                break;
+            }
+            nStreamedBricks += newStreams;
+
+            for (unsigned int i = 0; i < 8; i++) {
+                unsigned int childBrickIndex = firstChild + i;
+
+                BrickSelection::SplitType childSplitType;
+                float childSplitPoints = splitPoints(childBrickIndex, childSplitType);
+                BrickSelection childSelection = bs.splitSpatially(i % 2, (i/2) % 2, i/4, childBrickIndex, childSplitType, childSplitPoints);
+
+                if (childSplitType != BrickSelection::SplitType::None) {
+                    priorityQueue.push_back(childSelection);
+                    std::push_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+                } else {
+                    leafSelections.push_back(childSelection);
+                }
+            }
+        }
+    }
+
+    if (nStreamedBricks < totalStreamingBudget) {
+        while (priorityQueue.size() > 0) {
+            BrickSelection bs = priorityQueue.back();
+            if (bs.splitType != BrickSelection::SplitType::Temporal) {
+                bs.splitType = BrickSelection::SplitType::Temporal;
+                bs.splitPoints = temporalSplitPoints(bs.brickIndex);
+            }
+            priorityQueue.pop_back();
+            if (bs.splitPoints > -1) {
+                temporalSplitQueue.push_back(bs);
+                std::push_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                deadEnds.push_back(bs);
+            }
+        }
+
+        // Keep splitting until it's not possible anymore
+        while (nStreamedBricks < totalStreamingBudget - 1 && temporalSplitQueue.size() > 0) {
+            std::pop_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            BrickSelection bs = temporalSplitQueue.back();
+            temporalSplitQueue.pop_back();
+
+            unsigned int brickIndex = bs.brickIndex;
+            int newStreams = std::pow(2, bs.nTemporalSplits);
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // The current best choice would make us exceed the streaming budget, try next instead.
+                deadEnds.push_back(bs);
+                continue;
+            }
+
+            nStreamedBricks += newStreams;
+            unsigned int childBrickIndex;
+            bool pickRightTimeChild = bs.timestepInRightChild(timestep);
+
+            if (pickRightTimeChild) {
+                childBrickIndex = _tsp->getBstRight(brickIndex);
+            } else {
+                childBrickIndex = _tsp->getBstLeft(brickIndex);
+            }
+
+            float childSplitPoints = temporalSplitPoints(childBrickIndex);
+
+            if (childSplitPoints > -1) {
+                BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, BrickSelection::SplitType::Temporal, childSplitPoints);
+                temporalSplitQueue.push_back(childSelection);
+                std::push_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, BrickSelection::SplitType::None, -1);
+                deadEnds.push_back(childSelection);
+            }
+        }
+    } else {
+        // Write selected inner nodes to brickSelection vector
+        for (const BrickSelection& bs : priorityQueue) {
+            writeSelection(bs, bricks);
+        } 
+    }
+
+    // Write selected inner nodes to brickSelection vector
+    for (const BrickSelection& bs : temporalSplitQueue) {
+        writeSelection(bs, bricks);
+    }
+    for (const BrickSelection& bs : deadEnds) {
+        writeSelection(bs, bricks);
+    }
+    // Write selected leaf nodes to brickSelection vector
+    for (const BrickSelection& bs : leafSelections) {
+        writeSelection(bs, bricks);
+    }
+
+    //std::cout << "Bricks in memory: " << nBricksInMemory << "/" << _memoryBudget << "___\t\t"
+    //          << "Streamed bricks:  " << nStreamedBricks << "/" << totalStreamingBudget << std::flush << "___\r";
+}
+
+float SimpleTfBrickSelector::temporalSplitPoints(unsigned int brickIndex) {
+    if (_tsp->isBstLeaf(brickIndex)) {
+        return -1;
+    }
+    return _brickImportances[brickIndex] * 0.5;
+}
+
+float SimpleTfBrickSelector::spatialSplitPoints(unsigned int brickIndex) {
+    if (_tsp->isOctreeLeaf(brickIndex)) {
+        return -1;
+    }
+    return _brickImportances[brickIndex] * 0.125;
+}
+
+float SimpleTfBrickSelector::splitPoints(unsigned int brickIndex, BrickSelection::SplitType& splitType) {
+    float temporalPoints = temporalSplitPoints(brickIndex);
+    float spatialPoints = spatialSplitPoints(brickIndex);
+    float splitPoints;
+
+    if (spatialPoints > 0 && spatialPoints > temporalPoints) {
+        splitPoints = spatialPoints;
+        splitType = BrickSelection::SplitType::Spatial;
+    } else if (temporalPoints > 0) {
+        splitPoints = temporalPoints;
+        splitType = BrickSelection::SplitType::Temporal;
+    } else {
+        splitPoints = -1;
+        splitType = BrickSelection::SplitType::None;
+    }
+
+    return splitPoints;
+}
+
+
+bool SimpleTfBrickSelector::calculateBrickImportances() {
+    TransferFunction *tf = _transferFunction;
+    if (!tf) return false;
+
+    float tfWidth = tf->width();
+    if (tfWidth <= 0) return false;
+
+    /*    std::vector<float> gradients(tfWidth - 1);
+    for (size_t offset = 0; offset < tfWidth - 1; offset++) {
+        glm::vec4 prevRgba = tf->sample(offset);
+        glm::vec4 nextRgba = tf->sample(offset + 1);
+
+        float colorDifference = glm::distance(prevRgba, nextRgba);
+        float alpha = (prevRgba.w + nextRgba.w) * 0.5;
+
+        gradients[offset] = colorDifference*alpha;
+        }*/
+
+    unsigned int nHistograms = _tsp->numTotalNodes();
+    _brickImportances = std::vector<float>(nHistograms);
+
+    for (unsigned int brickIndex = 0; brickIndex < nHistograms; brickIndex++) {
+        const Histogram* histogram = _histogramManager->getHistogram(brickIndex);
+        if (!histogram->isValid()) {
+            return false;
+        }
+
+        float dotProduct = 0;
+        for (int i = 0; i < tf->width(); i++) {
+            float x = float(i) / tfWidth;
+            float sample = histogram->interpolate(x);
+
+            assert(sample >= 0);
+            dotProduct += sample * tf->sample(i).w;
+        }
+        _brickImportances[brickIndex] = dotProduct;
+    }
+
+    LINFO("Updated brick importances");
+
+    return true;
+}
+
+int SimpleTfBrickSelector::linearCoords(int x, int y, int z) {
+    const TSP::Header &header = _tsp->header();
+    return x + (header.xNumBricks_ * y) + (header.xNumBricks_ * header.yNumBricks_ * z);
+}
+
+void SimpleTfBrickSelector::writeSelection(BrickSelection brickSelection, std::vector<int>& bricks) {
+    BrickCover coveredBricks = brickSelection.cover;
+    for (int z = coveredBricks.lowZ; z < coveredBricks.highZ; z++) {
+        for (int y = coveredBricks.lowY; y < coveredBricks.highY; y++) {
+            for (int x = coveredBricks.lowX; x < coveredBricks.highX; x++) {
+                bricks[linearCoords(x, y, z)] = brickSelection.brickIndex;
+            }
+        }
+    }
+}
+
+
+} // namespace openspace

modules/multiresvolume/rendering/simpletfbrickselector.h

@@ -0,0 +1,71 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __SIMPLETFBRICKSELECTOR_H__
+#define __SIMPLETFBRICKSELECTOR_H__
+
+#include <vector>
+#include <modules/multiresvolume/rendering/brickselection.h>
+#include <modules/multiresvolume/rendering/brickselector.h>
+#include <modules/multiresvolume/rendering/brickcover.h>
+
+
+namespace openspace {
+
+class TSP;
+class HistogramManager;
+class TransferFunction;
+
+class SimpleTfBrickSelector : public BrickSelector {
+public:
+    SimpleTfBrickSelector(TSP* tsp, HistogramManager* hm, TransferFunction* tf, int memoryBudget, int streamingBudget);
+    ~SimpleTfBrickSelector();
+
+    virtual bool initialize();
+
+    void selectBricks(int timestep, std::vector<int>& bricks);
+    void setMemoryBudget(int memoryBudget);
+    void setStreamingBudget(int streamingBudget);
+    bool calculateBrickImportances();
+ private:
+
+    TSP* _tsp;
+    HistogramManager* _histogramManager;
+    TransferFunction* _transferFunction;
+    std::vector<float> _brickImportances;
+    float spatialSplitPoints(unsigned int brickIndex);
+    float temporalSplitPoints(unsigned int brickIndex);
+    float splitPoints(unsigned int brickIndex, BrickSelection::SplitType& splitType);
+
+    int linearCoords(int x, int y, int z);
+    void writeSelection(BrickSelection coveredBricks, std::vector<int>& bricks);
+
+    int _memoryBudget;
+    int _streamingBudget;
+};
+
+} // namespace openspace
+
+#endif // __SIMPLETFBRICKSELECTOR_H__
+

modules/multiresvolume/rendering/tfbrickselector.cpp

@@ -0,0 +1,395 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <modules/multiresvolume/rendering/tsp.h>
+#include <modules/multiresvolume/rendering/tfbrickselector.h>
+#include <modules/multiresvolume/rendering/errorhistogrammanager.h>
+#include <modules/multiresvolume/rendering/histogram.h>
+#include <openspace/rendering/transferfunction.h>
+#include <algorithm>
+#include <cassert>
+
+namespace {
+    const std::string _loggerCat = "TfBrickSelector";
+}
+
+namespace openspace {
+
+TfBrickSelector::TfBrickSelector(TSP* tsp, ErrorHistogramManager* hm, TransferFunction* tf, int memoryBudget, int streamingBudget)
+    : _tsp(tsp)
+    , _histogramManager(hm)
+    , _transferFunction(tf)
+    , _memoryBudget(memoryBudget)
+    , _streamingBudget(streamingBudget) {}
+
+TfBrickSelector::~TfBrickSelector() {}
+
+bool TfBrickSelector::initialize() {
+    return calculateBrickErrors();
+}
+
+void TfBrickSelector::setMemoryBudget(int memoryBudget) {
+    _memoryBudget = memoryBudget;
+}
+
+void TfBrickSelector::setStreamingBudget(int streamingBudget) {
+    _streamingBudget = streamingBudget;
+}
+
+void TfBrickSelector::selectBricks(int timestep, std::vector<int>& bricks) {
+    int numTimeSteps = _tsp->header().numTimesteps_;
+    int numBricksPerDim = _tsp->header().xNumBricks_;
+
+    unsigned int rootNode = 0;
+    BrickSelection::SplitType splitType;
+    float rootSplitPoints = splitPoints(rootNode, splitType);
+
+
+    BrickSelection brickSelection = BrickSelection(numBricksPerDim, numTimeSteps, splitType, rootSplitPoints);
+
+    std::vector<BrickSelection> priorityQueue;
+    std::vector<BrickSelection> leafSelections;
+    std::vector<BrickSelection> temporalSplitQueue;
+    std::vector<BrickSelection> deadEnds;
+
+    if (splitType != BrickSelection::SplitType::None) {
+        priorityQueue.push_back(brickSelection);
+    } else {
+        leafSelections.push_back(brickSelection);
+    }
+
+    int memoryBudget = _memoryBudget;
+    int totalStreamingBudget = _streamingBudget * numTimeSteps;
+    int nBricksInMemory = 1;
+    int nStreamedBricks = 1;
+
+    // First loop: While neither the memory nor the streaming budget is reached,
+    // try to optimize for visual quality vs memory.
+
+    while (nBricksInMemory <= memoryBudget - 7 && priorityQueue.size() > 0) {
+
+        std::pop_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+
+        BrickSelection bs = priorityQueue.back();
+        unsigned int brickIndex = bs.brickIndex;
+        priorityQueue.pop_back();
+
+        // TODO: handle edge case when we can only afford temporal splits or no split (only 1 spot left)
+        if (bs.splitType == BrickSelection::SplitType::Temporal) {
+            unsigned int childBrickIndex;
+            bool pickRightTimeChild = bs.timestepInRightChild(timestep);
+	    //	    assert(!pickRightTimeChild && "picked right child");
+
+            // On average on the whole time period, splitting this spatial brick in two time steps
+            // would generate twice as much streaming. Current number of streams of this spatial brick
+            // is 2^nTemporalSplits over the whole time period.
+
+            int newStreams = std::pow(2, bs.nTemporalSplits);
+            //std::cout << "preparing for " << newStreams << " new streams" << std::endl;
+
+            // Refining this one more step would require the double amount of streams
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                //std::cout << "Reached streaming budget when splitting temporally! Breaking" << std::endl;
+                // Reached dead end (streaming budget would be exceeded)
+                deadEnds.push_back(bs);
+                break;
+            }
+            nStreamedBricks += newStreams;
+
+            if (pickRightTimeChild) {
+                childBrickIndex = _tsp->getBstRight(brickIndex);
+            } else {
+                childBrickIndex = _tsp->getBstLeft(brickIndex);
+            }
+
+            BrickSelection::SplitType childSplitType;
+            float childSplitPoints = splitPoints(childBrickIndex, childSplitType);
+            BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, childSplitType, childSplitPoints);
+            if (childSplitType != BrickSelection::SplitType::None) {
+                priorityQueue.push_back(childSelection);
+                std::push_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                leafSelections.push_back(childSelection);
+            }
+        } else if (bs.splitType == BrickSelection::SplitType::Spatial) {
+            nBricksInMemory += 7; // Remove one and add eight.
+            unsigned int firstChild = _tsp->getFirstOctreeChild(brickIndex);
+
+            // On average on the whole time period, splitting this spatial brick into eight spatial bricks
+            // would generate eight times as much streaming. Current number of streams of this spatial brick
+            // is 2^nTemporalStreams over the whole time period.
+            int newStreams = 7*std::pow(2, bs.nTemporalSplits);
+
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // Reached dead end (streaming budget would be exceeded)
+                // However, temporal split might be possible
+                //std::cout << "Reached streaming budget when splitting spatially! Breaking" << std::endl;
+                if (bs.splitType != BrickSelection::SplitType::Temporal) {
+                    bs.splitType = BrickSelection::SplitType::Temporal;
+                    bs.splitPoints = temporalSplitPoints(bs.brickIndex);
+                }
+                if (bs.splitPoints > -1) {
+                    temporalSplitQueue.push_back(bs);
+                } else {
+                    deadEnds.push_back(bs);
+                }
+                break;
+            }
+            nStreamedBricks += newStreams;
+
+            for (unsigned int i = 0; i < 8; i++) {
+                unsigned int childBrickIndex = firstChild + i;
+
+                BrickSelection::SplitType childSplitType;
+                float childSplitPoints = splitPoints(childBrickIndex, childSplitType);
+                //std::cout << "Splitting spatially." << std::endl;
+                BrickSelection childSelection = bs.splitSpatially(i % 2, (i/2) % 2, i/4, childBrickIndex, childSplitType, childSplitPoints);
+
+                if (childSplitType != BrickSelection::SplitType::None) {
+                    priorityQueue.push_back(childSelection);
+                    std::push_heap(priorityQueue.begin(), priorityQueue.end(), BrickSelection::compareSplitPoints);
+                } else {
+                    leafSelections.push_back(childSelection);
+                }
+            }
+        }
+    }
+
+    if (nBricksInMemory <= memoryBudget - 7) {
+        //std::cout << "memory budget not reached. " << nBricksInMemory << " out of " << memoryBudget << std::endl;
+    }
+
+    // Is it possible that we may stream more bricks?
+    if (nStreamedBricks < totalStreamingBudget - 1) {
+        //std::cout << "streaming budget not reached. " << nStreamedBricks << " out of " << totalStreamingBudget << std::endl;
+        //std::cout << "there are " << priorityQueue.size() << " elements left in priority queue." << std::endl;
+
+        while (priorityQueue.size() > 0) {
+            BrickSelection bs = priorityQueue.back();
+            if (bs.splitType != BrickSelection::SplitType::Temporal) {
+                bs.splitType = BrickSelection::SplitType::Temporal;
+                bs.splitPoints = temporalSplitPoints(bs.brickIndex);
+            }
+            priorityQueue.pop_back();
+            if (bs.splitPoints > -1) {
+                temporalSplitQueue.push_back(bs);
+                std::push_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                deadEnds.push_back(bs);
+            }
+        }
+
+        // Keep splitting until it's not possible anymore
+        while (nStreamedBricks < totalStreamingBudget - 1 && temporalSplitQueue.size() > 0) {
+            std::pop_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            BrickSelection bs = temporalSplitQueue.back();
+            temporalSplitQueue.pop_back();
+
+            unsigned int brickIndex = bs.brickIndex;
+            int newStreams = std::pow(2, bs.nTemporalSplits);
+            if (nStreamedBricks + newStreams > totalStreamingBudget) {
+                // The current best choice would make us exceed the streaming budget, try next instead.
+                deadEnds.push_back(bs);
+                //std::cout << "Dead end trying to split " << brickIndex << ". Streamed would be " << (nStreamedBricks + newStreams) << std::endl;
+                continue;
+            }
+
+            nStreamedBricks += newStreams;
+            unsigned int childBrickIndex;
+            bool pickRightTimeChild = bs.timestepInRightChild(timestep);
+
+            if (pickRightTimeChild) {
+                childBrickIndex = _tsp->getBstRight(brickIndex);
+            } else {
+                childBrickIndex = _tsp->getBstLeft(brickIndex);
+            }
+
+            float childSplitPoints = temporalSplitPoints(childBrickIndex);
+
+            if (childSplitPoints > -1) {
+                BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, BrickSelection::SplitType::Temporal, childSplitPoints);
+                temporalSplitQueue.push_back(childSelection);
+                std::push_heap(temporalSplitQueue.begin(), temporalSplitQueue.end(), BrickSelection::compareSplitPoints);
+            } else {
+                BrickSelection childSelection = bs.splitTemporally(pickRightTimeChild, childBrickIndex, BrickSelection::SplitType::None, -1);
+                deadEnds.push_back(childSelection);
+            }
+        }
+    } else {
+        // Write selected inner nodes to brickSelection vector
+        //std::cout << "priority queue: " << priorityQueue.size() << std::endl;
+        for (const BrickSelection& bs : priorityQueue) {
+            writeSelection(bs, bricks);
+        }
+    }
+
+    //std::cout << "temporal split queue: " << temporalSplitQueue.size() << std::endl;
+    // Write selected inner nodes to brickSelection vector
+    for (const BrickSelection& bs : temporalSplitQueue) {
+        writeSelection(bs, bricks);
+    }
+    //std::cout << "dead ends: " << deadEnds.size() << std::endl;
+    for (const BrickSelection& bs : deadEnds) {
+        writeSelection(bs, bricks);
+    }
+    // Write selected leaf nodes to brickSelection vector
+    //std::cout << "leaf selections: " << leafSelections.size() << std::endl;
+    for (const BrickSelection& bs : leafSelections) {
+        writeSelection(bs, bricks);
+    }
+
+    //std::cout << "Bricks in memory: " << nBricksInMemory << "/" << _memoryBudget << "___\t\t"
+    //          << "Streamed bricks:  " << nStreamedBricks << "/" << totalStreamingBudget << std::flush << "___\r";
+}
+
+float TfBrickSelector::temporalSplitPoints(unsigned int brickIndex) {
+    if (_tsp->isBstLeaf(brickIndex)) {
+        return -1;
+    }
+
+    unsigned int leftChild = _tsp->getBstLeft(brickIndex);
+    unsigned int rightChild = _tsp->getBstRight(brickIndex);
+
+    float currentError = _brickErrors[brickIndex];
+    float splitError = _brickErrors[leftChild] + _brickErrors[rightChild];
+
+    /*if (currentError + 0.001 < splitError) {
+        std::cout << "Warning! (TEMPORAL SPLIT) Current error " << currentError << " is smaller than split error " << splitError << "." << std::endl;
+	}*/
+    float diff = currentError - splitError;
+    if (diff < 0.0) {
+	//std::cout << "local temporal split minimum for brick " << brickIndex << std::endl;
+	diff = -diff;
+    }
+    return diff * 0.5;
+}
+
+float TfBrickSelector::spatialSplitPoints(unsigned int brickIndex) {
+    if (_tsp->isOctreeLeaf(brickIndex)) {
+        return -1;
+    }
+
+
+    float currentError = _brickErrors[brickIndex];
+    float splitError = 0;
+
+    unsigned int firstChild = _tsp->getFirstOctreeChild(brickIndex);
+    for (unsigned int i = 0; i < 8; i++) {
+        unsigned int child = firstChild + i;
+        splitError += _brickErrors[child];
+    }
+
+    /*if (currentError + 0.001 < splitError) {
+        std::cout << "Warning! (SPATIAL SPLIT) Current error " << currentError << " is smaller than split error " << splitError << "." << std::endl;
+	}*/
+
+    float diff = currentError - splitError;
+    if (diff < 0.0) {
+	//std::cout << "local spatial split minimum for brick " << brickIndex << std::endl;
+	diff = -diff;
+    }
+
+    return diff * 0.125;
+}
+
+float TfBrickSelector::splitPoints(unsigned int brickIndex, BrickSelection::SplitType& splitType) {
+    float temporalPoints = temporalSplitPoints(brickIndex);
+    float spatialPoints = spatialSplitPoints(brickIndex);
+    float splitPoints;
+
+    if (spatialPoints > 0 && spatialPoints > temporalPoints) {
+        splitPoints = spatialPoints;
+        splitType = BrickSelection::SplitType::Spatial;
+    } else if (temporalPoints > 0) {
+        splitPoints = temporalPoints;
+        splitType = BrickSelection::SplitType::Temporal;
+    } else {
+        splitPoints = -1;
+        splitType = BrickSelection::SplitType::None;
+    }
+    return splitPoints;
+}
+
+
+bool TfBrickSelector::calculateBrickErrors() {
+    TransferFunction *tf = _transferFunction;
+    if (!tf) return false;
+
+    size_t tfWidth = tf->width();
+    if (tfWidth <= 0) return false;
+
+    std::vector<float> gradients(tfWidth - 1);
+    for (size_t offset = 0; offset < tfWidth - 1; offset++) {
+        glm::vec4 prevRgba = tf->sample(offset);
+        glm::vec4 nextRgba = tf->sample(offset + 1);
+
+        float colorDifference = glm::distance(prevRgba, nextRgba);
+        float alpha = (prevRgba.w + nextRgba.w) * 0.5;
+
+        gradients[offset] = colorDifference*alpha;
+    }
+
+
+
+    unsigned int nHistograms = _tsp->numTotalNodes();
+    _brickErrors = std::vector<float>(nHistograms);
+
+    for (unsigned int brickIndex = 0; brickIndex < nHistograms; brickIndex++) {
+        if (_tsp->isBstLeaf(brickIndex) && _tsp->isOctreeLeaf(brickIndex)) {
+            _brickErrors[brickIndex] = 0;
+        } else {
+            const Histogram* histogram = _histogramManager->getHistogram(brickIndex);
+            float error = 0;
+            for (int i = 0; i < gradients.size(); i++) {
+                float x = (i + 0.5) / tfWidth;
+                float sample = histogram->interpolate(x);
+                assert(sample >= 0);
+                assert(gradients[i] >= 0);
+                error += sample * gradients[i];
+            }
+            _brickErrors[brickIndex] = error;
+        }
+    }
+
+    return true;
+}
+
+int TfBrickSelector::linearCoords(int x, int y, int z) {
+    const TSP::Header &header = _tsp->header();
+    return x + (header.xNumBricks_ * y) + (header.xNumBricks_ * header.yNumBricks_ * z);
+}
+
+void TfBrickSelector::writeSelection(BrickSelection brickSelection, std::vector<int>& bricks) {
+    BrickCover coveredBricks = brickSelection.cover;
+    for (int z = coveredBricks.lowZ; z < coveredBricks.highZ; z++) {
+        for (int y = coveredBricks.lowY; y < coveredBricks.highY; y++) {
+            for (int x = coveredBricks.lowX; x < coveredBricks.highX; x++) {
+                bricks[linearCoords(x, y, z)] = brickSelection.brickIndex;
+            }
+        }
+    }
+}
+
+
+} // namespace openspace

modules/multiresvolume/rendering/tfbrickselector.h

@@ -0,0 +1,73 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2015                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __TFBRICKSELECTOR_H__
+#define __TFBRICKSELECTOR_H__
+
+#include <vector>
+#include <modules/multiresvolume/rendering/brickselection.h>
+#include <modules/multiresvolume/rendering/brickselector.h>
+#include <modules/multiresvolume/rendering/brickcover.h>
+
+
+namespace openspace {
+
+class TSP;
+class ErrorHistogramManager;
+class TransferFunction;
+
+class TfBrickSelector : public BrickSelector {
+public:
+    TfBrickSelector(TSP* tsp, ErrorHistogramManager* hm, TransferFunction* tf, int memoryBudget, int streamingBudget);
+
+    ~TfBrickSelector();
+
+    virtual bool initialize();
+
+    void selectBricks(int timestep, std::vector<int>& bricks);
+    void setMemoryBudget(int memoryBudget);
+    void setStreamingBudget(int streamingBudget);
+    bool calculateBrickErrors();
+ private:
+
+    TSP* _tsp;
+    ErrorHistogramManager* _histogramManager;
+    TransferFunction* _transferFunction;
+    std::vector<float> _brickErrors;
+    float spatialSplitPoints(unsigned int brickIndex);
+    float temporalSplitPoints(unsigned int brickIndex);
+    float splitPoints(unsigned int brickIndex, BrickSelection::SplitType& splitType);
+
+    int linearCoords(int x, int y, int z);
+    void writeSelection(BrickSelection coveredBricks, std::vector<int>& bricks);
+
+    int _memoryBudget;
+    int _streamingBudget;
+    
+};
+
+} // namespace openspace
+
+#endif // __TFBRICKSELECTOR_H__
+

modules/multiresvolume/rendering/tsp.cpp

@@ -0,0 +1,753 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#include <sgct.h>
+#include <modules/multiresvolume/rendering/tsp.h>
+
+// ghoul
+#include <ghoul/filesystem/filesystem.h>
+#include <ghoul/filesystem/cachemanager.h>
+#include <ghoul/logging/logmanager.h>
+
+// std
+#include <algorithm>
+#include <queue>
+
+namespace {
+	const std::string _loggerCat = "TSP";
+}
+
+namespace openspace {
+
+TSP::TSP(const std::string& filename)
+	: _filename(filename)
+	, _dataSSBO(0)
+	, paddedBrickDim_(0)
+	, numTotalNodes_(0)
+	, numBSTLevels_(0)
+	, numBSTNodes_(0)
+	, numOTLevels_(0)
+	, numOTNodes_(0)
+	, minSpatialError_(0.0f)
+	, maxSpatialError_(0.0f)
+	, medianSpatialError_(0.0f)
+	, minTemporalError_(0.0f)
+	, maxTemporalError_(0.0f)
+	, medianTemporalError_(0.0f)
+{
+	_file.open(_filename, std::ios::in | std::ios::binary);
+}
+
+TSP::~TSP() {
+	if (_file.is_open())
+		_file.close();
+}
+
+bool TSP::load() {
+	if (!readHeader()) {
+		LERROR("Could not read header");
+		return false;
+	}
+
+	if (readCache()) {
+	//if (false) {
+		LINFO("Using cache");
+	}
+	else {
+		if (!construct()) {
+			LERROR("Could not construct");
+			return false;
+		}
+
+		if (false) {
+			if (!calculateSpatialError()) {
+				LERROR("Could not calculate spatial error");
+				return false;
+			}
+			if (!calculateTemporalError()) {
+				LERROR("Could not calculate temporal error");
+				return false;
+			}
+			if (!writeCache()) {
+				LERROR("Could not write cache");
+				return false;
+			}
+		}
+	}
+	initalizeSSO();
+
+	return true;
+}
+
+bool TSP::readHeader() {
+
+	if (!_file.good())
+		return false;
+
+	_file.seekg(_file.beg);
+
+	_file.read(reinterpret_cast<char*>(&_header), sizeof(Header));
+	/*
+	file.read(reinterpret_cast<char*>(&gridType_),			sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&numOrigTimesteps_),	sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&numTimesteps_),		sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&xBrickDim_),			sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&yBrickDim_),			sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&zBrickDim_),			sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&xNumBricks_),		sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&yNumBricks_),		sizeof(unsigned int));
+	file.read(reinterpret_cast<char*>(&zNumBricks_),		sizeof(unsigned int));
+	*/
+
+	LDEBUG("Grid type: " << _header.gridType_);
+	LDEBUG("Brick dimensions: " << _header.xBrickDim_ << " " << _header.yBrickDim_ << " " << _header.zBrickDim_);
+	LDEBUG("Num bricks: " << _header.xNumBricks_ << " " << _header.yNumBricks_ << " " << _header.zNumBricks_);
+
+	paddedBrickDim_ = _header.xBrickDim_ + 2 * paddingWidth_;
+	// TODO support dimensions of different size
+	numOTLevels_ = static_cast<unsigned int>(log((int)_header.xNumBricks_) / log(2) + 1);
+	numOTNodes_ = static_cast<unsigned int>((pow(8, numOTLevels_) - 1) / 7);
+	numBSTLevels_ = static_cast<unsigned int>(log((int)_header.numTimesteps_) / log(2) + 1);
+	numBSTNodes_ = static_cast<unsigned int>(_header.numTimesteps_ * 2 - 1);
+	numTotalNodes_ = numOTNodes_ * numBSTNodes_;
+
+	LDEBUG("Num OT levels: " << numOTLevels_);
+	LDEBUG("Num OT nodes: " << numOTNodes_);
+	LDEBUG("Num BST levels: " << numBSTLevels_);
+	LDEBUG("Num BST nodes: " << numBSTNodes_);
+	LDEBUG("NUm total nodes: " << numTotalNodes_);
+
+	// Allocate space for TSP structure
+	data_.resize(numTotalNodes_*NUM_DATA);
+	LDEBUG("data size: " << data_.size());
+
+	return true;
+}
+
+bool TSP::construct() {
+
+	LDEBUG("Constructing TSP tree");
+
+	// Loop over the OTs (one per BST node)
+	for (unsigned int OT = 0; OT<numBSTNodes_; ++OT) {
+
+		// Start at the root of each OT  
+		unsigned int OTNode = OT*numOTNodes_;
+
+		// Calculate BST level (first level is level 0)
+		unsigned int BSTLevel = static_cast<unsigned int>(log(OT + 1) / log(2));
+
+		// Traverse OT
+		unsigned int OTChild = 1;
+		unsigned int OTLevel = 0;
+		while (OTLevel < numOTLevels_) {
+
+			unsigned int OTNodesInLevel = static_cast<unsigned int>(pow(8, OTLevel));
+			for (unsigned int i = 0; i<OTNodesInLevel; ++i) {
+
+				// Brick index
+				data_[OTNode*NUM_DATA + BRICK_INDEX] = (int)OTNode;
+
+				// Error metrics
+				//int localOTNode = (OTNode - OT*numOTNodes_);
+				data_[OTNode*NUM_DATA + TEMPORAL_ERR] = static_cast<int>(numBSTLevels_ - 1 - BSTLevel);
+				data_[OTNode*NUM_DATA + SPATIAL_ERR] = static_cast<int>(numOTLevels_ - 1 - OTLevel);
+
+				if (BSTLevel == 0) {
+					// Calculate OT child index (-1 if node is leaf)
+					int OTChildIndex =
+						(OTChild < numOTNodes_) ? static_cast<int>(OT*numOTNodes_ + OTChild) : -1;
+					data_[OTNode*NUM_DATA + CHILD_INDEX] = OTChildIndex;
+				}
+				else {
+					// Calculate BST child index (-1 if node is BST leaf)
+
+					// First BST node of current level
+					int firstNode =
+						static_cast<unsigned int>((2 * pow(2, BSTLevel - 1) - 1)*numOTNodes_);
+					// First BST node of next level
+					int firstChild =
+						static_cast<unsigned int>((2 * pow(2, BSTLevel) - 1)*numOTNodes_);
+					// Difference between first nodes between levels
+					int levelGap = firstChild - firstNode;
+					// How many nodes away from the first node are we?
+					int offset = (OTNode - firstNode) / numOTNodes_;
+
+					// Use level gap and offset to calculate child index
+					int BSTChildIndex =
+						(BSTLevel < numBSTLevels_ - 1) ?
+						static_cast<int>(OTNode + levelGap + (offset*numOTNodes_)) : -1;
+
+					data_[OTNode*NUM_DATA + CHILD_INDEX] = BSTChildIndex;
+				}
+
+				OTNode++;
+				OTChild += 8;
+
+			}
+
+			OTLevel++;
+		}
+	}
+	return true;
+}
+
+bool TSP::initalizeSSO() {
+
+	if (!_dataSSBO)
+		glGenBuffers(1, &_dataSSBO);
+
+	const size_t size = sizeof(GLint)*data_.size();
+	glBindBuffer(GL_SHADER_STORAGE_BUFFER, _dataSSBO);
+	//glBufferData(GL_SHADER_STORAGE_BUFFER, size, data_.data(), GL_DYNAMIC_READ);
+	glBufferData(GL_SHADER_STORAGE_BUFFER, size, data_.data(), GL_STATIC_DRAW);
+	glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
+	glFinish();
+	return true;
+}
+
+const TSP::Header& TSP::header() const {
+	return _header;
+}
+
+long long TSP::dataPosition() {
+	return sizeof(Header);
+}
+
+std::ifstream& TSP::file() {
+	return _file;
+}
+
+unsigned int TSP::numTotalNodes() const { 
+	return numTotalNodes_; 
+}
+
+unsigned int TSP::numValuesPerNode() const { 
+	return NUM_DATA; 
+}
+
+unsigned int TSP::numBSTNodes() const { 
+	return numBSTNodes_; 
+}
+
+unsigned int TSP::numBSTLevels() const {
+	return numBSTLevels_;
+}
+
+unsigned int TSP::numOTNodes() const { 
+	return numOTNodes_; 
+}
+
+unsigned int TSP::numOTLevels() const { 
+	return numOTLevels_; 
+}
+
+unsigned int TSP::brickDim() const {
+	return _header.xBrickDim_;
+}
+
+unsigned int TSP::paddedBrickDim() const {
+	return paddedBrickDim_;
+}
+
+unsigned int TSP::numBricksPerAxis() const {
+	return _header.xNumBricks_;
+}
+
+GLuint TSP::ssbo() const {
+	return _dataSSBO;
+}
+
+bool TSP::calculateSpatialError() {
+	unsigned int numBrickVals = paddedBrickDim_*paddedBrickDim_*paddedBrickDim_;
+
+	if (!_file.is_open())
+		return false;
+
+	std::vector<float> buffer(numBrickVals);
+	std::vector<float> averages(numTotalNodes_);
+	std::vector<float> stdDevs(numTotalNodes_);
+
+	// First pass: Calculate average color for each brick
+	LDEBUG("Calculating spatial error, first pass");
+	for (unsigned int brick = 0; brick<numTotalNodes_; ++brick) {
+
+		// Offset in file
+		std::streampos offset = dataPosition() + static_cast<long long>(brick*numBrickVals*sizeof(float));
+		_file.seekg(offset);
+
+		_file.read(reinterpret_cast<char*>(&buffer[0]),
+			static_cast<size_t>(numBrickVals)*sizeof(float));
+
+		double average = 0.0;
+		for (auto it = buffer.begin(); it != buffer.end(); ++it) {
+			average += *it;
+		}
+
+		averages[brick] = average / static_cast<double>(numBrickVals);
+	}
+
+	// Spatial SNR stats
+	float minError = 1e20f;
+	float maxError = 0.f;
+	std::vector<float> medianArray(numTotalNodes_);
+
+	// Second pass: For each brick, compare the covered leaf voxels with
+	// the brick average
+	LDEBUG("Calculating spatial error, second pass");
+	for (unsigned int brick = 0; brick<numTotalNodes_; ++brick) {
+
+		// Fetch mean intensity 
+		float brickAvg = averages[brick];
+
+		// Sum  for std dev computation
+		float stdDev = 0.f;
+
+		// Get a list of leaf bricks that the current brick covers
+		std::list<unsigned int> coveredLeafBricks =
+			CoveredLeafBricks(brick);
+
+		// If the brick is already a leaf, assign a negative error.
+		// Ad hoc "hack" to distinguish leafs from other nodes that happens
+		// to get a zero error due to rounding errors or other reasons.
+		if (coveredLeafBricks.size() == 1) {
+			stdDev = -0.1f;
+		}
+		else {
+
+			// Calculate "standard deviation" corresponding to leaves
+			for (auto lb = coveredLeafBricks.begin();
+				lb != coveredLeafBricks.end(); ++lb) {
+
+				// Read brick
+				std::streampos offset = dataPosition() + static_cast<long long>((*lb)*numBrickVals*sizeof(float));
+				_file.seekg(offset);
+
+				_file.read(reinterpret_cast<char*>(&buffer[0]),
+					static_cast<size_t>(numBrickVals)*sizeof(float));
+
+				// Add to sum
+				for (auto v = buffer.begin(); v != buffer.end(); ++v) {
+					stdDev += pow(*v - brickAvg, 2.f);
+				}
+
+
+			}
+
+			// Finish calculation
+			if (sizeof(float) != sizeof(int)) {
+				LERROR("Float and int sizes don't match, can't reintepret");
+				return false;
+			}
+
+			stdDev /= static_cast<float>(coveredLeafBricks.size()*numBrickVals);
+			stdDev = sqrt(stdDev);
+
+		} // if not leaf
+
+		if (stdDev < minError) {
+			minError = stdDev;
+		}
+		else if (stdDev > maxError) {
+			maxError = stdDev;
+		}
+
+		stdDevs[brick] = stdDev;
+		medianArray[brick] = stdDev;
+
+	}
+
+	std::sort(medianArray.begin(), medianArray.end());
+	//float medError = medianArray[medianArray.size()/2];
+
+	/*
+	LDEBUG("\nMin spatial std dev: " << minError);
+	LDEBUG("Max spatial std dev: " << maxError);
+	LDEBUG("Median spatial std dev: " << medError);
+	LDEBUG("");
+	*/
+
+	// "Normalize" errors
+	float minNorm = 1e20f;
+	float maxNorm = 0.f;
+	for (unsigned int i = 0; i<numTotalNodes_; ++i) {
+		//float normalized = (stdDevs[i]-minError)/(maxError-minError);
+		if (stdDevs[i] > 0.f) {
+			stdDevs[i] = pow(stdDevs[i], 0.5f);
+		}
+		//data_[i*NUM_DATA + SPATIAL_ERR] = *reinterpret_cast<int*>(&stdDevs[i]);
+		data_[i*NUM_DATA + SPATIAL_ERR] = glm::floatBitsToInt(stdDevs[i]);
+		if (stdDevs[i] < minNorm) {
+			minNorm = stdDevs[i];
+		}
+		else if (stdDevs[i] > maxNorm) {
+			maxNorm = stdDevs[i];
+		}
+	}
+
+	std::sort(stdDevs.begin(), stdDevs.end());
+	float medNorm = stdDevs[stdDevs.size() / 2];
+
+	minSpatialError_ = minNorm;
+	maxSpatialError_ = maxNorm;
+	medianSpatialError_ = medNorm;
+
+	LDEBUG("Min normalized spatial std dev: " << minNorm);
+	LDEBUG("Max normalized spatial std dev: " << maxNorm);
+	LDEBUG("Median normalized spatial std dev: " << medNorm);
+
+	return true;
+}
+
+bool TSP::calculateTemporalError() {
+
+	if (!_file.is_open())
+		return false;
+
+	LDEBUG("Calculating temporal error");
+
+	// Statistics
+	//float minErr = 1e20f;
+	//float maxErr = 0.f;
+	std::vector<float> meanArray(numTotalNodes_);
+
+	// Save errors
+	std::vector<float> errors(numTotalNodes_);
+
+	// Calculate temporal error for one brick at a time
+	for (unsigned int brick = 0; brick<numTotalNodes_; ++brick) {
+
+		unsigned int numBrickVals =
+			paddedBrickDim_*paddedBrickDim_*paddedBrickDim_;
+
+		// Save the individual voxel's average over timesteps. Because the
+		// BSTs are built by averaging leaf nodes, we only need to sample
+		// the brick at the correct coordinate.
+		std::vector<float> voxelAverages(numBrickVals);
+		std::vector<float> voxelStdDevs(numBrickVals);
+
+		// Read the whole brick to fill the averages
+		std::streampos offset = dataPosition() + static_cast<long long>(brick*numBrickVals*sizeof(float));
+		_file.seekg(offset);
+
+		_file.read(reinterpret_cast<char*>(&voxelAverages[0]),
+			static_cast<size_t>(numBrickVals)*sizeof(float));
+
+		// Build a list of the BST leaf bricks (within the same octree level) that
+		// this brick covers
+		std::list<unsigned int> coveredBricks = CoveredBSTLeafBricks(brick);
+
+		// If the brick is at the lowest BST level, automatically set the error 
+		// to -0.1 (enables using -1 as a marker for "no error accepted");
+		// Somewhat ad hoc to get around the fact that the error could be
+		// 0.0 higher up in the tree
+		if (coveredBricks.size() == 1) {
+			errors[brick] = -0.1f;
+		} else {
+			// Calculate standard deviation per voxel, average over brick
+			float avgStdDev = 0.f;
+			for (unsigned int voxel = 0; voxel<numBrickVals; ++voxel) {
+
+				float stdDev = 0.f;
+				for (auto leaf = coveredBricks.begin();
+					leaf != coveredBricks.end(); ++leaf) {
+
+					// Sample the leaves at the corresponding voxel position
+
+					std::streampos offset = dataPosition() + static_cast<long long>((*leaf*numBrickVals + voxel)*sizeof(float));
+					_file.seekg(offset);
+
+					float sample;
+					_file.read(reinterpret_cast<char*>(&sample), sizeof(float));
+
+					stdDev += pow(sample - voxelAverages[voxel], 2.f);
+				}
+				stdDev /= static_cast<float>(coveredBricks.size());
+				stdDev = sqrt(stdDev);
+
+				avgStdDev += stdDev;
+			} // for voxel
+
+			avgStdDev /= static_cast<float>(numBrickVals);
+			meanArray[brick] = avgStdDev;
+			errors[brick] = avgStdDev;
+
+		}
+
+		/*
+		if (avgStdDev < minErr) {
+		minErr = avgStdDev;
+		} else if (avgStdDev > maxErr) {
+		maxErr = avgStdDev;
+		}
+		*/
+
+	} // for all bricks
+
+	std::sort(meanArray.begin(), meanArray.end());
+	//float medErr = meanArray[meanArray.size()/2];
+
+	/*
+	LDEBUG("\nMin temporal error: " << minErr);
+	LDEBUG("Max temporal error: " << maxErr);
+	LDEBUG("Median temporal error: " << medErr);
+	*/
+
+	// Adjust errors using user-provided exponents
+	float minNorm = 1e20f;
+	float maxNorm = 0.f;
+	for (unsigned int i = 0; i<numTotalNodes_; ++i) {
+		if (errors[i] > 0.f) {
+			errors[i] = pow(errors[i], 0.25f);
+		}
+		//data_[i*NUM_DATA + TEMPORAL_ERR] = *reinterpret_cast<int*>(&errors[i]);
+		data_[i*NUM_DATA + TEMPORAL_ERR] = glm::floatBitsToInt(errors[i]);
+		if (errors[i] < minNorm) {
+			minNorm = errors[i];
+		}
+		else if (errors[i] > maxNorm) {
+			maxNorm = errors[i];
+		}
+	}
+
+	std::sort(errors.begin(), errors.end());
+	float medNorm = errors[errors.size() / 2];
+
+	minTemporalError_ = minNorm;
+	maxTemporalError_ = maxNorm;
+	medianTemporalError_ = medNorm;
+
+	LDEBUG("Min normalized temporal std dev: " << minNorm);
+	LDEBUG("Max normalized temporal std dev: " << maxNorm);
+	LDEBUG("Median normalized temporal std dev: " << medNorm);
+
+	return true;
+}
+
+
+bool TSP::readCache() {
+
+	if (!FileSys.cacheManager())
+		return false;
+
+	ghoul::filesystem::File f = _filename;
+	std::string cacheFilename = FileSys.cacheManager()->cachedFilename(
+        f.baseName(), "", ghoul::filesystem::CacheManager::Persistent::Yes);
+
+	std::ifstream file(cacheFilename, std::ios::in | std::ios::binary);
+	if (!file.is_open()) {
+		LWARNING("Failed to open " << cacheFilename);
+		return false;
+	}
+
+
+	file.read(reinterpret_cast<char*>(&minSpatialError_), sizeof(float));
+	file.read(reinterpret_cast<char*>(&maxSpatialError_), sizeof(float));
+	file.read(reinterpret_cast<char*>(&medianSpatialError_), sizeof(float));
+	file.read(reinterpret_cast<char*>(&minTemporalError_), sizeof(float));
+	file.read(reinterpret_cast<char*>(&maxTemporalError_), sizeof(float));
+	file.read(reinterpret_cast<char*>(&medianTemporalError_), sizeof(float));
+	size_t dataSize = static_cast<size_t>(numTotalNodes_*NUM_DATA)*sizeof(int);
+	file.read(reinterpret_cast<char*>(&data_[0]), dataSize);
+	file.close();
+
+	LDEBUG("Cached errors:");
+	LDEBUG("Min spatial error: " << minSpatialError_);
+	LDEBUG("Max spatial error: " << maxSpatialError_);
+	LDEBUG("Median spatial error: " << medianSpatialError_);
+	LDEBUG("Min temporal error: " << minTemporalError_);
+	LDEBUG("Max temporal error: " << maxTemporalError_);
+	LDEBUG("Median temporal error: " << medianTemporalError_);
+
+	return true;
+}
+
+bool TSP::writeCache() {
+
+	if (!FileSys.cacheManager())
+		return false;
+
+	ghoul::filesystem::File f = _filename;
+	std::string cacheFilename = FileSys.cacheManager()->cachedFilename(
+        f.baseName(), "", ghoul::filesystem::CacheManager::Persistent::Yes);
+
+	std::ofstream file(cacheFilename, std::ios::out | std::ios::binary);
+	if (!file.is_open()) {
+		LWARNING("Failed to open " << cacheFilename);
+		return false;
+	}
+	LINFO("Writing cache to " << cacheFilename);
+
+
+	file.write(reinterpret_cast<char*>(&minSpatialError_), sizeof(float));
+	file.write(reinterpret_cast<char*>(&maxSpatialError_), sizeof(float));
+	file.write(reinterpret_cast<char*>(&medianSpatialError_), sizeof(float));
+	file.write(reinterpret_cast<char*>(&minTemporalError_), sizeof(float));
+	file.write(reinterpret_cast<char*>(&maxTemporalError_), sizeof(float));
+	file.write(reinterpret_cast<char*>(&medianTemporalError_), sizeof(float));
+	file.write(reinterpret_cast<char*>(&data_[0]), data_.size()*sizeof(float));
+
+	file.close();
+
+	/*
+	LDEBUG("\nData:");
+	for (unsigned i=0; i<data_.size()/NUM_DATA; ++i) {
+	LDEBUG("Brick nr " << i);
+	LDEBUG("Brick index " << data_[i*NUM_DATA + BRICK_INDEX]);
+	LDEBUG("Child index " << data_[i*NUM_DATA + CHILD_INDEX]);
+	LDEBUG("Spatial err " << *reinterpret_cast<float*>(&data_[i*NUM_DATA + SPATIAL_ERR]));
+	LDEBUG("Temporal err " << *reinterpret_cast<float*>(&data_[i*NUM_DATA + TEMPORAL_ERR]));
+	}
+	*/
+
+	return true;
+}
+
+float TSP::getSpatialError(unsigned int _brickIndex) {
+	return reinterpret_cast<float &>(data_[_brickIndex*NUM_DATA + SPATIAL_ERR]);
+}
+
+float TSP::getTemporalError(unsigned int _brickIndex) {
+	return reinterpret_cast<float &>(data_[_brickIndex*NUM_DATA + TEMPORAL_ERR]);
+}
+
+unsigned int TSP::getFirstOctreeChild(unsigned int _brickIndex) {
+    unsigned int otNode = _brickIndex % numOTNodes_;
+    unsigned int bstOffset = _brickIndex - otNode;
+
+    unsigned int depth = log(7 * otNode + 1) / log(8);
+    unsigned int firstInLevel = (pow(8, depth) - 1) / 7;
+    unsigned int levelOffset = otNode - firstInLevel;
+    unsigned int firstInChildLevel = (pow(8, depth + 1) - 1) / 7;
+    unsigned int childIndex = firstInChildLevel + 8*levelOffset;
+
+    return bstOffset + childIndex;
+}
+
+unsigned int TSP::getBstLeft(unsigned int _brickIndex) {
+    unsigned int bstNode = _brickIndex / numOTNodes_;
+    unsigned int otOffset = _brickIndex % numOTNodes_;
+    unsigned int depth = log(bstNode + 1) / log(2);
+    unsigned int firstInLevel = pow(2, depth) - 1;
+    unsigned int levelOffset = bstNode - firstInLevel;
+    unsigned int firstInChildLevel = pow(2, depth + 1) - 1;
+    unsigned int childIndex = firstInChildLevel + 2*levelOffset;
+    return otOffset + childIndex * numOTNodes_;
+}
+
+unsigned int TSP::getBstRight(unsigned int _brickIndex) {
+    return getBstLeft(_brickIndex) + numOTNodes_;
+}
+
+bool TSP::isBstLeaf(unsigned int _brickIndex) {
+    unsigned int bstNode = _brickIndex / numOTNodes_;
+    return bstNode >= numBSTNodes_ / 2;
+}
+
+bool TSP::isOctreeLeaf(unsigned int _brickIndex) {
+    unsigned int otNode = _brickIndex % numOTNodes_;
+    unsigned int depth = log(7 * otNode + 1) / log(8);
+    return depth == numOTLevels_ - 1;
+}
+
+
+std::list<unsigned int> TSP::CoveredLeafBricks(unsigned int _brickIndex) {
+	std::list<unsigned int> out;
+
+	// Find what octree skeleton node the index belongs to
+	unsigned int OTNode = _brickIndex % numOTNodes_;
+
+	// Find what BST node the index corresponds to using int division
+	unsigned int BSTNode = _brickIndex / numOTNodes_;
+	// Calculate BST offset (to translate to root octree)
+	unsigned int BSTOffset = BSTNode * numOTNodes_;
+
+	// Traverse root octree structure to leaves
+	// When visiting the leaves, translate back to correct BST level and save
+	std::queue<unsigned int> queue;
+	queue.push(OTNode);
+	do {
+
+		// Get front of queue and pop it
+		unsigned int toVisit = queue.front();
+		queue.pop();
+
+		// See if the node has children
+		int child = data_[toVisit*NUM_DATA + CHILD_INDEX];
+		if (child == -1) {
+			// Translate back and save
+			out.push_back(toVisit + BSTOffset);
+		}
+		else {
+			// Queue the eight children
+			for (int i = 0; i<8; ++i) {
+				queue.push(child + i);
+			}
+		}
+
+	} while (!queue.empty());
+
+	return out;
+}
+
+std::list<unsigned int> TSP::CoveredBSTLeafBricks(unsigned int _brickIndex) {
+	std::list<unsigned int> out;
+
+	// Traverse the BST children until we are at root
+	std::queue<unsigned int> queue;
+	queue.push(_brickIndex);
+	do {
+
+		unsigned int toVisit = queue.front();
+		queue.pop();
+
+		bool BSTRoot = toVisit < numOTNodes_;
+		if (BSTRoot) {
+			if (numBSTLevels_ == 1) {
+				out.push_back(toVisit);
+			}
+			else {
+				queue.push(toVisit + numOTNodes_);
+				queue.push(toVisit + numOTNodes_ * 2);
+			}
+		}
+		else {
+			int child = data_[toVisit*NUM_DATA + CHILD_INDEX];
+			if (child == -1) {
+				// Save leaf brick to list
+				out.push_back(toVisit);
+			}
+			else {
+				// Queue children
+				queue.push(child);
+				queue.push(child + numOTNodes_);
+			}
+		}
+
+	} while (!queue.empty());
+
+	return out;
+}
+
+}

modules/multiresvolume/rendering/tsp.h

@@ -0,0 +1,148 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+#ifndef __TSP_H__
+#define __TSP_H__
+
+// std includes
+#include <string>
+#include <vector>
+#include <list>
+#include <iostream>
+#include <fstream>
+
+// ghoul includes
+#include <ghoul/opengl/ghoul_gl.h>
+
+namespace openspace {
+class TSP {
+public:
+
+	struct Header {
+		unsigned int gridType_;
+		unsigned int numOrigTimesteps_;
+		unsigned int numTimesteps_;
+		unsigned int xBrickDim_;
+		unsigned int yBrickDim_;
+		unsigned int zBrickDim_;
+		unsigned int xNumBricks_;
+		unsigned int yNumBricks_;
+		unsigned int zNumBricks_;
+
+	};
+
+	enum NodeData {
+		BRICK_INDEX = 0,
+		CHILD_INDEX,
+		SPATIAL_ERR,
+		TEMPORAL_ERR,
+		NUM_DATA
+	};
+
+	TSP(const std::string& filename);
+	~TSP();
+
+	// load performs readHeader, readCache, writeCache and construct 
+	// in the correct sequence
+	bool load();
+
+	bool readHeader();
+	bool readCache();
+	bool writeCache();
+	bool construct();
+	bool initalizeSSO();
+
+	const Header& header() const;
+	static long long dataPosition();
+	std::ifstream& file();
+	unsigned int numTotalNodes() const;
+	unsigned int numValuesPerNode() const;
+	unsigned int numBSTNodes() const;
+	unsigned int numBSTLevels() const;
+	unsigned int numOTNodes() const;
+	unsigned int numOTLevels() const;
+	unsigned int brickDim() const;
+	unsigned int paddedBrickDim() const;
+	unsigned int numBricksPerAxis() const;
+	GLuint ssbo() const;
+
+	bool calculateSpatialError();
+	bool calculateTemporalError();
+
+	float getSpatialError(unsigned int _brickIndex);
+	float getTemporalError(unsigned int _brickIndex);
+	unsigned int getFirstOctreeChild(unsigned int _brickIndex);
+
+	unsigned int getBstLeft(unsigned int _brickIndex);
+	unsigned int getBstRight(unsigned int _brickIndex);
+
+	bool isBstLeaf(unsigned int _brickIndex);
+	bool isOctreeLeaf(unsigned int _brickIndex);
+
+private:
+	// Returns a list of the octree leaf nodes that a given input 
+	// brick covers. If the input is already a leaf, the list will
+	// only contain that one index.
+	std::list<unsigned int> CoveredLeafBricks(unsigned int _brickIndex);
+
+	// Returns a list of the BST leaf nodes that a given input brick
+	// covers (at the same spatial subdivision level).
+	std::list<unsigned int> CoveredBSTLeafBricks(unsigned int _brickIndex);
+
+	// Return a list of eight children brick incices given a brick index
+	std::list<unsigned int> ChildBricks(unsigned int _brickIndex);
+
+	std::string _filename;
+	std::ifstream _file;
+	std::streampos _dataOffset;
+
+	// Holds the actual structure
+	std::vector<int> data_;
+	GLuint _dataSSBO;
+
+	// Data from file
+	Header _header;
+
+	// Additional metadata
+	unsigned int paddedBrickDim_;
+	unsigned int numTotalNodes_;
+	unsigned int numBSTLevels_;
+	unsigned int numBSTNodes_;
+	unsigned int numOTLevels_;
+	unsigned int numOTNodes_;
+
+	const unsigned int paddingWidth_ = 1;
+
+	// Error stats
+	float minSpatialError_;
+	float maxSpatialError_;
+	float medianSpatialError_;
+	float minTemporalError_;
+	float maxTemporalError_;
+	float medianTemporalError_;
+
+}; // class TSP
+}  // namespace openspace
+
+#endif

modules/multiresvolume/shaders/boundsFs.glsl

@@ -0,0 +1,40 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2016                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+in vec3 vPosition;
+in vec4 worldPosition;
+
+#include "PowerScaling/powerScaling_fs.hglsl"
+#include "fragment.glsl"
+
+Fragment getFragment() {
+	vec4 fragColor = vec4(vPosition+0.5, 1.0);
+	vec4 position = worldPosition;
+	float depth = pscDepth(position);
+
+    Fragment frag;
+    frag.color = fragColor;
+    frag.depth = depth;
+    return frag;
+}

modules/multiresvolume/shaders/boundsVs.glsl

@@ -0,0 +1,47 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2016                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+ 
+#version __CONTEXT__
+
+layout(location = 0) in vec4 vertPosition;
+
+uniform mat4 viewProjection;
+uniform mat4 modelTransform;
+
+out vec3 vPosition;
+out vec4 worldPosition;
+
+#include "PowerScaling/powerScaling_vs.hglsl"
+
+void main() {
+	vPosition = vertPosition.xyz;
+	worldPosition = modelTransform*vertPosition;
+	
+	vec4 position = pscTransform(worldPosition, mat4(1.0));
+    
+	// project the position to view space
+	gl_Position = viewProjection * position;
+
+    gl_Position.z = 1.0;
+}

modules/multiresvolume/shaders/helper.glsl

@@ -0,0 +1,31 @@
+#define MULTIRES_PI      3.14159265358979323846  /* pi */
+#define MULTIRES_SQRT1_3 0.57735026919           /* 1/sqrt(3) */
+#define MULTIRES_OPACITY_THRESHOLD 0.01
+vec3 multires_cartesianToSpherical(vec3 _cartesian) {
+    // Put cartesian in [-1..1] range first
+    vec3 cartesian = vec3(-1.0,-1.0,-1.0) + _cartesian * 2.0f;
+
+    float r = length(cartesian);
+    float theta, phi;
+
+    if (r == 0.0) {
+        theta = phi = 0.0;
+    } else {
+        theta = acos(cartesian.z/r) / MULTIRES_PI;
+        phi = (MULTIRES_PI + atan(cartesian.y, cartesian.x)) / (2.0*MULTIRES_PI );
+    }
+    r *= MULTIRES_SQRT1_3;
+    return vec3(r, theta, phi);
+}
+
+int multires_intCoord(ivec3 vec3Coords, ivec3 spaceDim) {
+    return vec3Coords.x + spaceDim.x*vec3Coords.y + spaceDim.x*spaceDim.y*vec3Coords.z;
+}
+
+vec3 multires_vec3Coords(uint intCoord, ivec3 spaceDim) {
+    vec3 coords = vec3(0.0);
+    coords.x = mod(intCoord, spaceDim.x);
+    coords.y = mod(intCoord / spaceDim.x, spaceDim.y);
+    coords.z = intCoord / spaceDim.x / spaceDim.y;
+    return coords;
+}

modules/multiresvolume/shaders/raycast.glsl

@@ -0,0 +1,123 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2016                                                                    *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+
+uniform float opacity_#{id};
+uniform sampler1D transferFunction_#{id};
+uniform sampler3D textureAtlas_#{id};
+uniform int gridType_#{id};
+uniform uint maxNumBricksPerAxis_#{id};
+uniform uint paddedBrickDim_#{id};
+uniform ivec3 nBricksInAtlas_#{id};
+uniform ivec3 atlasSize_#{id};
+uniform float stepSizeCoefficient_#{id} = 1.0;
+
+layout (shared) buffer atlasMapBlock_#{id} {
+    uint atlasMap_#{id}[];
+};
+
+void atlasMapDataFunction_#{id}(ivec3 brickCoords, inout uint atlasIntCoord, inout uint level) {
+    int linearBrickCoord = multires_intCoord(brickCoords, ivec3(maxNumBricksPerAxis_#{id}));
+    uint mapData = atlasMap_#{id}[linearBrickCoord];
+    level = mapData >> 28;
+    atlasIntCoord = mapData & 0x0FFFFFFF;
+}
+
+vec3 atlasCoordsFunction_#{id}(vec3 position) {
+    uint maxNumBricksPerAxis = maxNumBricksPerAxis_#{id};
+    uint paddedBrickDim = paddedBrickDim_#{id};
+
+    ivec3 brickCoords = ivec3(position * maxNumBricksPerAxis);
+    uint atlasIntCoord, level;
+    atlasMapDataFunction_#{id}(brickCoords, atlasIntCoord, level);
+
+    float levelDim = float(maxNumBricksPerAxis) / pow(2.0, level);
+    vec3 inBrickCoords = mod(position*levelDim, 1.0);
+
+    float scale = float(paddedBrickDim) - 2.0;
+    vec3 paddedInBrickCoords = (1.0 + inBrickCoords * scale) / paddedBrickDim;
+
+    ivec3 numBricksInAtlas = ivec3(vec3(atlasSize_#{id}) / paddedBrickDim);
+    vec3 atlasOffset = multires_vec3Coords(atlasIntCoord, numBricksInAtlas);
+    return (atlasOffset + paddedInBrickCoords) / vec3(numBricksInAtlas);
+}
+
+float stepSize#{id}(vec3 samplePos, vec3 dir){
+    return 0.01;
+    if (true /*opacity_#{id} >= MULTIRES_OPACITY_THRESHOLD*/) {
+        return stepSizeCoefficient_#{id}/float(maxNumBricksPerAxis_#{id})/float(paddedBrickDim_#{id});
+    } else {
+        // return a number that is garantueed to be bigger than the whole volume
+        return 2.0;
+    }
+}
+
+vec4 sample#{id}(vec3 samplePos, vec3 dir, vec4 foregroundColor, inout float maxStepSize) {
+    //return vec4(1.0, 1.0, 1.0, 1.0);
+    
+    if (true /*opacity_#{id} >= MULTIRES_OPACITY_THRESHOLD*/) {
+        if (gridType_#{id} == 1) {
+            samplePos = multires_cartesianToSpherical(samplePos);
+        }
+        vec3 sampleCoords = atlasCoordsFunction_#{id}(samplePos);
+        //return vec4(sampleCoords, 1.0);
+        //sampleCoords = vec3(1.0,0.0, 0.0);
+        float intensity = texture(textureAtlas_#{id}, sampleCoords).x;
+        //intensity = sampleCoords;
+
+        //return vec4(vec3(intensity), 1.0);
+        vec4 contribution = texture(transferFunction_#{id}, intensity);
+
+        //contribution = vec4(sampleCoords, 1.0);
+        //vec4 contribution = vec4(vec3(intensity), 1.0);
+        //contribution.a *= 0.3;
+        //contribution = vec4(1.0, 1.0, 1.0, intensity * 1000000.0);
+        //contribution = vec4(1.0, 1.0, 1.0, 1.0);
+        maxStepSize = stepSizeCoefficient_#{id}/float(maxNumBricksPerAxis_#{id})/float(paddedBrickDim_#{id});
+        //contribution.a *= opacity_#{id};
+
+        //maxStepSize = 0.01;
+        return contribution;
+    } else {
+        maxStepSize = 2.0;
+        return vec4(0.0);
+    }
+}
+
+
+
+/*uniform vec4 color#{id};
+uniform float time#{id};
+uniform float maxStepSize#{id} = 0.02;
+
+vec4 sample#{id}(vec3 samplePos, vec3 dir, vec4 foregroundColor, inout float maxStepSize) {
+    maxStepSize = 0.01;
+    return vec4(1.0, 0.0, 0.0, 0.9999999);
+}
+
+float stepSize#{id}(vec3 samplePos, vec3 dir) {
+    return 0.01;
+}
+*/
+

modules/toyvolume/rendering/renderabletoyvolume.cpp

@@ -77,7 +77,7 @@ namespace openspace {
 RenderableToyVolume::~RenderableToyVolume() {}
 
 bool RenderableToyVolume::initialize() {
-    _raycaster = std::make_unique<ToyVolumeRaycaster>(ToyVolumeRaycaster(_color));
+    _raycaster = std::make_unique<ToyVolumeRaycaster>(_color);
     _raycaster->initialize();
 
     OsEng.renderEngine().raycasterManager().attachRaycaster(*_raycaster.get());

shaders/abuffer/abufferfragment.glsl

@@ -26,29 +26,31 @@
 #define _ABUFFERFRAGMENT_GLSL_
 
 struct ABufferFragment {
-    uint rgba;
-    uint depth;
-    uint data;
-    uint composition;
+    uint channel0;
+    uint channel1;
+    uint channel2;
+    uint channel3;
 };
 
+// An abuffer frame
 // Values stored in abuffer:
-// -------RGBA--------
-// r            8 bits
-// g            8 bits
-// b            8 bits
-// a            8 bits
-// -------DEPTH-------
-// depth       32 bits 
-// -------DATA--------
+//            -----CHANNEL 0-----
 // type         8 bits (signed char) 0: geometry, >0: volume entry, <0: volume exit
-// msaa         8 bits 
-// reserved    16 bits
-// ----COMPOSITION----
-// reserved     4 bits (may be suitable for blend modes)
-// next        28 bits
-// -------------------
-// in total: 16 + 4 = 20 reserved bits for future use.
+// msaa         8 bits
+// --- if geometry --- | ---- if volume ----
+// reserved    16 bits | x           16 bits
+//            -----CHANNEL 1-----
+// r            8 bits | y           16 bits
+// g            8 bits |
+// b            8 bits | z           16 bits
+// a            8 bits |
+//            -----CHANNEL 2-----
+//            depth       32 bits
+//            -----CHANNEL 3-----
+//            blend        4 bits (0 = normal, 1 = additive, 2-7 reserved)
+//            next        28 bits
+//            -------------------
+
 
 const uint mask_1 = uint(1);
 const uint mask_8 = uint(255);
@@ -58,14 +60,17 @@ const uint mask_28 = uint(268435455);
 const uint mask_31 = uint(2147483647);
 const uint mask_32 = uint(4294967295);
 
+const uint mask_red    = mask_16;
+const uint shift_red   = 0;
+
 const uint mask_type   = mask_32 - mask_24;
 const uint shift_type  = 24;
 
 const uint mask_msaa   = mask_24 - mask_16;
 const uint shift_msaa  = 16;
 
-const uint mask_blend  = mask_32 - mask_31;
-const uint shift_blend = 31;
+const uint mask_blend  = mask_32 - mask_28;
+const uint shift_blend = 28;
 
 const uint mask_next   = mask_28;
 const uint shift_next  = 0;
@@ -82,22 +87,37 @@ uint bitextract(in uint pack, uint mask, uint shift) {
  * Color
  */
 void _color_(inout ABufferFragment frag, vec4 color) {
-    frag.rgba = packUnorm4x8(color);
+    frag.channel1 = packUnorm4x8(color);
 }
 
 vec4 _color_(ABufferFragment frag) {
-    return unpackUnorm4x8(frag.rgba);
+    return unpackUnorm4x8(frag.channel1);
+}
+
+/**
+ * Position
+ */
+void _position_(inout ABufferFragment frag, vec3 pos) {
+    uint red = uint(round(clamp(pos.x, 0.0, 1.0) * 65535.0));
+    bitinsert(frag.channel0, red, mask_red, shift_red);
+    frag.channel1 = packUnorm2x16(pos.yz);
+}
+
+vec3 _position_(ABufferFragment frag) {
+    float x = float(bitextract(frag.channel0, mask_red, shift_red)) / 65535.0;
+    vec2 yz = unpackUnorm2x16(frag.channel1);
+    return vec3(x, yz);
 }
 
 /**
  * Depth
  */
 void _depth_(inout ABufferFragment frag, float depth) {
-    frag.depth = floatBitsToUint(depth);
+    frag.channel2 = floatBitsToUint(depth);
 }
 
 float _depth_(ABufferFragment frag) {
-    return uintBitsToFloat(frag.depth);
+    return uintBitsToFloat(frag.channel2);
 }
 
 /**
@@ -110,11 +130,11 @@ void _type_(inout ABufferFragment frag, int type) {
     } else {
         val = type;
     }
-    bitinsert(frag.data, val, mask_type, shift_type);
+    bitinsert(frag.channel0, val, mask_type, shift_type);
 }
 
 int _type_(ABufferFragment frag) {
-    uint val = bitextract(frag.data, mask_type, shift_type);
+    uint val = bitextract(frag.channel0, mask_type, shift_type);
     if (val > 127) {
         return 128 - int(val);
     } else {
@@ -127,11 +147,11 @@ int _type_(ABufferFragment frag) {
  */
 void _msaa_(inout ABufferFragment frag, int type) {
     uint val = uint(type);
-    bitinsert(frag.data, val, mask_msaa, shift_msaa);
+    bitinsert(frag.channel0, val, mask_msaa, shift_msaa);
 }
 
 int _msaa_(ABufferFragment frag) {
-    uint val = bitextract(frag.data, mask_msaa, shift_msaa);
+    uint val = bitextract(frag.channel0, mask_msaa, shift_msaa);
     return int(val);
 }
 
@@ -139,11 +159,23 @@ int _msaa_(ABufferFragment frag) {
  * Next
  */
 void _next_(inout ABufferFragment frag, uint val) {
-    bitinsert(frag.composition, val, mask_next, shift_next);
+    bitinsert(frag.channel3, val, mask_next, shift_next);
 }
 
 uint _next_(ABufferFragment frag) {
-    uint val = bitextract(frag.composition, mask_next, shift_next);
+    uint val = bitextract(frag.channel3, mask_next, shift_next);
+    return val;
+}
+
+/**
+ * Blend
+ */
+void _blend_(inout ABufferFragment frag, uint val) {
+    bitinsert(frag.channel3, val, mask_blend, shift_blend);
+}
+
+uint _blend_(ABufferFragment frag) {
+    uint val = bitextract(frag.channel3, mask_blend, shift_blend);
     return val;
 }
 
@@ -151,14 +183,14 @@ uint _next_(ABufferFragment frag) {
  * Raw data
  */
 void _raw_(inout ABufferFragment frag, uvec4 raw) {
-    frag.rgba = raw.x;
-    frag.depth = raw.y;
-    frag.data = raw.z;
-    frag.composition = raw.w;
+    frag.channel0 = raw.x;
+    frag.channel1 = raw.y;
+    frag.channel2 = raw.z;
+    frag.channel3 = raw.w;
 }
 
 uvec4 _raw_(inout ABufferFragment frag) {
-    return uvec4(frag.rgba, frag.depth, frag.data, frag.composition);
+    return uvec4(frag.channel0, frag.channel1, frag.channel2, frag.channel3);
 }
 
 #endif

shaders/abuffer/boundsabuffer.frag

@@ -41,8 +41,9 @@ void main() {
     uint prevHead = imageAtomicExchange(anchorPointerTexture, ivec2(gl_FragCoord.xy), newHead);
 
     ABufferFragment aBufferFrag;
-    _color_(aBufferFrag, frag.color);
+    _position_(aBufferFrag, frag.color.rgb);
     _depth_(aBufferFrag, frag.depth);
+    _blend_(aBufferFrag, frag.blend);
 
     int fragmentType = #{fragmentType};
 

shaders/abuffer/raycasterdata.glsl

@@ -30,6 +30,7 @@ struct RaycasterData {
     vec3 direction;
     float scale;
     float previousJitterDistance;
+    uint blend;
 }
 
 #endif

shaders/abuffer/renderabuffer.frag

@@ -43,6 +43,8 @@ void main() {
     ABufferFragment aBufferFrag;
     _color_(aBufferFrag, frag.color);
     _depth_(aBufferFrag, frag.depth);
+    _blend_(aBufferFrag, frag.blend);
+
     _type_(aBufferFrag, 0); // 0 = geometry type
     _msaa_(aBufferFrag, gl_SampleMaskIn[0]);
     _next_(aBufferFrag, prevHead);

shaders/abuffer/resolveabuffer.frag

@@ -39,7 +39,7 @@ uniform int nAaSamples;
 #define RAYCASTING_ENABLED #{raycastingEnabled}
 #define N_RAYCASTERS #{nRaycasters}
 #define ALPHA_LIMIT 0.99
-#define RAYCAST_MAX_STEPS 1000
+#define RAYCAST_MAX_STEPS 10000
 #define INT_MAX 2147483647
 
 /////////////////////////
@@ -132,19 +132,21 @@ void retrieveRaycasterData(uint nFrags) {
     }
     for (int i = 0; i < nFrags; i++) {
         int type = _type_(fragments[i]); // - 1;
-        vec4 color = _color_(fragments[i]);
+        vec3 position = _position_(fragments[i]);
         float depth = _depth_(fragments[i]);
+        uint blend = _blend_(fragments[i]);
         if (type > 0) { // enter raycaster
             int raycasterId = type - 1;
             if (entryDepths[raycasterId] < 0) { // first entry
-                raycasterData[raycasterId].position = color.rgb;
+                raycasterData[raycasterId].position = position;
                 raycasterData[raycasterId].previousJitterDistance = 0;
+                raycasterData[raycasterId].blend = blend;
                 entryDepths[raycasterId] = depth;
                 raycasterData[raycasterId].scale = -1;
             }
         } else if (type < 0) { // exit raycaster
             int raycasterId = -type - 1;
-            vec3 localDirection = color.xyz - raycasterData[raycasterId].position;
+            vec3 localDirection = position - raycasterData[raycasterId].position;
             raycasterData[raycasterId].direction = safeNormalize(localDirection);
             raycasterData[raycasterId].scale = safeLength(localDirection) / (depth - entryDepths[raycasterId]);
         }
@@ -183,6 +185,7 @@ void raycast(float raycastDepth, uint raycasterMask, inout vec4 finalColor) {
         nextStepSize = raycastDepth - currentDepth;
 
 #for index, raycaster in raycasters
+
         if ((raycasterMask & #{raycaster.bitmask}) != 0) {
             RaycasterData data = raycasterData[#{raycaster.id}];
             float stepSizeLocal = currentStepSize * data.scale;
@@ -195,8 +198,13 @@ void raycast(float raycastDepth, uint raycasterMask, inout vec4 finalColor) {
 
             vec4 raycasterContribution = sample#{raycaster.id}(jitteredPosition, data.direction, finalColor, maxStepSizeLocal);
             float sampleDistance = jitteredStepSizeLocal + data.previousJitterDistance;
+            uint blend = raycasterData[#{raycaster.id}].blend;
 
-            blendStep(finalColor, raycasterContribution, sampleDistance);
+            if (blend == BLEND_MODE_NORMAL) {
+                normalBlendStep(finalColor, raycasterContribution, sampleDistance);
+            } else if (blend == BLEND_MODE_ADDITIVE) {
+                additiveBlendStep(finalColor, raycasterContribution, sampleDistance);
+            }
 
             raycasterData[#{raycaster.id}].previousJitterDistance = stepSizeLocal - jitteredStepSizeLocal;
             float maxStepSize = maxStepSizeLocal/data.scale;
@@ -224,10 +232,15 @@ void main() {
         ABufferFragment frag = fragments[i];
 
         int type = _type_(frag);
+        uint blend = _blend_(frag);
 
         if (type == 0) { // geometry fragment
             vec4 color = _color_(frag);
-            blend(finalColor, color);
+            if (blend == BLEND_MODE_NORMAL) {
+                normalBlend(finalColor, color);
+            } else if (blend == BLEND_MODE_ADDITIVE) {
+                additiveBlend(finalColor, color);
+            }
         }
 #if RAYCASTING_ENABLED
         else if (type > 0) { // enter volume

shaders/blending.glsl

@@ -25,6 +25,52 @@
 #ifndef _BLENDING_GLSL_
 #define _BLENDING_GLSL_
 
+
+/**
+ * Blend in src behind dst using normal blending
+ * dst is premultiplied
+ * src is expressed in straight RGBA
+ */
+void normalBlend(inout vec4 dst, vec4 src) {
+    dst.rgb = dst.rgb + (1.0 - dst.a) * src.a * src.rgb;
+    dst.a = dst.a + (1.0 - dst.a) * src.a;
+}
+
+/**
+ * Blend in src behind dst using additive blending
+ * dst is premultiplied
+ * src is expressed in straight RGBA
+ */
+void additiveBlend(inout vec4 dst, vec4 src) {
+    dst.rgb = dst.rgb + (1.0 - dst.a) * src.a * src.rgb;//    dst.rgb = dst.rgb + src.a * src.rgb;
+    //dst.a = dst.a + src.a;
+}
+
+
+/**
+ * Blend in src behind dst using normal blending
+ * dst is premultiplied
+ * src is expressed in straight RGBA
+ * stepSize = 0: alpha becomes 0
+ * stepSize = 1: alpha becomes src.a
+ */
+void normalBlendStep(inout vec4 dst, vec4 src, float stepSize) {
+    src.a = 1.0 - pow(1.0 - src.a, stepSize);
+    normalBlend(dst, src);
+}
+
+/**
+ * Blend in src behind dst using addivtive blending
+ * dst is premultiplied
+ * src is expressed in straight RGBA
+ * stepSize = 0: alpha becomes 0
+ * stepSize = 1: alpha becomes src.a
+ */
+void additiveBlendStep(inout vec4 dst, vec4 src, float stepSize) {
+    src.a = 1.0 - pow(1.0 - src.a, stepSize);
+    additiveBlend(dst, src);
+}
+
 /**
  * Blend in src behind dst
  * dst is premultiplied
@@ -47,4 +93,4 @@ void blendStep(inout vec4 dst, vec4 src, float stepSize) {
     blend(dst, src);
 }
 
-#endif
+#endif

shaders/fragment.glsl

@@ -25,10 +25,13 @@
 #ifndef _FRAGMENT_GLSL_
 #define _FRAGMENT_GLSL_
 
+#define BLEND_MODE_NORMAL 0
+#define BLEND_MODE_ADDITIVE 1
+
 struct Fragment {
     vec4 color;
     float depth;
-    int blend;
+    uint blend;
 };
 
 #endif    

src/CMakeLists.txt

@@ -64,6 +64,7 @@ set(OPENSPACE_SOURCE
     ${OPENSPACE_BASE_DIR}/src/rendering/renderable.cpp
     ${OPENSPACE_BASE_DIR}/src/rendering/renderengine.cpp
     ${OPENSPACE_BASE_DIR}/src/rendering/renderengine_lua.inl
+    ${OPENSPACE_BASE_DIR}/src/rendering/transferfunction.cpp
     ${OPENSPACE_BASE_DIR}/src/scene/ephemeris.cpp
     ${OPENSPACE_BASE_DIR}/src/scene/scene.cpp
     ${OPENSPACE_BASE_DIR}/src/scene/scene_lua.inl
@@ -135,6 +136,7 @@ set(OPENSPACE_HEADER
     ${OPENSPACE_BASE_DIR}/include/openspace/rendering/renderer.h
     ${OPENSPACE_BASE_DIR}/include/openspace/rendering/renderengine.h
     ${OPENSPACE_BASE_DIR}/include/openspace/rendering/volumeraycaster.h
+    ${OPENSPACE_BASE_DIR}/include/openspace/rendering/transferfunction.h
     ${OPENSPACE_BASE_DIR}/include/openspace/scene/ephemeris.h
     ${OPENSPACE_BASE_DIR}/include/openspace/scene/scene.h
     ${OPENSPACE_BASE_DIR}/include/openspace/scene/scenegraph.h

src/rendering/framebufferrenderer.cpp

@@ -220,12 +220,12 @@ void FramebufferRenderer::updateResolution() {
     glTexImage2D(
         GL_TEXTURE_2D,
         0,
-        GL_RGBA,
+        GL_RGBA16,
         GLsizei(_resolution.x),
         GLsizei(_resolution.y),
         0,
         GL_RGBA,
-        GL_BYTE,
+        GL_UNSIGNED_SHORT,
         nullptr);
 
     glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

src/rendering/transferfunction.cpp

@@ -0,0 +1,227 @@
+/*****************************************************************************************
+ *                                                                                       *
+ * OpenSpace                                                                             *
+ *                                                                                       *
+ * Copyright (c) 2014-2015                                                               *
+ *                                                                                       *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this  *
+ * software and associated documentation files (the "Software"), to deal in the Software *
+ * without restriction, including without limitation the rights to use, copy, modify,    *
+ * merge, publish, distribute, sublicense, and/or sell copies of the Software, and to    *
+ * permit persons to whom the Software is furnished to do so, subject to the following   *
+ * conditions:                                                                           *
+ *                                                                                       *
+ * The above copyright notice and this permission notice shall be included in all copies *
+ * or substantial portions of the Software.                                              *
+ *                                                                                       *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   *
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A         *
+ * PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT    *
+ * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF  *
+ * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE  *
+ * OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.                                         *
+ ****************************************************************************************/
+
+// open space includes
+#include <openspace/rendering/transferfunction.h>
+
+// ghoul includes
+#include <ghoul/opengl/texture.h>
+#include <ghoul/io/texture/texturereader.h>
+#include <ghoul/filesystem/filesystem.h>
+#include <ghoul/filesystem/cachemanager.h>
+
+#include <ghoul/logging/logmanager.h>
+#include <fstream>
+#include <string>
+#include <cstring>
+#include <iterator>
+
+#include <iostream>
+
+
+namespace {
+    const std::string _loggerCat = "TransferFunction";
+
+    ghoul::opengl::Texture::FilterMode filtermode = ghoul::opengl::Texture::FilterMode::Linear;
+    ghoul::opengl::Texture::WrappingMode wrappingmode = ghoul::opengl::Texture::WrappingMode::ClampToEdge;
+
+    bool hasExtension (std::string const &filepath, std::string const &extension)
+    {
+        std::string ending = "." + extension;
+        if (filepath.length() > ending.length()) {
+            return (0 == filepath.compare (filepath.length() - ending.length(), ending.length(), ending));
+        } else {
+            return false;
+        }
+    }
+
+}
+
+
+namespace openspace {
+TransferFunction::TransferFunction(const std::string& filepath, TfChangedCallback tfChangedCallback) : _filepath(filepath) {
+    setPath(filepath);
+    setCallback(tfChangedCallback);
+}
+
+void TransferFunction::setPath(const std::string& filepath) {
+    if (_file) {
+        _file = nullptr;
+    }
+    std::string f = absPath(filepath);
+    if (!FileSys.fileExists(f)) {
+        LERROR("Could not find transfer function file.");
+        _file = nullptr;
+    }
+    _filepath = f;
+    _file = std::make_unique<ghoul::filesystem::File>(filepath, ghoul::filesystem::File::RawPath::Yes);
+    _needsUpdate = true;
+    _file->setCallback([this](const ghoul::filesystem::File& file) {
+        _needsUpdate = true;
+    });
+}
+
+ghoul::opengl::Texture& TransferFunction::getTexture() {
+    ghoul_assert(_texture != nullptr, "Transfer function is null");
+    update();
+    return *_texture.get();
+}
+
+void TransferFunction::update() {
+    if (_needsUpdate) {
+        if (hasExtension(_filepath, "txt")) {
+            setTextureFromTxt();
+        } else {
+            setTextureFromImage();
+        }
+        _texture->uploadTexture();
+        _needsUpdate = false;
+        if (_tfChangedCallback) {
+            _tfChangedCallback(*this);
+        }
+    }
+}
+
+void TransferFunction::setCallback(TfChangedCallback callback) {
+    _tfChangedCallback = std::move(callback);
+}
+
+void TransferFunction::setTextureFromTxt() {
+    std::ifstream in;
+    in.open(_filepath.c_str());
+
+    if (!in.is_open()) {
+        LERROR("Could not open file " << _filepath);
+        return;
+    }
+
+    int width = 512;
+    float lower = 0.0f;
+    float upper = 1.0f;
+
+    std::vector<MappingKey> mappingKeys;
+
+    std::string line;
+
+    while (std::getline(in, line)) {
+        std::istringstream iss(line);
+        std::string key;
+        iss >> key;
+
+        if(key == "width") {
+            iss >> width;
+        } else if(key == "lower") {
+            iss >> lower;
+            lower = glm::clamp(lower, 0.0f, 1.0f);
+        } else if(key == "upper") {
+            iss >> upper;
+            upper = glm::clamp(upper, lower, 1.0f);
+        } else if(key == "mappingkey") {
+            float intensity;
+            glm::vec4 rgba = glm::vec4(0.0f);
+            iss >> intensity;
+            for(int i = 0; i < 4; ++i) {
+                iss >> rgba[i];
+            }
+            mappingKeys.push_back({intensity, rgba});
+        }
+    }
+    in.close();
+
+    if (mappingKeys.size() < 1) {
+        return;
+    }
+
+    if (mappingKeys.front().position > lower) {
+        mappingKeys.insert(mappingKeys.begin(), {lower,mappingKeys.front().color});
+    }
+
+    if (mappingKeys.back().position < upper) {
+        mappingKeys.push_back({upper,mappingKeys.back().color});
+    }
+
+    // allocate new float array with zeros
+    float* transferFunction = new float[width*4]();
+    for (int i = 0; i < 4*width; ++i) {
+        transferFunction[i] = 0.0f;
+    }
+
+    size_t lowerIndex = static_cast<size_t>(floorf(lower*static_cast<float>(width-1)));
+    size_t upperIndex = static_cast<size_t>(floorf(upper*static_cast<float>(width-1)));
+
+    auto prevKey = mappingKeys.begin();
+    auto currentKey = prevKey + 1;
+    auto lastKey = mappingKeys.end() -1;
+
+    for (size_t i=lowerIndex; i<=upperIndex; i++) {
+        float fpos = static_cast<float>(i)/static_cast<float>(width-1);
+        if (fpos > (*currentKey).position) {
+            prevKey = currentKey;
+            currentKey++;
+            if (currentKey == mappingKeys.end()) {
+                currentKey = lastKey;
+            }
+        }
+
+        float dist = fpos-(*prevKey).position;
+        float weight = dist/((*currentKey).position-(*prevKey).position);
+
+        for (size_t channel=0; channel<4; ++channel) {
+            size_t position = 4*i + channel;
+            // Interpolate linearly between prev and next mapping key
+
+            float value =
+                ((*prevKey).color[channel]*(1.f-weight) + (*currentKey).color[channel]*weight)/255.f;
+            transferFunction[position] = value;
+        }
+    }
+
+    // no need to deallocate transferFunction. Ownership is transferred to the Texture.
+    _texture = std::make_unique<ghoul::opengl::Texture>(transferFunction,
+                                            glm::size3_t(width,1,1),ghoul::opengl::Texture::Format::RGBA,
+                                            GL_RGBA, GL_FLOAT, filtermode, wrappingmode);
+}
+
+void TransferFunction::setTextureFromImage() {
+    _texture = ghoul::io::TextureReader::ref().loadTexture(_filepath);
+    _texture->setWrapping(wrappingmode);
+}
+
+glm::vec4 TransferFunction::sample(size_t offset) {
+    if (!_texture) return glm::vec4(0.0);
+
+    int nPixels = _texture->width();
+
+    // Clamp to range.
+    if (offset >= nPixels) offset = nPixels - 1;
+    if (offset < 0) offset = 0;
+
+    return _texture->texelAsFloat(offset);
+}
+
+size_t TransferFunction::width() {
+    update();
+    return _texture->width();
+}
+}