Starting from 2021-07-01, all LRZ GitLab users will be required to explicitly accept the GitLab Terms of Service. Please see the detailed information at https://doku.lrz.de/display/PUBLIC/GitLab and make sure that your projects conform to the requirements.

similaritymeasure.cpp 21.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
// ================================================================================================
// 
// 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.
// 
// ================================================================================================

#include "similaritymeasure.h"
#include "tgt/logmanager.h"
#include "tgt/shadermanager.h"
#include "tgt/textureunit.h"
34
#include "tgt/openglgarbagecollector.h"
35 36 37 38 39 40 41 42 43

#include "core/datastructures/facegeometry.h"
#include "core/datastructures/imagedata.h"
#include "core/datastructures/imagerepresentationgl.h"
#include "core/datastructures/renderdata.h"
#include "core/pipeline/processordecoratorbackground.h"

#include "core/classification/simpletransferfunction.h"

44
#include "core/tools/glreduction.h"
45 46 47
#include "core/tools/quadrenderer.h"

namespace campvis {
48 49 50 51 52
    static const GenericOption<nlopt::algorithm> optimizers[3] = {
        GenericOption<nlopt::algorithm>("cobyla", "COBYLA", nlopt::LN_COBYLA),
        GenericOption<nlopt::algorithm>("newuoa", "NEWUOA", nlopt::LN_NEWUOA),
        GenericOption<nlopt::algorithm>("neldermead", "Nelder-Mead Simplex", nlopt::LN_NELDERMEAD)
    };
53

54
    static const GenericOption<std::string> metrics[5] = {
55
        GenericOption<std::string>("SUM", "Sum"),
56
        GenericOption<std::string>("SAD", "SAD"),
57
        GenericOption<std::string>("SSD", "SSD"),
58 59
        GenericOption<std::string>("NCC", "NCC"),
        GenericOption<std::string>("SNR", "SNR")
60 61
    };

62 63 64 65
    const std::string SimilarityMeasure::loggerCat_ = "CAMPVis.modules.vis.SimilarityMeasure";

    SimilarityMeasure::SimilarityMeasure()
        : VisualizationProcessor(0)
66
        , p_referenceId("ReferenceId", "Reference Image", "", DataNameProperty::READ, AbstractProcessor::INVALID_PROPERTIES)
67
        , p_movingId("MovingId", "Moving Image", "", DataNameProperty::READ, AbstractProcessor::VALID)
68 69 70 71
        , p_clipX("clipX", "X Axis Clip Coordinates", tgt::ivec2(0), tgt::ivec2(0), tgt::ivec2(0))
        , p_clipY("clipY", "Y Axis Clip Coordinates", tgt::ivec2(0), tgt::ivec2(0), tgt::ivec2(0))
        , p_clipZ("clipZ", "Z Axis Clip Coordinates", tgt::ivec2(0), tgt::ivec2(0), tgt::ivec2(0))
        , p_applyMask("ApplyMask", "Apply Mask", true)
72 73
        , p_translation("Translation", "Moving Image Translation", tgt::vec3(0.f), tgt::vec3(-100.f), tgt::vec3(100.f), tgt::vec3(1.f), tgt::vec3(5.f))
        , p_rotation("Rotation", "Moving Image Rotation", tgt::vec3(0.f), tgt::vec3(-tgt::PIf), tgt::vec3(tgt::PIf), tgt::vec3(.01f), tgt::vec3(7.f))
74
        , p_viewportSize("ViewportSize", "Viewport Size", tgt::ivec2(1), tgt::ivec2(1), tgt::ivec2(1000), tgt::ivec2(1), AbstractProcessor::VALID)
75
        , p_metric("Metric", "Similarity Metric", metrics, 5)
76 77 78
        , p_computeSimilarity("ComputeSimilarity", "Compute Similarity")
        , p_optimizer("Optimizer", "Optimizer", optimizers, 3)
        , p_performOptimization("PerformOptimization", "Perform Optimization", AbstractProcessor::INVALID_RESULT | PERFORM_OPTIMIZATION)
79 80
        , p_differenceImageId("DifferenceImageId", "Difference Image", "difference", DataNameProperty::WRITE, AbstractProcessor::VALID)
        , p_computeDifferenceImage("ComputeDifferenceImage", "Compute Difference Image", AbstractProcessor::INVALID_RESULT | COMPUTE_DIFFERENCE_IMAGE)
81
        , p_forceStop("Force Stop", "Force Stop", AbstractProcessor::VALID)
82 83 84
        , _sadssdCostFunctionShader(0)
        , _nccCostFunctionShader(0)
        , _differenceShader(0)
85
        , _glr(0)
86
        , _opt(0)
87 88 89
    {
        addProperty(&p_referenceId);
        addProperty(&p_movingId);
90 91 92 93 94 95

        addProperty(&p_clipX);
        addProperty(&p_clipY);
        addProperty(&p_clipZ);
        addProperty(&p_applyMask);

96 97
        addProperty(&p_translation);
        addProperty(&p_rotation);
98
        addProperty(&p_metric);
99
        addProperty(&p_computeSimilarity);
100

101 102
        addProperty(&p_differenceImageId);
        addProperty(&p_computeDifferenceImage);
103

104 105
        addProperty(&p_optimizer);
        addProperty(&p_performOptimization);
106 107 108
        addProperty(&p_forceStop);

        p_forceStop.s_clicked.connect(this, &SimilarityMeasure::forceStop);
109 110 111 112 113 114 115 116 117 118

        _viewportSizeProperty = &p_viewportSize;
    }

    SimilarityMeasure::~SimilarityMeasure() {

    }

    void SimilarityMeasure::init() {
        VisualizationProcessor::init();
119 120 121 122 123 124 125
        _sadssdCostFunctionShader = ShdrMgr.loadSeparate("core/glsl/passthrough.vert", "modules/registration/glsl/similaritymeasuresadssd.frag", "", false);
        _sadssdCostFunctionShader->setAttributeLocation(0, "in_Position");
        _sadssdCostFunctionShader->setAttributeLocation(1, "in_TexCoord");

        _nccCostFunctionShader = ShdrMgr.loadSeparate("core/glsl/passthrough.vert", "modules/registration/glsl/similaritymeasurencc.frag", "", false);
        _nccCostFunctionShader->setAttributeLocation(0, "in_Position");
        _nccCostFunctionShader->setAttributeLocation(1, "in_TexCoord");
126 127 128 129

        _differenceShader = ShdrMgr.loadSeparate("core/glsl/passthrough.vert", "modules/registration/glsl/differenceimage.frag", "", false);
        _differenceShader->setAttributeLocation(0, "in_Position");
        _differenceShader->setAttributeLocation(1, "in_TexCoord");
130

131
        _glr = new GlReduction(GlReduction::PLUS);
132 133 134 135
    }

    void SimilarityMeasure::deinit() {
        VisualizationProcessor::deinit();
136 137 138
        ShdrMgr.dispose(_sadssdCostFunctionShader);
        ShdrMgr.dispose(_nccCostFunctionShader);
        ShdrMgr.dispose(_differenceShader);
139 140 141

        delete _glr;
        _glr = 0;
142 143 144

        delete _opt;
        _opt = 0;
145 146 147 148 149 150 151
    }

    void SimilarityMeasure::process(DataContainer& data) {
        ImageRepresentationGL::ScopedRepresentation referenceImage(data, p_referenceId.getValue());
        ImageRepresentationGL::ScopedRepresentation movingImage(data, p_movingId.getValue());

        if (referenceImage != 0 && movingImage != 0) {
152 153 154
            if (getInvalidationLevel() & PERFORM_OPTIMIZATION) {
                performOptimization(referenceImage, movingImage);
            }
155

156
            float similarity = computeSimilarity(referenceImage, movingImage, p_translation.getValue(), p_rotation.getValue());
157
            LDEBUG("Similarity Measure: " << similarity);
158 159 160

            if (getInvalidationLevel() & COMPUTE_DIFFERENCE_IMAGE) 
                generateDifferenceImage(&data, referenceImage, movingImage, p_translation.getValue(), p_rotation.getValue());
161 162 163 164 165 166 167 168 169 170 171 172
        }
        else {
            LERROR("No suitable input image found.");
        }

        validate(INVALID_RESULT);
    }

    void SimilarityMeasure::updateProperties(DataContainer& dc) {
        ScopedTypedData<ImageData> referenceImage(dc, p_referenceId.getValue());
        if (referenceImage != 0) {
            p_viewportSize.setValue(referenceImage->getSize().xy());
173 174 175 176 177 178 179 180

            p_clipX.setMaxValue(tgt::ivec2(static_cast<int>(referenceImage->getSize().x), static_cast<int>(referenceImage->getSize().x)));
            p_clipY.setMaxValue(tgt::ivec2(static_cast<int>(referenceImage->getSize().y), static_cast<int>(referenceImage->getSize().y)));
            p_clipZ.setMaxValue(tgt::ivec2(static_cast<int>(referenceImage->getSize().z), static_cast<int>(referenceImage->getSize().z)));

            p_clipX.setValue(tgt::ivec2(0, static_cast<int>(referenceImage->getSize().x)));
            p_clipY.setValue(tgt::ivec2(0, static_cast<int>(referenceImage->getSize().y)));
            p_clipZ.setValue(tgt::ivec2(0, static_cast<int>(referenceImage->getSize().z)));
181 182 183 184
        }

        validate(AbstractProcessor::INVALID_PROPERTIES);
    }
185 186 187 188 189 190 191

    void SimilarityMeasure::performOptimization(const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage) {
        tgtAssert(referenceImage != 0, "Reference Image must not be 0.");
        tgtAssert(movingImage != 0, "Moving Image must not be 0.");

        MyFuncData_t mfd = { this, referenceImage, movingImage, 0 };

192
        _opt = new nlopt::opt(p_optimizer.getOptionValue(), 6);
193
        if (p_metric.getOptionValue() == "NCC" || p_metric.getOptionValue() == "SNR") {
194 195 196 197 198
            _opt->set_max_objective(&SimilarityMeasure::optimizerFunc, &mfd);
        }
        else {
            _opt->set_min_objective(&SimilarityMeasure::optimizerFunc, &mfd);
        }
199
        _opt->set_xtol_rel(1e-4);
200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215

        std::vector<double> x(6);
        x[0] = p_translation.getValue().x;
        x[1] = p_translation.getValue().y;
        x[2] = p_translation.getValue().z;
        x[3] = p_rotation.getValue().x;
        x[4] = p_rotation.getValue().y;
        x[5] = p_rotation.getValue().z;

        std::vector<double> stepSize(6);
        stepSize[0] = 8.0;
        stepSize[1] = 8.0;
        stepSize[2] = 8.0;
        stepSize[3] = 0.5;
        stepSize[4] = 0.5;
        stepSize[5] = 0.5;
216
        _opt->set_initial_step(stepSize);
217 218 219

        double minf;

220 221 222 223 224 225 226 227 228
        nlopt::result result = nlopt::SUCCESS;
        try {
            result = _opt->optimize(x, minf);
        }
        catch (std::exception& e) {
            LERROR("Excpetion during optimization: " << e.what());
        }
        
        if (result >= nlopt::SUCCESS || result <= nlopt::ROUNDOFF_LIMITED) {
229 230 231 232 233
            LDEBUG("Optimization successful, took " << mfd._count << " steps.");
            p_translation.setValue(tgt::vec3(x[0], x[1], x[2]));
            p_rotation.setValue(tgt::vec3(x[3], x[4], x[5]));
        }

234 235 236
        delete _opt;
        _opt = 0;

237 238 239 240 241 242 243 244 245 246 247 248 249 250
        validate(PERFORM_OPTIMIZATION);
    }

    float SimilarityMeasure::computeSimilarity(const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const tgt::vec3& translation, const tgt::vec3& rotation) {
        tgtAssert(referenceImage != 0, "Reference Image must not be 0.");
        tgtAssert(movingImage != 0, "Moving Image must not be 0.");

        const tgt::svec3& size = referenceImage->getSize();
        p_viewportSize.setValue(size.xy());

        // reserve texture units
        tgt::TextureUnit referenceUnit, movingUnit;
        referenceUnit.activate();

251 252 253 254 255
        // create temporary texture(s) for result
        tgt::Shader* leShader = _sadssdCostFunctionShader;
        tgt::Texture* similarityTex = 0;
        tgt::Texture* similarityTex2 = 0;
        similarityTex = new tgt::Texture(0, tgt::ivec3(p_viewportSize.getValue(), 1), GL_RGBA, GL_RGBA32F, GL_FLOAT, tgt::Texture::NEAREST);
256 257
        similarityTex->uploadTexture();
        similarityTex->setWrapping(tgt::Texture::CLAMP);
258
        if (p_metric.getOptionValue() == "NCC" || p_metric.getOptionValue() == "SNR") {
259 260 261 262 263
            similarityTex2 = new tgt::Texture(0, tgt::ivec3(p_viewportSize.getValue(), 1), GL_RGBA, GL_RGBA32F, GL_FLOAT, tgt::Texture::NEAREST);
            similarityTex2->uploadTexture();
            similarityTex2->setWrapping(tgt::Texture::CLAMP);
            leShader = _nccCostFunctionShader;
        }
264 265 266 267 268 269

        // activate FBO and attach texture
        _fbo->activate();
        const tgt::ivec2& windowSize = p_viewportSize.getValue();
        glViewport(0, 0, static_cast<GLsizei>(windowSize.x), static_cast<GLsizei>(windowSize.y));
        _fbo->attachTexture(similarityTex);
270
        if (p_metric.getOptionValue() == "NCC" || p_metric.getOptionValue() == "SNR")
271
            _fbo->attachTexture(similarityTex2, GL_COLOR_ATTACHMENT1);
272 273 274
        LGL_ERROR;

        // bind input images
275 276 277 278 279 280 281
        leShader->activate();
        leShader->setUniform("_applyMask", p_applyMask.getValue());
        leShader->setUniform("_xClampRange", tgt::vec2(p_clipX.getValue()) / static_cast<float>(size.x));
        leShader->setUniform("_yClampRange", tgt::vec2(p_clipY.getValue()) / static_cast<float>(size.y));
        leShader->setUniform("_zClampRange", tgt::vec2(p_clipZ.getValue()) / static_cast<float>(size.z));
        referenceImage->bind(leShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
        movingImage->bind(leShader, movingUnit, "_movingTexture", "_movingTextureParams");
282

283
        tgt::mat4 registrationMatrix = tgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
284 285 286
        tgt::mat4 registrationInverse;
        if (! registrationMatrix.invert(registrationInverse))
            tgtAssert(false, "Could not invert registration matrix. This should not happen!");
287 288 289

        const tgt::mat4& w2t = movingImage->getParent()->getMappingInformation().getWorldToTextureMatrix();
        const tgt::mat4& t2w = referenceImage->getParent()->getMappingInformation().getTextureToWorldMatrix();
290
        registrationInverse = w2t * registrationInverse * t2w;
291 292 293


        // render quad to compute similarity measure by shader
294
        leShader->setUniform("_registrationInverse", registrationInverse);
295
        if (p_metric.getOptionValue() == "NCC" || p_metric.getOptionValue() == "SNR") {
296 297 298 299 300 301 302 303 304 305
            static const GLenum buffers[] = { GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1};
            glDrawBuffers(2, buffers);
            QuadRdr.renderQuad();
            glDrawBuffers(1, buffers);
        }
        else {
            QuadRdr.renderQuad();
        }
        
        leShader->deactivate();
306 307 308 309 310 311 312

        // detach texture and reduce it
        //data.addData("All glory to the HYPNOTOAD!", new RenderData(_fbo));
        _fbo->detachAll();
        _fbo->deactivate();

        // reduce the juice
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
        float toReturn = 0.f;
        if (p_metric.getOptionValue() == "NCC") {
            std::vector<float> similarities = _glr->reduce(similarityTex);
            std::vector<float> similarities2 = _glr->reduce(similarityTex2);

            if (similarities.size() >= 3 && similarities2.size() >= 3) {
                float countRCP = 1.f / similarities[0];
                float meanFixed = similarities[1] * countRCP;
                float meanMoving = similarities[2] * countRCP;
                float varFixed	= (similarities2[1] - (similarities[2] * similarities[2]) * countRCP) * countRCP;
                float _varMoving = (similarities2[0] - (similarities[1] * similarities[1]) * countRCP) * countRCP;

                float correlation = 0.0f;
                if (varFixed > 0.0f && _varMoving > 0.0f)
                {
                    correlation = (similarities2[2] - (similarities[1] * similarities[2]) * countRCP) * countRCP;
                    toReturn =  correlation / sqrt(varFixed * _varMoving);
                }
            }
332 333 334 335 336 337 338 339 340
        }
        else if (p_metric.getOptionValue() == "SNR") {
            std::vector<float> similarities = _glr->reduce(similarityTex);
            std::vector<float> similarities2 = _glr->reduce(similarityTex2);

            if (similarities.size() >= 4 && similarities2.size() >= 4) {
                float countRCP = 1.f / similarities[0];
                float signal = similarities[3] * countRCP;
                float noise = sqrt(similarities2[3] * countRCP);
341

342 343
                toReturn = signal/noise;
            }
344 345 346 347 348 349 350 351 352 353 354
        }
        else {
            std::vector<float> similarities = _glr->reduce(similarityTex);

            if (p_metric.getOptionValue() == "SUM")
                toReturn = similarities[0];
            else if (p_metric.getOptionValue() == "SAD")
                toReturn = similarities[1];
            else if (p_metric.getOptionValue() == "SSD")
                toReturn = similarities[2];
        }
355

356 357
        delete similarityTex;
        delete similarityTex2;
358 359 360
        tgt::TextureUnit::setZeroUnit();
        LGL_ERROR;

361
        return toReturn;
362 363 364 365 366 367 368 369 370 371 372
    }

    double SimilarityMeasure::optimizerFunc(const std::vector<double>& x, std::vector<double>& grad, void* my_func_data) {
        tgtAssert(x.size() == 6, "Must have 6 values in x.");
        tgtAssert(my_func_data != 0, "my_func_data must not be 0");

        MyFuncData_t* mfd = static_cast<MyFuncData_t*>(my_func_data);
        ++mfd->_count;
        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);
373
        GLGC.deleteGarbage();
374

375
        LDEBUG(translation << rotation << " : " << similarity);
376 377 378 379

        return similarity;
    }

380
    tgt::mat4 SimilarityMeasure::euleranglesToMat4(const tgt::vec3& eulerAngles) {
381 382 383 384 385 386
        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);
387 388

        tgt::mat4 toReturn(cosY * cosZ,   cosZ * sinX * sinY - cosX * sinZ,   sinX * sinZ + cosX * cosZ * sinY,   0.f,
389 390
                           cosY * sinZ,   sinX * sinY * sinZ + cosX * cosZ,   cosX * sinY * sinZ - cosZ * sinX,   0.f,
                           (-1) * sinY,   cosY * sinX,                        cosX * cosY,                        0.f,
391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407
                           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
408
        tgt::Texture* similarityTex = new tgt::Texture(0, tgt::ivec3(size), GL_RED, GL_R32F, GL_FLOAT, tgt::Texture::NEAREST);
409 410 411 412 413 414 415 416 417 418
        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();
419 420 421 422
        _differenceShader->setUniform("_applyMask", p_applyMask.getValue());
        _differenceShader->setUniform("_xClampRange", tgt::vec2(p_clipX.getValue()) / static_cast<float>(size.x));
        _differenceShader->setUniform("_yClampRange", tgt::vec2(p_clipY.getValue()) / static_cast<float>(size.y));
        _differenceShader->setUniform("_zClampRange", tgt::vec2(p_clipZ.getValue()) / static_cast<float>(size.z));
423 424 425 426
        referenceImage->bind(_differenceShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
        movingImage->bind(_differenceShader, movingUnit, "_movingTexture", "_movingTextureParams");

        tgt::mat4 registrationMatrix = tgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
427 428 429
        tgt::mat4 registrationInverse;
        if (! registrationMatrix.invert(registrationInverse))
            tgtAssert(false, "Could not invert registration matrix. This should not happen!");
430 431 432

        const tgt::mat4& w2t = movingImage->getParent()->getMappingInformation().getWorldToTextureMatrix();
        const tgt::mat4& t2w = referenceImage->getParent()->getMappingInformation().getTextureToWorldMatrix();
433
        registrationInverse = w2t * registrationInverse * t2w;
434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459


        // 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);
    }

460 461 462 463 464 465 466
    void SimilarityMeasure::forceStop() {
        if (_opt != 0)
            _opt->force_stop();

        validate(PERFORM_OPTIMIZATION);
    }

467

468
}