Further work on registration module:

 * added computation of difference image to SimilarityMeasure
 * fixed possible out of bounds access in datacontainerinspectorwidget.cpp
 * added assertions to ImageRepresentationGL constructor

application/gui/datacontainerinspectorwidget.cpp

@@ -315,7 +315,7 @@ namespace campvis {
         size_t index = 0;
         size_t remainder = 0;
 
-        while (numBytes > 1024 && index < 5) {
+        while (numBytes > 1024 && index < 4) {
             remainder = numBytes % 1024;
             numBytes /= 1024;
             ++index;

core/datastructures/imagerepresentationgl.cpp

@@ -62,6 +62,8 @@ namespace campvis {
         , _texture(texture)
     {
         tgtAssert(texture != 0, "Given texture must not be 0.");
+        tgtAssert(parent->getDimensionality() >= 3 || texture->getDimensions().z == 1, "Dimensionality of Parent and texture mismatch!");
+        tgtAssert(parent->getDimensionality() >= 2 || texture->getDimensions().y == 1, "Dimensionality of Parent and texture mismatch!");
     }
 
     ImageRepresentationGL::ImageRepresentationGL(ImageData* parent, const WeaklyTypedPointer& wtp) 

modules/registration/glsl/differenceimage.frag

+// ================================================================================================
+// 
+// This file is part of the CAMPVis Software Framework.
+// 
+// If not explicitly stated otherwise: Copyright (C) 2012, all rights reserved,
+//      Christian Schulte zu Berge <christian.szb@in.tum.de>
+//      Chair for Computer Aided Medical Procedures
+//      Technische Universität München
+//      Boltzmannstr. 3, 85748 Garching b. München, Germany
+// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
+// 
+// The licensing of this softare is not yet resolved. Until then, redistribution in source or
+// binary forms outside the CAMP chair is not permitted, unless explicitly stated in legal form.
+// However, the names of the original authors and the above copyright notice must retain in its
+// original state in any case.
+// 
+// Legal disclaimer provided by the BSD license:
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
+// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY 
+// AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR 
+// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
+// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
+// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
+// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+// 
+// ================================================================================================
+
+in vec3 ex_TexCoord;
+out vec4 out_Color;
+
+#include "tools/texture3d.frag"
+
+uniform sampler3D _referenceTexture;
+uniform TextureParameters3D _referenceTextureParams;
+
+uniform sampler3D _movingTexture;
+uniform TextureParameters3D _movingTextureParams;
+
+uniform float _zTex;
+uniform mat4 _registrationInverse;
+
+void main() {
+    // compute lookup coordinates
+    vec3 referenceLookupTexCoord = vec3(ex_TexCoord.xy, _zTex);
+    vec4 movingLookupTexCoord = _registrationInverse * vec4(referenceLookupTexCoord, 1.0);
+
+    // fetch texels
+    float referenceValue = texture(_referenceTexture, referenceLookupTexCoord).a;
+    float movingValue = texture(_movingTexture, movingLookupTexCoord.xyz).a;
+
+    // compute differences
+    float difference = referenceValue - movingValue;
+    float sad = abs(difference);
+    //float ssd = difference * difference;
+
+    // write output color
+    out_Color = vec4(sad, sad, sad, sad);
+}

modules/registration/pipelines/reducertest.cpp

@@ -41,18 +41,17 @@ namespace campvis {
 
     ReducerTest::ReducerTest(DataContainer* dc)
         : AutoEvaluationPipeline(dc)
-        , _imageReader()
-        , _sliceExtractor(&_canvasSize)
+        , _referenceReader()
+        , _movingReader()
         , _sm()
-        , _wheelHandler(&_sliceExtractor.p_zSliceNumber)
-        , _tfWindowingHandler(&_sliceExtractor.p_transferFunction)
+        , _ve(&_canvasSize)
     {
-        addProcessor(&_imageReader);
-        addProcessor(&_sliceExtractor);
+        addProcessor(&_referenceReader);
+        addProcessor(&_movingReader);
         addProcessor(&_sm);
+        addProcessor(&_ve);
 
-        addEventListenerToBack(&_wheelHandler);
-        addEventListenerToBack(&_tfWindowingHandler);
+        addEventListenerToBack(&_ve);
     }
 
     ReducerTest::~ReducerTest() {
@@ -61,54 +60,23 @@ namespace campvis {
     void ReducerTest::init() {
         AutoEvaluationPipeline::init();
 
-        _imageReader.p_url.setValue("D:\\Medical Data\\smallHeart.mhd");
-        _imageReader.p_targetImageID.setValue("reader.output");
-        _imageReader.p_targetImageID.addSharedProperty(&_sliceExtractor.p_sourceImageID);
-        _imageReader.p_targetImageID.addSharedProperty(&_sm.p_referenceId);
-        _imageReader.p_targetImageID.addSharedProperty(&_sm.p_movingId);
-        _imageReader.s_validated.connect(this, &ReducerTest::onProcessorValidated);
+        _referenceReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S3/Volume_0.mhd");
+        _referenceReader.p_targetImageID.setValue("Reference Image");
+        _referenceReader.p_targetImageID.addSharedProperty(&_sm.p_referenceId);
 
-        _sliceExtractor.p_xSliceNumber.setValue(0);
-        _sliceExtractor.s_validated.connect(this, &ReducerTest::onProcessorValidated);
+        _movingReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S3/Volume_1.mhd");
+        _movingReader.p_targetImageID.setValue("Moving Image");
+        _movingReader.p_targetImageID.addSharedProperty(&_sm.p_movingId);
 
-        // TODO: replace this hardcoded domain by automatically determined from image min/max values
-        Geometry1DTransferFunction* tf = new Geometry1DTransferFunction(128, tgt::vec2(0.f, .08f));
-        tf->addGeometry(TFGeometry1D::createQuad(tgt::vec2(0.f, 1.f), tgt::col4(0, 0, 0, 0), tgt::col4(255, 255, 255, 255)));
-        _sliceExtractor.p_transferFunction.replaceTF(tf);
 
-        _renderTargetID.setValue("renderTarget");
-        _renderTargetID.addSharedProperty(&(_sliceExtractor.p_targetImageID));
-    }
+        _sm.p_differenceImageId.addSharedProperty(&_ve.p_inputVolume);
 
-    void ReducerTest::keyEvent(tgt::KeyEvent* e) {
-        if (e->pressed()) {
-            switch (e->keyCode()) {
-                case tgt::KeyEvent::K_UP:
-                    _sliceExtractor.p_xSliceNumber.increment();
-                    break;
-                case tgt::KeyEvent::K_DOWN:
-                    _sliceExtractor.p_xSliceNumber.decrement();
-                    break;
-                default:
-                    break;
-            }
-        }
+        _ve.p_outputImage.setValue("renderTarget");
+        _renderTargetID.setValue("renderTarget");
     }
 
     void ReducerTest::onProcessorValidated(AbstractProcessor* processor) {
-        if (processor == &_imageReader) {
-            ScopedTypedData<ImageData> img(*_data, _imageReader.p_targetImageID.getValue());
-            if (img != 0) {
-                _sliceExtractor.p_transferFunction.getTF()->setImageHandle(img.getDataHandle());
-            }
-        }
-        if (processor == &_sliceExtractor) {
-            ScopedTypedData<RenderData> rd(*_data, _sliceExtractor.p_targetImageID.getValue());
-            if (rd != 0 && rd->getColorTexture() != 0) {
-//                 GlReduction reducer;
-//                 reducer.reduce(rd->getColorTexture());
-            }
-        }
+
     }
 
 }

modules/registration/pipelines/reducertest.h

@@ -35,7 +35,7 @@
 #include "core/eventhandlers/transfuncwindowingeventlistener.h"
 #include "core/pipeline/autoevaluationpipeline.h"
 #include "modules/io/processors/mhdimagereader.h"
-#include "modules/vis/processors/sliceextractor.h"
+#include "modules/vis/processors/volumeexplorer.h"
 #include "modules/preprocessing/processors/gradientvolumegenerator.h"
 #include "modules/preprocessing/processors/lhhistogram.h"
 
@@ -62,8 +62,6 @@ namespace campvis {
         static const std::string getId() { return "ReducerTest"; };
 
 
-        virtual void keyEvent(tgt::KeyEvent* e);
-
     protected:
         /**
          * Slot getting called when one of the observed processors got validated.
@@ -72,12 +70,11 @@ namespace campvis {
          */
         virtual void onProcessorValidated(AbstractProcessor* processor);
 
-        MhdImageReader _imageReader;
-        SliceExtractor _sliceExtractor;
+        MhdImageReader _referenceReader;
+        MhdImageReader _movingReader;
         SimilarityMeasure _sm;
 
-        MWheelToNumericPropertyEventListener _wheelHandler;
-        TransFuncWindowingEventListener _tfWindowingHandler;
+        VolumeExplorer _ve;
 
     };
 

modules/registration/processors/similaritymeasure.cpp

@@ -31,6 +31,7 @@
 #include "tgt/logmanager.h"
 #include "tgt/shadermanager.h"
 #include "tgt/textureunit.h"
+#include "tgt/openglgarbagecollector.h"
 
 #include "core/datastructures/facegeometry.h"
 #include "core/datastructures/imagedata.h"
@@ -62,7 +63,9 @@ namespace campvis {
         , p_computeSimilarity("ComputeSimilarity", "Compute Similarity")
         , p_optimizer("Optimizer", "Optimizer", optimizers, 3)
         , p_performOptimization("PerformOptimization", "Perform Optimization", AbstractProcessor::INVALID_RESULT | PERFORM_OPTIMIZATION)
-        , _shader(0)
+        , p_differenceImageId("DifferenceImageId", "Difference Image", "difference", DataNameProperty::WRITE, AbstractProcessor::VALID)
+        , p_computeDifferenceImage("ComputeDifferenceImage", "Compute Difference Image", AbstractProcessor::INVALID_RESULT | COMPUTE_DIFFERENCE_IMAGE)
+        , _costFunctionShader(0)
         , _glr(0)
     {
         addProperty(&p_referenceId);
@@ -70,6 +73,8 @@ namespace campvis {
         addProperty(&p_translation);
         addProperty(&p_rotation);
         addProperty(&p_computeSimilarity);
+        addProperty(&p_differenceImageId);
+        addProperty(&p_computeDifferenceImage);
         addProperty(&p_optimizer);
         addProperty(&p_performOptimization);
 
@@ -82,16 +87,20 @@ namespace campvis {
 
     void SimilarityMeasure::init() {
         VisualizationProcessor::init();
-        _shader = ShdrMgr.loadSeparate("core/glsl/passthrough.vert", "modules/registration/glsl/similaritymeasure.frag", "", false);
-        _shader->setAttributeLocation(0, "in_Position");
-        _shader->setAttributeLocation(1, "in_TexCoord");
+        _costFunctionShader = ShdrMgr.loadSeparate("core/glsl/passthrough.vert", "modules/registration/glsl/similaritymeasure.frag", "", false);
+        _costFunctionShader->setAttributeLocation(0, "in_Position");
+        _costFunctionShader->setAttributeLocation(1, "in_TexCoord");
+
+        _differenceShader = ShdrMgr.loadSeparate("core/glsl/passthrough.vert", "modules/registration/glsl/differenceimage.frag", "", false);
+        _differenceShader->setAttributeLocation(0, "in_Position");
+        _differenceShader->setAttributeLocation(1, "in_TexCoord");
 
         _glr = new GlReduction();
     }
 
     void SimilarityMeasure::deinit() {
         VisualizationProcessor::deinit();
-        ShdrMgr.dispose(_shader);
+        ShdrMgr.dispose(_costFunctionShader);
 
         delete _glr;
         _glr = 0;
@@ -108,6 +117,9 @@ namespace campvis {
 
             float similarity = computeSimilarity(referenceImage, movingImage, p_translation.getValue(), p_rotation.getValue());
             LDEBUG("Similarity Measure: " << similarity);
+
+            if (getInvalidationLevel() & COMPUTE_DIFFERENCE_IMAGE) 
+                generateDifferenceImage(&data, referenceImage, movingImage, p_translation.getValue(), p_rotation.getValue());
         }
         else {
             LERROR("No suitable input image found.");
@@ -189,14 +201,11 @@ namespace campvis {
         LGL_ERROR;
 
         // bind input images
-        _shader->activate();
-        referenceImage->bind(_shader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
-        movingImage->bind(_shader, movingUnit, "_movingTexture", "_movingTextureParams");
+        _costFunctionShader->activate();
+        referenceImage->bind(_costFunctionShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
+        movingImage->bind(_costFunctionShader, movingUnit, "_movingTexture", "_movingTextureParams");
 
-        tgt::mat4 registrationMatrix = tgt::mat4::createTranslation(translation) 
-            * tgt::mat4::createRotationZ(rotation.z)
-            * tgt::mat4::createRotationY(rotation.y)
-            * tgt::mat4::createRotationX(rotation.x);
+        tgt::mat4 registrationMatrix = tgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
 
         const tgt::mat4& w2t = movingImage->getParent()->getMappingInformation().getWorldToTextureMatrix();
         const tgt::mat4& t2w = referenceImage->getParent()->getMappingInformation().getTextureToWorldMatrix();
@@ -207,10 +216,9 @@ namespace campvis {
             tgtAssert(false, "Could not invert registration matrix. This should not happen!");
 
         // render quad to compute similarity measure by shader
-        //_shader->setUniform("_registrationMatrix", registrationMatrix);
-        _shader->setUniform("_registrationInverse", registrationInverse);
+        _costFunctionShader->setUniform("_registrationInverse", registrationInverse);
         QuadRdr.renderQuad();
-        _shader->deactivate();
+        _costFunctionShader->deactivate();
 
         // detach texture and reduce it
         //data.addData("All glory to the HYPNOTOAD!", new RenderData(_fbo));
@@ -237,11 +245,88 @@ namespace campvis {
         tgt::vec3 translation(x[0], x[1], x[2]);
         tgt::vec3 rotation(x[3], x[4], x[5]);
         float similarity = mfd->_object->computeSimilarity(mfd->_reference, mfd->_moving, translation, rotation);
+        GLGC.deleteGarbage();
 
         LDEBUG(translation << rotation << " : " << similarity);
 
         return similarity;
     }
 
+    tgt::mat4 SimilarityMeasure::euleranglesToMat4(const tgt::vec3& eulerAngles) {
+	    float sinX = sin(eulerAngles.x);
+	    float cosX = cos(eulerAngles.x);
+	    float sinY = sin(eulerAngles.y);
+	    float cosY = cos(eulerAngles.y);
+	    float sinZ = sin(eulerAngles.z);
+	    float cosZ = cos(eulerAngles.z);
+
+        tgt::mat4 toReturn(cosY * cosZ,   cosZ * sinX * sinY - cosX * sinZ,   sinX * sinZ + cosX * cosZ * sinY,   0.f,
+	                       cosY * sinZ,   sinX * sinY * sinZ + cosX * cosZ,   cosX * sinY * sinZ - cosZ * sinX,   0.f,
+	                       (-1) * sinY,   cosY * sinX,                        cosX * cosY,                        0.f,
+                           0.f,           0.f,                                0.f,                                1.f);
+        return toReturn;
+    }
+
+    void SimilarityMeasure::generateDifferenceImage(DataContainer* dc, const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const tgt::vec3& translation, const tgt::vec3& rotation) {
+        tgtAssert(dc != 0, "DataContainer must not be 0.");
+        tgtAssert(referenceImage != 0, "Reference Image must not be 0.");
+        tgtAssert(movingImage != 0, "Moving Image must not be 0.");
+
+        const tgt::svec3& size = referenceImage->getSize();
+        tgt::ivec2 viewportSize = size.xy();
+
+        // reserve texture units
+        tgt::TextureUnit referenceUnit, movingUnit;
+        referenceUnit.activate();
+
+        // create temporary texture for result
+        tgt::Texture* similarityTex = new tgt::Texture(0, tgt::ivec3(size), GL_ALPHA, GL_ALPHA32F_ARB, GL_FLOAT, tgt::Texture::NEAREST);
+        similarityTex->uploadTexture();
+        similarityTex->setWrapping(tgt::Texture::CLAMP);
+
+        // activate FBO and attach texture
+        _fbo->activate();
+        glViewport(0, 0, static_cast<GLsizei>(viewportSize.x), static_cast<GLsizei>(viewportSize.y));
+        LGL_ERROR;
+
+        // bind input images
+        _differenceShader->activate();
+        referenceImage->bind(_differenceShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
+        movingImage->bind(_differenceShader, movingUnit, "_movingTexture", "_movingTextureParams");
+
+        tgt::mat4 registrationMatrix = tgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
+
+        const tgt::mat4& w2t = movingImage->getParent()->getMappingInformation().getWorldToTextureMatrix();
+        const tgt::mat4& t2w = referenceImage->getParent()->getMappingInformation().getTextureToWorldMatrix();
+        registrationMatrix = w2t * registrationMatrix * t2w;
+
+        tgt::mat4 registrationInverse;
+        if (! registrationMatrix.invert(registrationInverse))
+            tgtAssert(false, "Could not invert registration matrix. This should not happen!");
+
+        // render quad to compute similarity measure by shader
+        _differenceShader->setUniform("_registrationInverse", registrationInverse);
+
+        for (int z = 0; z < size.z; ++z) {
+            float texZ = static_cast<float>(z)/static_cast<float>(size.z) + .5f/static_cast<float>(size.z);
+            _differenceShader->setUniform("_zTex", texZ);
+            _fbo->attachTexture(similarityTex, GL_COLOR_ATTACHMENT0, 0, z);
+            QuadRdr.renderQuad();
+        }
+        _differenceShader->deactivate();
+
+        // detach texture and reduce it
+        ImageData* id = new ImageData(3, size, 1);
+        ImageRepresentationGL::create(id, similarityTex);
+        id->setMappingInformation(referenceImage->getParent()->getMappingInformation());
+        dc->addData(p_differenceImageId.getValue(), id);
+        _fbo->deactivate();
+
+        tgt::TextureUnit::setZeroUnit();
+        LGL_ERROR;
+
+        validate(COMPUTE_DIFFERENCE_IMAGE);
+    }
+
 
 }

modules/registration/processors/similaritymeasure.h

@@ -33,6 +33,7 @@
 #include <string>
 
 #include "tgt/buffer.h"
+#include "tgt/matrix.h"
 #include "tgt/vertexarrayobject.h"
 
 #include "core/pipeline/abstractprocessordecorator.h"
@@ -63,7 +64,8 @@ namespace campvis {
     class SimilarityMeasure : public VisualizationProcessor {
     public:
         enum AdditionalInvalidationLevels {
-            PERFORM_OPTIMIZATION = 1U << 6
+            PERFORM_OPTIMIZATION = 1U << 6,
+            COMPUTE_DIFFERENCE_IMAGE = 1U << 7
         };
 
         /**
@@ -101,6 +103,9 @@ namespace campvis {
 
         ButtonProperty p_computeSimilarity;
 
+        DataNameProperty p_differenceImageId;
+        ButtonProperty p_computeDifferenceImage;
+
         GenericOptionProperty<nlopt::algorithm> p_optimizer;
         ButtonProperty p_performOptimization;
 
@@ -119,12 +124,16 @@ namespace campvis {
 
         float computeSimilarity(const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const tgt::vec3& translation, const tgt::vec3& rotation);
 
+        void generateDifferenceImage(DataContainer* dc, const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const tgt::vec3& translation, const tgt::vec3& rotation);
+
+        static tgt::mat4 euleranglesToMat4(const tgt::vec3& eulerAngles);
 
         static double optimizerFunc(const std::vector<double>& x, std::vector<double>& grad, void* my_func_data);
 
         IVec2Property p_viewportSize;
 
-        tgt::Shader* _shader;                           ///< Shader for slice rendering
+        tgt::Shader* _costFunctionShader;                           ///< Shader for slice rendering
+        tgt::Shader* _differenceShader;                 ///< Shader for computing the difference image
         GlReduction* _glr;
 
         static const std::string loggerCat_;