Currently job artifacts in CI/CD pipelines on LRZ GitLab never expire. Starting from Wed 26.1.2022 the default expiration time will be 30 days (GitLab default). Currently existing artifacts in already completed jobs will not be affected by the change. The latest artifacts for all jobs in the latest successful pipelines will be kept. More information: https://gitlab.lrz.de/help/user/admin_area/settings/continuous_integration.html#default-artifacts-expiration

similaritymeasure.cpp 16.7 KB
Newer Older
1
2
3
4
// ================================================================================================
// 
// This file is part of the CAMPVis Software Framework.
// 
5
// If not explicitly stated otherwise: Copyright (C) 2012-2015, all rights reserved,
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
//      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 "similaritymeasure.h"
26
27
28
#include "cgt/logmanager.h"
#include "cgt/shadermanager.h"
#include "cgt/textureunit.h"
29
30
31
32
33
34
35
36
37
38
39
40
41

#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"

#include "core/tools/glreduction.h"
#include "core/tools/quadrenderer.h"

namespace campvis {
42
namespace registration {
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

    static const GenericOption<std::string> metrics[5] = {
        GenericOption<std::string>("SUM", "Sum"),
        GenericOption<std::string>("SAD", "SAD"),
        GenericOption<std::string>("SSD", "SSD"),
        GenericOption<std::string>("NCC", "NCC"),
        GenericOption<std::string>("SNR", "SNR")
    };

    const std::string SimilarityMeasure::loggerCat_ = "CAMPVis.modules.vis.SimilarityMeasure";

    SimilarityMeasure::SimilarityMeasure()
        : VisualizationProcessor(0)
        , p_referenceId("ReferenceId", "Reference Image", "", DataNameProperty::READ)
        , p_movingId("MovingId", "Moving Image", "", DataNameProperty::READ)
58
59
60
        , p_clipX("clipX", "X Axis Clip Coordinates", cgt::ivec2(0), cgt::ivec2(0), cgt::ivec2(0))
        , p_clipY("clipY", "Y Axis Clip Coordinates", cgt::ivec2(0), cgt::ivec2(0), cgt::ivec2(0))
        , p_clipZ("clipZ", "Z Axis Clip Coordinates", cgt::ivec2(0), cgt::ivec2(0), cgt::ivec2(0))
61
        , p_applyMask("ApplyMask", "Apply Mask", true)
62
63
        , p_translation("Translation", "Moving Image Translation", cgt::vec3(0.f), cgt::vec3(-100.f), cgt::vec3(100.f), cgt::vec3(1.f), cgt::vec3(5.f))
        , p_rotation("Rotation", "Moving Image Rotation", cgt::vec3(0.f), cgt::vec3(-cgt::PIf), cgt::vec3(cgt::PIf), cgt::vec3(.01f), cgt::vec3(7.f))
64
65
66
67
        , p_metric("Metric", "Similarity Metric", metrics, 5)
        , p_computeSimilarity("ComputeSimilarity", "Compute Similarity")
        , p_differenceImageId("DifferenceImageId", "Difference Image", "difference", DataNameProperty::WRITE)
        , p_computeDifferenceImage("ComputeDifferenceImage", "Compute Difference Image")
68
        , p_viewportSize("ViewportSize", "Viewport Size", cgt::ivec2(1), cgt::ivec2(1), cgt::ivec2(1000), cgt::ivec2(1))
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
        , _sadssdCostFunctionShader(0)
        , _nccsnrCostFunctionShader(0)
        , _differenceShader(0)
        , _glr(0)
    {
        addProperty(p_referenceId, INVALID_RESULT | INVALID_PROPERTIES);
        addProperty(p_movingId, VALID);

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

        addProperty(p_translation);
        addProperty(p_rotation);
        addProperty(p_metric);
        addProperty(p_computeSimilarity);

        addProperty(p_differenceImageId, VALID);
        addProperty(p_computeDifferenceImage, INVALID_RESULT | COMPUTE_DIFFERENCE_IMAGE);


        _viewportSizeProperty = &p_viewportSize;
    }

    SimilarityMeasure::~SimilarityMeasure() {

    }

    void SimilarityMeasure::init() {
        VisualizationProcessor::init();
        _sadssdCostFunctionShader = ShdrMgr.load("core/glsl/passthrough.vert", "modules/registration/glsl/similaritymeasuresadssd.frag", "");
        _nccsnrCostFunctionShader = ShdrMgr.load("core/glsl/passthrough.vert", "modules/registration/glsl/similaritymeasurenccsnr.frag", "");
        _differenceShader = ShdrMgr.load("core/glsl/passthrough.vert", "modules/registration/glsl/differenceimage.frag", "");

        _glr = new GlReduction(GlReduction::PLUS);
    }

    void SimilarityMeasure::deinit() {
        ShdrMgr.dispose(_sadssdCostFunctionShader);
        ShdrMgr.dispose(_nccsnrCostFunctionShader);
        ShdrMgr.dispose(_differenceShader);

        delete _glr;
        _glr = 0;

        VisualizationProcessor::deinit();
    }

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

        if (referenceImage != 0 && movingImage != 0) {
            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.");
        }
    }

    void SimilarityMeasure::updateProperties(DataContainer& dc) {
        ScopedTypedData<ImageData> referenceImage(dc, p_referenceId.getValue());
        if (referenceImage != 0) {
            p_viewportSize.setValue(referenceImage->getSize().xy());

139
140
141
            p_clipX.setMaxValue(cgt::ivec2(static_cast<int>(referenceImage->getSize().x), static_cast<int>(referenceImage->getSize().x)));
            p_clipY.setMaxValue(cgt::ivec2(static_cast<int>(referenceImage->getSize().y), static_cast<int>(referenceImage->getSize().y)));
            p_clipZ.setMaxValue(cgt::ivec2(static_cast<int>(referenceImage->getSize().z), static_cast<int>(referenceImage->getSize().z)));
142

143
144
145
            p_clipX.setValue(cgt::ivec2(0, static_cast<int>(referenceImage->getSize().x)));
            p_clipY.setValue(cgt::ivec2(0, static_cast<int>(referenceImage->getSize().y)));
            p_clipZ.setValue(cgt::ivec2(0, static_cast<int>(referenceImage->getSize().z)));
146
147
148
        }
    }

149
150
151
    float SimilarityMeasure::computeSimilarity(const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const cgt::vec3& translation, const cgt::vec3& rotation) {
        cgtAssert(referenceImage != 0, "Reference Image must not be 0.");
        cgtAssert(movingImage != 0, "Moving Image must not be 0.");
152

153
        const cgt::svec3& size = referenceImage->getSize();
154
155
156
        p_viewportSize.setValue(size.xy());

        // reserve texture units
157
        cgt::TextureUnit referenceUnit, movingUnit;
158
159
160
        referenceUnit.activate();

        // create temporary texture(s) for result
161
162
163
        cgt::Shader* leShader = _sadssdCostFunctionShader;
        cgt::Texture* similarityTex = 0;
        cgt::Texture* similarityTex2 = 0;
164
        similarityTex = new cgt::Texture(GL_TEXTURE_2D, cgt::ivec3(p_viewportSize.getValue(), 1), GL_RGBA32F, cgt::Texture::NEAREST);
165
        similarityTex->setWrapping(cgt::Texture::CLAMP_TO_EDGE);
166
167
        // NCC and SNR need a second texture and a different shader...
        if (p_metric.getOptionValue() == "NCC" || p_metric.getOptionValue() == "SNR") {
168
            similarityTex2 = new cgt::Texture(GL_TEXTURE_2D, cgt::ivec3(p_viewportSize.getValue(), 1), GL_RGBA32F, cgt::Texture::NEAREST);
169
            similarityTex2->setWrapping(cgt::Texture::CLAMP_TO_EDGE);
170
171
172
173
174
            leShader = _nccsnrCostFunctionShader;
        }

        // activate FBO and attach texture
        _fbo->activate();
175
        const cgt::ivec2& windowSize = p_viewportSize.getValue();
176
177
178
179
180
181
182
183
184
        glViewport(0, 0, static_cast<GLsizei>(windowSize.x), static_cast<GLsizei>(windowSize.y));
        _fbo->attachTexture(similarityTex);
        if (p_metric.getOptionValue() == "NCC" || p_metric.getOptionValue() == "SNR")
            _fbo->attachTexture(similarityTex2, GL_COLOR_ATTACHMENT1);
        LGL_ERROR;

        // bind input images
        leShader->activate();
        leShader->setUniform("_applyMask", p_applyMask.getValue());
185
186
187
        leShader->setUniform("_xClampRange", cgt::vec2(p_clipX.getValue()) / static_cast<float>(size.x));
        leShader->setUniform("_yClampRange", cgt::vec2(p_clipY.getValue()) / static_cast<float>(size.y));
        leShader->setUniform("_zClampRange", cgt::vec2(p_clipZ.getValue()) / static_cast<float>(size.z));
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
        referenceImage->bind(leShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
        movingImage->bind(leShader, movingUnit, "_movingTexture", "_movingTextureParams");

        // render quad to compute similarity measure by shader
        leShader->setUniform("_registrationInverse", computeRegistrationMatrix(referenceImage, movingImage, translation, rotation));
        if (p_metric.getOptionValue() == "NCC" || p_metric.getOptionValue() == "SNR") {
            static const GLenum buffers[] = { GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1};
            glDrawBuffers(2, buffers);
            QuadRdr.renderQuad();
            glDrawBuffers(1, buffers);
        }
        else {
            QuadRdr.renderQuad();
        }
        
        leShader->deactivate();

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

        // reduce the juice
        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 varFixed	= (similarities2[1] - (similarities[2] * similarities[2]) * countRCP) * countRCP;
                float _varMoving = (similarities2[0] - (similarities[1] * similarities[1]) * countRCP) * countRCP;

                if (varFixed > 0.0f && _varMoving > 0.0f) {
                    float correlation = (similarities2[2] - (similarities[1] * similarities[2]) * countRCP) * countRCP;
                    toReturn =  correlation / sqrt(varFixed * _varMoving);
                }
            }
        }
        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);

                toReturn = signal/noise;
            }
        }
        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];
        }

        delete similarityTex;
        delete similarityTex2;
252
        cgt::TextureUnit::setZeroUnit();
253
254
255
256
257
258
        LGL_ERROR;

        return toReturn;
    }


259
    cgt::mat4 SimilarityMeasure::euleranglesToMat4(const cgt::vec3& eulerAngles) {
260
261
262
263
264
265
266
        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);

267
        cgt::mat4 toReturn(cosY * cosZ,   cosZ * sinX * sinY - cosX * sinZ,   sinX * sinZ + cosX * cosZ * sinY,   0.f,
268
269
270
271
272
273
                           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;
    }

274
275
276
277
    void SimilarityMeasure::generateDifferenceImage(DataContainer* dc, const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const cgt::vec3& translation, const cgt::vec3& rotation) {
        cgtAssert(dc != 0, "DataContainer must not be 0.");
        cgtAssert(referenceImage != 0, "Reference Image must not be 0.");
        cgtAssert(movingImage != 0, "Moving Image must not be 0.");
278

279
280
        const cgt::svec3& size = referenceImage->getSize();
        cgt::ivec2 viewportSize = size.xy();
281
282

        // reserve texture units
283
        cgt::TextureUnit referenceUnit, movingUnit;
284
285
286
        referenceUnit.activate();

        // create temporary texture for result
287
        cgt::Texture* differenceImage = new cgt::Texture(GL_TEXTURE_3D, cgt::ivec3(size), GL_R32F, cgt::Texture::LINEAR);
288
        differenceImage->setWrapping(cgt::Texture::CLAMP_TO_EDGE);
289
290
291
292

        // bind input images
        _differenceShader->activate();
        _differenceShader->setUniform("_applyMask", p_applyMask.getValue());
293
294
295
        _differenceShader->setUniform("_xClampRange", cgt::vec2(p_clipX.getValue()) / static_cast<float>(size.x));
        _differenceShader->setUniform("_yClampRange", cgt::vec2(p_clipY.getValue()) / static_cast<float>(size.y));
        _differenceShader->setUniform("_zClampRange", cgt::vec2(p_clipZ.getValue()) / static_cast<float>(size.z));
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
        referenceImage->bind(_differenceShader, referenceUnit, "_referenceTexture", "_referenceTextureParams");
        movingImage->bind(_differenceShader, movingUnit, "_movingTexture", "_movingTextureParams");
        _differenceShader->setUniform("_registrationInverse", computeRegistrationMatrix(referenceImage, movingImage, translation, rotation));

        // activate FBO and attach texture
        _fbo->activate();
        glViewport(0, 0, static_cast<GLsizei>(viewportSize.x), static_cast<GLsizei>(viewportSize.y));

        // render quad to compute difference measure by shader
        for (int z = 0; z < static_cast<int>(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(differenceImage, GL_COLOR_ATTACHMENT0, 0, z);
            QuadRdr.renderQuad();
        }
        _differenceShader->deactivate();
        _fbo->deactivate();

        // put difference image into DataContainer
        ImageData* id = new ImageData(3, size, 1);
        ImageRepresentationGL::create(id, differenceImage);
        id->setMappingInformation(referenceImage->getParent()->getMappingInformation());
        dc->addData(p_differenceImageId.getValue(), id);

320
        cgt::TextureUnit::setZeroUnit();
321
322
323
324
325
        LGL_ERROR;

        validate(COMPUTE_DIFFERENCE_IMAGE);
    }

326
327
328
    cgt::mat4 SimilarityMeasure::computeRegistrationMatrix(const ImageRepresentationGL* referenceImage, const ImageRepresentationGL* movingImage, const cgt::vec3& translation, const cgt::vec3& rotation) {
        cgt::mat4 registrationMatrix = cgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
        cgt::mat4 registrationInverse;
329
        if (! registrationMatrix.invert(registrationInverse))
330
            cgtAssert(false, "Could not invert registration matrix. This should not happen!");
331

332
333
334
335
336
        cgt::Bounds movingBounds = movingImage->getParent()->getWorldBounds();
        cgt::vec3 halfDiagonal = movingBounds.getLLF() + (movingBounds.diagonal() / 2.f);
        const cgt::mat4& w2t = movingImage->getParent()->getMappingInformation().getWorldToTextureMatrix();
        const cgt::mat4& t2w = referenceImage->getParent()->getMappingInformation().getTextureToWorldMatrix();
        return w2t * cgt::mat4::createTranslation(halfDiagonal) * registrationInverse * cgt::mat4::createTranslation(-halfDiagonal) * t2w;
337
338
    }

339
}
340
}