Fixed computation of SSIM in GlStructuralSimilarity processor. Extended...

Fixed computation of SSIM in GlStructuralSimilarity processor. Extended SsimDemo pipeline to compute batch-wise SSIM of an entire US sequence.

modules/preprocessing/glsl/glstructuralsimilarity.frag

@@ -33,36 +33,37 @@ uniform TextureParameters2D _image1Params;
 uniform sampler2D _image2;
 uniform TextureParameters2D _image2Params;
 
-#define WINDOW_SIZE 4
-#define WINDOW_SIZE_SQ 16
+#define HALF_WINDOW_SIZE 4
+#define WINDOW_SIZE 9
+#define WINDOW_SIZE_SQ 81
 
 
 void main() {
     vec4 toReturn = vec4(0.0);
 
-    const float L = 255;
+    const float L = 128;
     const float c1 = (0.01 * L) * (0.01 * L);
     const float c2 = (0.03 * L) * (0.03 * L);
 
     vec4 xMean = vec4(0);
     vec4 yMean = vec4(0);
-    vec4 xM2 = vec4(0);
-    vec4 yM2 = vec4(0);
-    vec4 cov = vec4(0);
+    vec4 xM2   = vec4(0);
+    vec4 yM2   = vec4(0);
+    vec4 cov   = vec4(0);
     float n = 0;
 
-    for (int i = -WINDOW_SIZE_SQ; i <= WINDOW_SIZE_SQ; ++i) {
-        vec2 iOffset = vec2(i % WINDOW_SIZE, i / WINDOW_SIZE) * _image1Params._sizeRCP;
-        vec4 xi = texture(_image1, ex_TexCoord.xy + iOffset);
-        vec4 yi = texture(_image2, ex_TexCoord.xy + iOffset);
+    vec4 mean1 = vec4(0);
+    vec4 mean2 = vec4(0);
+    vec4 M12   = vec4(0);
+    float NN = WINDOW_SIZE_SQ;
 
-        for (int j = i+1; j <= WINDOW_SIZE_SQ; ++j) {
-            vec2 jOffset = vec2(j % WINDOW_SIZE, j / WINDOW_SIZE) * _image1Params._sizeRCP;
-            vec4 xj = texture(_image1, ex_TexCoord.xy + jOffset);
-            vec4 yj = texture(_image2, ex_TexCoord.xy + jOffset);
+    vec2 shift = _image1Params._sizeRCP * HALF_WINDOW_SIZE;
+    float mult = _image1Params._sizeRCP;
 
-            cov += (xi - xj) * (yi - yj);
-        }
+    for (int i = 0; i <= WINDOW_SIZE_SQ; ++i) {
+        vec2 iOffset = (vec2(i % WINDOW_SIZE, i / WINDOW_SIZE) * mult) - shift;
+        vec4 xi = texture(_image1, ex_TexCoord.xy + iOffset);
+        vec4 yi = texture(_image2, ex_TexCoord.xy + iOffset);
 
         // compute mean and variance with Knuth's algorithm 
         // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
@@ -75,11 +76,22 @@ void main() {
         vec4 yDelta = yi - yMean;
         yMean += yDelta/n;
         yM2 += yDelta * (yi - yMean);
+
+
+
+        // compute the covariance
+        vec4 delta1 = (xi - mean1) / n;
+        mean1 += delta1;
+        vec4 delta2 = (yi - mean2) / n;
+        mean2 += delta2;
+        M12 += n * delta1 * delta2 - M12 / n;
+
     }
 
+    cov = NN / (NN - 1) * M12;
     vec4 xVar = xM2 / (n-1);
     vec4 yVar = yM2 / (n-1);
     vec4 ssim = ((2*xMean*yMean + c1) * (2*cov + c2)) / ((xMean * xMean + yMean * yMean + c1) * (xVar + yVar + c2));
 
-    out_Color = vec4(ssim);
+    out_Color = min(vec4(1.0), vec4(ssim));
 }

modules/preprocessing/pipelines/ssimdemo.cpp

@@ -24,20 +24,40 @@
 
 #include "ssimdemo.h"
 
+#include "cgt/glcontextmanager.h"
 #include "cgt/shadermanager.h"
+#include "core/datastructures/imagerepresentationgl.h"
 
+#include <iomanip>
+#include <fstream>
+#include <sstream>
 
 namespace campvis {
 
     SsimDemo::SsimDemo(DataContainer& dc)
         : AutoEvaluationPipeline(dc, getId())
-        , _imageReader()
-        , _gaussian(&_canvasSize)
+        , _imageReader1()
+        , _imageReader2()
         , _ssim(&_canvasSize)
+        , _fanRenderer(&_canvasSize)
+        , _sumReduction(nullptr)
+        , _minReduction(nullptr)
+        , _currentlyBatchProcessing(false)
+        , p_sourcePath1("SourcePath1", "Source Path 1", "Z:/cm_stuff/cmCuda_512_1000", StringProperty::DIRECTORY)
+        , p_sourcePath2("SourcePath2", "Source Path 2", "Z:/cm_stuff/cmCuda_128_1000", StringProperty::DIRECTORY)
+        , p_range("Range", "range", cgt::ivec2(0, 399), cgt::ivec2(0), cgt::ivec2(1000))
+        , p_execute("Excecute", "Execute Batch Process")
     {
-        addProcessor(&_imageReader);
-        addProcessor(&_gaussian);
+        addProcessor(&_imageReader1);
+        addProcessor(&_imageReader2);
         addProcessor(&_ssim);
+        addProcessor(&_fanRenderer);
+
+        addProperty(p_sourcePath1);
+        addProperty(p_sourcePath2);
+        addProperty(p_range);
+        addProperty(p_execute);
+        p_execute.s_clicked.connect(this, &SsimDemo::executeBatchProcess);
     }
 
     SsimDemo::~SsimDemo() {
@@ -46,15 +66,99 @@ namespace campvis {
     void SsimDemo::init() {
         AutoEvaluationPipeline::init();
 
-        _ssim.p_outputImage.setValue("result");
-        _renderTargetID.setValue("result");
+        _sumReduction = new GlReduction(GlReduction::PLUS);
+        _minReduction = new GlReduction(GlReduction::MIN);
+
+        _renderTargetID.setValue("us.fan");
+
+        _imageReader1.p_url.setValue("Z:/cm_stuff/cmCuda_512_1000/export0026.bmp");
+        _imageReader1.p_importType.setValue(1);
+        _imageReader1.p_targetImageID.setValue("reader1.output");
+
+        _imageReader2.p_url.setValue("Z:/cm_stuff/cmCuda_128_1000/export0026.bmp");
+        _imageReader2.p_importType.setValue(1);
+        _imageReader2.p_targetImageID.setValue("reader2.output");
+
+        _ssim.p_inputImage1.setValue("reader1.output");
+        _ssim.p_inputImage2.setValue("reader2.output");
+        _ssim.s_validated.connect(this, &SsimDemo::onProcessorValidated);
+        _ssim.p_outputImage.setValue("ssim");
+
+        _fanRenderer.p_inputImage.setValue("ssim");
+        _fanRenderer.p_renderTargetID.setValue("us.fan");
 
-        _imageReader.p_url.setValue(ShdrMgr.completePath("/modules/"));
-        _imageReader.p_targetImageID.setValue("reader.output");
-        _imageReader.p_targetImageID.addSharedProperty(&_gaussian.p_inputImage);
+    }
+
+    void SsimDemo::deinit() {
+        delete _minReduction;
+        delete _sumReduction;
+        AutoEvaluationPipeline::deinit();
+    }
 
-        _gaussian.p_outputImage.setValue("blurred");
+    void SsimDemo::onProcessorValidated(AbstractProcessor* p) {
+        if (_currentlyBatchProcessing)
+            return;
+
+        if (p == &_ssim) {
+            ImageRepresentationGL::ScopedRepresentation ssim(getDataContainer(), _ssim.p_outputImage.getValue());
+            if (ssim) {
+                cgt::GLContextScopedLock lock(this->_canvas);
+                auto sums = _sumReduction->reduce(ssim->getTexture());
+                auto mins = _minReduction->reduce(ssim->getTexture());
+
+                LINFO("Structured Similarity, Average: " << (sums[0] / ssim->getNumElements()) << ", Minimum: " << mins[0]);
+            }
+        }
     }
 
+    void SsimDemo::executeBatchProcess() {
+        if (p_range.getValue().x > p_range.getValue().y)
+            return;
+
+        _currentlyBatchProcessing = true;
+        cgt::GLContextScopedLock lock(_canvas);
+
+        float sumsum = 0.f;
+        float minmin = 1.f;
+
+        std::string foo = p_sourcePath1.getValue().substr(p_sourcePath1.getValue().find_last_of("/\\") + 1);
+        std::string bar = p_sourcePath2.getValue().substr(p_sourcePath2.getValue().find_last_of("/\\") + 1);
+
+        LINFO("Comparing " << p_sourcePath1.getValue() << " to " << p_sourcePath2.getValue() << ":");
+
+        std::ofstream csvFile("C:\\temp\\" + foo + "---" + bar + ".csv");
+        csvFile << "index, AvgSSIM, MinSSIM\n";
+
+        for (int i = p_range.getValue().x; i < p_range.getValue().y; ++i) {
+            // set up processors:
+            std::stringstream ss;
+            ss << "export" << std::setfill('0') << std::setw(4) << i << ".bmp";
+            std::string fileName = ss.str();
+
+            // read image
+            _imageReader1.p_url.setValue(p_sourcePath1.getValue() + "\\" + fileName);
+            _imageReader2.p_url.setValue(p_sourcePath2.getValue() + "\\" + fileName);
+            forceExecuteProcessor(&_imageReader1);
+            forceExecuteProcessor(&_imageReader2);
+            forceExecuteProcessor(&_ssim);
+
+            ImageRepresentationGL::ScopedRepresentation ssim(getDataContainer(), _ssim.p_outputImage.getValue());
+            if (ssim) {
+                auto sums = _sumReduction->reduce(ssim->getTexture());
+                auto mins = _minReduction->reduce(ssim->getTexture());
+
+                sumsum += sums[0] / ssim->getNumElements();
+                minmin = std::min(minmin, mins[0]);
+
+                csvFile << i << ", " << (sums[0] / ssim->getNumElements()) << ", " << mins[0] << "\n";
+            }
+
+            if (i % 100 == 0) 
+                LINFO(i);
+        }
+
+        LINFO("Structured Similarity, Averaged average: " << (sumsum / (p_range.getValue().y - p_range.getValue().x)) << ", Minimum of the Minima: " << minmin);
+        _currentlyBatchProcessing = false;
+    }
 
 }
\ No newline at end of file

modules/preprocessing/pipelines/ssimdemo.h

@@ -25,12 +25,14 @@
 #ifndef SSIMDEMO_H__
 #define SSIMDEMO_H__
 
+#include "core/tools/glreduction.h"
 #include "core/pipeline/autoevaluationpipeline.h"
 
 #include "modules/modulesapi.h"
 #include "modules/devil/processors/devilimagereader.h"
 #include "modules/preprocessing/processors/glgaussianfilter.h"
 #include "modules/preprocessing/processors/glstructuralsimilarity.h"
+#include "modules/vis/processors/usfanrenderer.h"
 
 namespace campvis {
     class CAMPVIS_MODULES_API SsimDemo : public AutoEvaluationPipeline {
@@ -47,16 +49,31 @@ namespace campvis {
          **/
         virtual ~SsimDemo();
 
-        /// \see AutoEvaluationPipeline::init()
-        virtual void init();
+        virtual void init() override;
+        virtual void deinit() override;
 
 
         static const std::string getId() { return "SsimDemo"; };
 
     protected:
-        DevilImageReader _imageReader;
-        GlGaussianFilter _gaussian;
+        void onProcessorValidated(AbstractProcessor* p);
+
+        void executeBatchProcess();
+
+        DevilImageReader _imageReader1;
+        DevilImageReader _imageReader2;
         GlStructuralSimilarity _ssim;
+        UsFanRenderer _fanRenderer;
+
+        GlReduction* _sumReduction;
+        GlReduction* _minReduction;
+        bool _currentlyBatchProcessing;
+
+        StringProperty p_sourcePath1;
+        StringProperty p_sourcePath2;
+
+        IVec2Property p_range;                          ///< Range for image iteration
+        ButtonProperty p_execute;                       ///< Button to start the batch process
     };
 }
 

modules/preprocessing/processors/glstructuralsimilarity.cpp

@@ -58,7 +58,7 @@ namespace campvis {
     void GlStructuralSimilarity::init() {
         VisualizationProcessor::init();
 
-        _shader2D = ShdrMgr.load("core/glsl/passthrough.vert", "modules/preprocessing/glsl/glstructuralsimilarity.frag", "");
+        _shader2D = ShdrMgr.loadWithCustomGlslVersion("core/glsl/passthrough.vert", "", "modules/preprocessing/glsl/glstructuralsimilarity.frag", "", "400");
     }
 
     void GlStructuralSimilarity::deinit() {
@@ -77,7 +77,7 @@ namespace campvis {
             img1Unit.activate();
 
             // create texture for result
-            cgt::Texture* resultTexture = new cgt::Texture(GL_TEXTURE_2D, originalSize, img1->getTexture()->getInternalFormat(), cgt::Texture::LINEAR);
+            cgt::Texture* resultTexture = new cgt::Texture(GL_TEXTURE_2D, originalSize, GL_R32F, cgt::Texture::LINEAR);
             LGL_ERROR;
 
             // activate shader and bind textures
@@ -99,7 +99,7 @@ namespace campvis {
             _shader2D->deactivate();
 
             // put resulting image into DataContainer
-            ImageData* id = new ImageData(img1->getParent()->getDimensionality(), originalSize, img1->getParent()->getNumChannels());
+            ImageData* id = new ImageData(img1->getParent()->getDimensionality(), originalSize, 1);
             ImageRepresentationGL::create(id, resultTexture);
             id->setMappingInformation(img1->getParent()->getMappingInformation());
             data.addData(p_outputImage.getValue(), id);