MedianFilter demonstrates use of compute shaders

modules/preprocessing/glsl/medianfilter.comp

+// ================================================================================================
+// 
+// This file is part of the CAMPVis Software Framework.
+// 
+// If not explicitly stated otherwise: Copyright (C) 2012-2015, all rights reserved,
+//      Christian Schulte zu Berge <christian.szb@in.tum.de>
+//      Chair for Computer Aided Medical Procedures
+//      Technische Universitaet Muenchen
+//      Boltzmannstr. 3, 85748 Garching b. Muenchen, Germany
+// 
+// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
+// 
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file 
+// except in compliance with the License. You may obtain a copy of the License at
+// 
+// http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software distributed under the 
+// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, 
+// either express or implied. See the License for the specific language governing permissions 
+// and limitations under the License.
+// 
+// ================================================================================================
+
+
+// ==================== Shader for median filtering. ====================
+// Expects the following dynamic defines:
+//
+// #define OUTPUT_TEXTURE_FORMAT
+//          The texture format of the output texture. Preferrably specified through Texture::calcMatchingWriteFormat()
+//
+// #define MEDIAN_WINDOW_X
+// #define MEDIAN_WINDOW_Y
+// #define MEDIAN_WINDOW_Z
+//          Size of the Median filter kernel in every dimension. All values to 1 performs results in a 3x3x3 filter window
+
+// The work group size can be overridden by dynamic defines
+#ifndef WORK_GROUP_SIZE_X
+#define WORK_GROUP_SIZE_X 1
+#endif
+#ifndef WORK_GROUP_SIZE_Y
+#define WORK_GROUP_SIZE_Y 1
+#endif
+#ifndef WORK_GROUP_SIZE_Z
+#define WORK_GROUP_SIZE_Z 1
+#endif
+
+#include "tools/texture3d.frag"
+
+uniform sampler3D _texture;
+uniform TextureParameters3D _textureParams;
+
+layout(OUTPUT_TEXTURE_FORMAT, binding = 0) uniform image3D _outputImage;
+layout(local_size_x = WORK_GROUP_SIZE_X, local_size_y = WORK_GROUP_SIZE_Y, local_size_z = WORK_GROUP_SIZE_Z) in;
+
+#define KERNEL_X (MEDIAN_WINDOW_X*2 + 1)
+#define KERNEL_Y (MEDIAN_WINDOW_Y*2 + 1)
+#define KERNEL_Z (MEDIAN_WINDOW_Z*2 + 1)
+#define KERNEL_S KERNEL_X*KERNEL_Y*KERNEL_Z
+
+float sorted[KERNEL_S];
+int end = 0;
+
+void sorted_insert2(float val) {
+    for(int i=0; i < end; i++) {
+        if(sorted[i].x < val.x) {
+            float tmp = val;
+            val = sorted[i];
+            sorted[i] = tmp;
+        }
+    }
+    sorted[end] = val;
+    end++;
+}
+
+void sorted_insert(float val) {
+    int i=end-1;
+    while( i >= 0) {
+        if(sorted[i].x < val.x) {
+            sorted[i+1] = sorted[i];
+        }
+        else {
+            break;
+        }
+        --i;
+    }
+    sorted[i+1] = val;
+    end++;
+}
+
+void sorted_insert3(float val) {
+    sorted[end] = val;
+    end++;
+}
+
+void main() {
+  // get index in global work group i.e x,y position
+  ivec3 pixel_coords = ivec3(gl_GlobalInvocationID.xyz);
+  
+  //vec4 pixel = texelFetch(_texture, pixel_coords, 0);
+  vec4 pixel = vec4(0);
+
+  const ivec3 wnd = ivec3(MEDIAN_WINDOW_X, MEDIAN_WINDOW_Y, MEDIAN_WINDOW_Z);
+
+  ivec3 lo = max(pixel_coords - wnd, ivec3(0));
+  ivec3 hi = min(pixel_coords + wnd, ivec3(_textureParams._size) - ivec3(1));
+
+  for(int x = lo.x; x <= hi.x; x++) {
+      for(int y = lo.y; y <= hi.y; y++) {
+          for(int z = lo.z; z <= hi.z; z++) {
+              sorted_insert(texelFetch(_texture, ivec3(x,y,z), 0).r);
+          }
+      }
+  }
+   
+  float result = sorted[end/2];
+
+  // output to a specific pixel in the image
+  imageStore(_outputImage, pixel_coords, vec4(result));
+}
\ No newline at end of file

modules/preprocessing/preprocessing.cmake

@@ -12,6 +12,7 @@ IF(${ModuleEnabled})
 	# Header files
 	FILE(GLOB ThisModHeaders RELATIVE ${ModulesDir}
 		modules/preprocessing/glsl/*.frag
+		modules/preprocessing/glsl/*.comp
 		modules/preprocessing/pipelines/*.h
 		modules/preprocessing/processors/*.h
 		modules/preprocessing/tools/*.h

modules/preprocessing/preprocessing.cpp

@@ -39,6 +39,7 @@
 #include "modules/preprocessing/processors/glstructuralsimilarity.h"
 #include "modules/preprocessing/processors/glvesselnessfilter.h"
 #include "modules/preprocessing/processors/gradientvolumegenerator.h"
+#include "modules/preprocessing/processors/medianfilter.h"
 
 namespace campvis {
 
@@ -60,5 +61,6 @@ namespace campvis {
     template class SmartProcessorRegistrar<GlStructuralSimilarity>;
     template class SmartProcessorRegistrar<GlVesselnessFilter>;
     template class SmartProcessorRegistrar<GradientVolumeGenerator>;
+    template class SmartProcessorRegistrar<MedianFilter>;
 
 }
\ No newline at end of file

modules/preprocessing/processors/medianfilter.cpp

+// ================================================================================================
+// 
+// This file is part of the CAMPVis Software Framework.
+// 
+// If not explicitly stated otherwise: Copyright (C) 2012-2015, all rights reserved,
+//      Christian Schulte zu Berge <christian.szb@in.tum.de>
+//      Chair for Computer Aided Medical Procedures
+//      Technische Universitaet Muenchen
+//      Boltzmannstr. 3, 85748 Garching b. Muenchen, Germany
+// 
+// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
+// 
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file 
+// except in compliance with the License. You may obtain a copy of the License at
+// 
+// http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software distributed under the 
+// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, 
+// either express or implied. See the License for the specific language governing permissions 
+// and limitations under the License.
+// 
+// ================================================================================================
+
+#include "medianfilter.h"
+
+#include "cgt/buffer.h"
+#include "cgt/imageunit.h"
+#include "cgt/logmanager.h"
+#include "cgt/shadermanager.h"
+#include "cgt/textureunit.h"
+#include "cgt/texture.h"
+
+#include "core/datastructures/imagedata.h"
+#include "core/datastructures/imagerepresentationgl.h"
+#include "core/datastructures/renderdata.h"
+
+#include "core/tools/quadrenderer.h"
+#include "core/tools/stringutils.h"
+
+namespace campvis {
+
+    const std::string MedianFilter::loggerCat_ = "CAMPVis.modules.preprocessing.MedianFilter";
+
+
+    MedianFilter::MedianFilter()
+        : AbstractProcessor()
+        , p_inputImage("InputImage", "Input Image", "", DataNameProperty::READ)
+        , p_outputImage("OutputImage", "Output Image", "GlGaussianFilter.out", DataNameProperty::WRITE)
+        , p_windowSize("WindowSize", "Window Size", cgt::ivec3(1, 1, 0), cgt::ivec3(0), cgt::ivec3(30))
+        , p_workgroupSize("WorkgroupSize", "Workgroup Size", cgt::ivec3(4), cgt::ivec3(1), cgt::ivec3(16))
+        , _shader(0)
+    {
+        addProperty(p_inputImage);
+        addProperty(p_outputImage);
+        addProperty(p_windowSize);
+        addProperty(p_workgroupSize);
+    }
+
+    MedianFilter::~MedianFilter() {
+
+    }
+
+    void MedianFilter::init() {
+        AbstractProcessor::init();
+   
+        LGL_ERROR;
+    }
+
+    void MedianFilter::deinit() {
+        ShdrMgr.dispose(_shader);
+        
+        AbstractProcessor::deinit();
+    }
+
+    void MedianFilter::updateResult(DataContainer& data) {
+        ImageRepresentationGL::ScopedRepresentation img(data, p_inputImage.getValue());
+
+        if (img != 0) {
+            if (img->getParent()->getDimensionality() > 1) {
+
+                generateShader(img);
+
+                if (_shader) {
+                    cgt::ivec3 size = img->getSize();
+                    cgt::ivec3 wgSize = p_workgroupSize.getValue();
+                    cgt::vec3 groups = cgt::ceil(cgt::vec3(size) / cgt::vec3(wgSize));
+
+                    cgt::TextureUnit inputUnit;
+                    inputUnit.activate();
+
+                    cgt::ImageUnit outputUnit;
+
+                    _shader->activate();
+
+                    // create texture for result
+                    std::unique_ptr<cgt::Texture> resultTexture(new cgt::Texture(img->getTexture()->getType(), size, img->getTexture()->getInternalFormat(), cgt::Texture::LINEAR));
+                    img->bind(_shader, inputUnit);
+                    LGL_ERROR;
+                    resultTexture->bindImage(outputUnit, GL_WRITE_ONLY);
+                    LGL_ERROR;
+                    _shader->setUniform("_outputImage", outputUnit.getUnitNumber());
+                    LGL_ERROR;
+                    
+                    LGL_ERROR;
+
+                    glDispatchCompute(groups.x, groups.y, groups.z);
+                    LGL_ERROR;
+
+                    _shader->deactivate();
+
+                    glFinish();
+
+                    // put resulting image into DataContainer
+                    std::unique_ptr<ImageData> id(new ImageData(img->getParent()->getDimensionality(), size, img->getParent()->getNumChannels()));
+                    ImageRepresentationGL::create(id.get(), resultTexture.release());
+                    id->setMappingInformation(img->getParent()->getMappingInformation());
+                    data.addData(p_outputImage.getValue(), id.release());
+
+                    cgt::TextureUnit::setZeroUnit();
+                    LGL_ERROR;
+                }
+            }
+            else {
+                LERROR("Supports only 2D and 3D Median Filtering.");
+            }
+        }
+        else {
+            LDEBUG("No suitable input image found.");
+        }
+    }
+
+    void MedianFilter::generateShader(const ImageRepresentationGL::ScopedRepresentation & img)
+    {
+        auto wnd = p_windowSize.getValue();
+        auto wg = p_workgroupSize.getValue();
+        std::stringstream ss;
+        ss << "#define MEDIAN_WINDOW_X " << wnd.x << std::endl;
+        ss << "#define MEDIAN_WINDOW_Y " << wnd.y << std::endl;
+        ss << "#define MEDIAN_WINDOW_Z " << wnd.z << std::endl;
+        ss << "#define WORK_GROUP_SIZE_X " << wg.x << std::endl;
+        ss << "#define WORK_GROUP_SIZE_Y " << wg.y << std::endl;
+        ss << "#define WORK_GROUP_SIZE_Z " << wg.z << std::endl;
+        ss << "#define OUTPUT_TEXTURE_FORMAT " << cgt::Texture::calcMatchingWriteFormat(img->getTexture()->getInternalFormat()) << std::endl;
+
+        // very simple caching to eliminate the glsl compiler as a bottleneck
+        static std::string headers;
+        if (!_shader || headers != ss.str()) {
+            headers = ss.str();
+            _shader = ShdrMgr.loadCompute("modules/preprocessing/glsl/medianfilter.comp", headers);
+        }
+    }
+}

modules/preprocessing/processors/medianfilter.h

+// ================================================================================================
+// 
+// This file is part of the CAMPVis Software Framework.
+// 
+// If not explicitly stated otherwise: Copyright (C) 2012-2015, all rights reserved,
+//      Christian Schulte zu Berge <christian.szb@in.tum.de>
+//      Chair for Computer Aided Medical Procedures
+//      Technische Universitaet Muenchen
+//      Boltzmannstr. 3, 85748 Garching b. Muenchen, Germany
+// 
+// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
+// 
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file 
+// except in compliance with the License. You may obtain a copy of the License at
+// 
+// http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software distributed under the 
+// License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, 
+// either express or implied. See the License for the specific language governing permissions 
+// and limitations under the License.
+// 
+// ================================================================================================
+
+#ifndef MEDIANFILTER_H__
+#define MEDIANFILTER_H__
+
+#include <string>
+
+#include "core/datastructures/imagerepresentationgl.h"
+#include "core/pipeline/abstractprocessordecorator.h"
+#include "core/pipeline/visualizationprocessor.h"
+#include "core/properties/datanameproperty.h"
+#include "core/properties/floatingpointproperty.h"
+#include "core/properties/optionproperty.h"
+
+#include "modules/modulesapi.h"
+
+namespace cgt {
+    class BufferObject;
+    class Shader;
+}
+
+
+namespace campvis {
+    /**
+     * Performs a median filtering on the input image using OpenGL Compute Shaders.
+     */
+    class CAMPVIS_MODULES_API MedianFilter : public AbstractProcessor {
+    public:
+        /**
+         * Constructs a new GlGaussianFilter Processor
+         **/
+        MedianFilter();
+
+        /**
+         * Destructor
+         **/
+        virtual ~MedianFilter();
+
+        /// \see AbstractProcessor::init
+        virtual void init();
+        /// \see AbstractProcessor::deinit
+        virtual void deinit();
+        
+        /** 
+         * To be used in ProcessorFactory static methods
+         */
+        static const std::string getId() { return "MedianFilter"; };
+        /// \see AbstractProcessor::getName()
+        virtual const std::string getName() const { return getId(); };
+        /// \see AbstractProcessor::getDescription()
+        virtual const std::string getDescription() const { return "Performs a median filtering on the input image using OpenGL Compute Shaders."; };
+        /// \see AbstractProcessor::getAuthor()
+        virtual const std::string getAuthor() const { return "Jakob Weiss <jakob.weiss@tum.de>"; };
+        /// \see AbstractProcessor::getProcessorState()
+        virtual ProcessorState getProcessorState() const { return AbstractProcessor::TESTING; };
+
+        DataNameProperty p_inputImage;      ///< ID for input volume
+        DataNameProperty p_outputImage;     ///< ID for output gradient volume
+
+        IVec3Property p_windowSize;              ///< Size for specifying filtering window
+        IVec3Property p_workgroupSize;  ///< Size of the workgroup
+
+    protected:
+        /// \see AbstractProcessor::updateResult
+        virtual void updateResult(DataContainer& dataContainer) override;
+
+        void generateShader(const ImageRepresentationGL::ScopedRepresentation& img);
+
+        cgt::Shader* _shader;             ///< Shader for performing median filtering
+
+        static const std::string loggerCat_;
+    };
+}
+
+#endif // MEDIANFILTER_H__