nloptregistration.cpp 10.4 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
//      Christian Schulte zu Berge <christian.szb@in.tum.de>
//      Chair for Computer Aided Medical Procedures
8
9
//      Technische Universitaet Muenchen
//      Boltzmannstr. 3, 85748 Garching b. Muenchen, Germany
10
// 
11
12
// For a full list of authors and contributors, please refer to the file "AUTHORS.txt".
// 
13
14
15
16
// 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
17
// 
18
19
20
21
// 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.
22
23
24
25
26
// 
// ================================================================================================

#include "nloptregistration.h"

27
#include "cgt/event/keyevent.h"
28
#include "cgt/opengljobprocessor.h"
29
#include "cgt/painter.h"
30
31
32
33

#include "core/classification/geometry1dtransferfunction.h"
#include "core/classification/tfgeometry1d.h"
#include "core/datastructures/renderdata.h"
34
#include "core/datastructures/transformdata.h"
35
36
37
#include "core/tools/glreduction.h"

namespace campvis {
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
namespace registration {
    namespace {
        cgt::mat4 euleranglesToMat4(const cgt::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);

            cgt::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;
        }
    }

56
57
58
59
60
61
    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)
    };

62
63
    NloptRegistration::NloptRegistration(DataContainer& dc)
        : AutoEvaluationPipeline(dc, getId())
64
65
        , p_optimizer("Optimizer", "Optimizer", optimizers, 3)
        , p_liveUpdate("LiveUpdate", "Live Update of Difference Image", false)
66
67
        , p_performOptimization("PerformOptimization", "Perform Optimization")
        , p_forceStop("Force Stop", "Force Stop")
68
        , p_translationStepSize("TranslationStepSize", "Initial Step Size Translation", 8.f, .1f, 100.f)
69
        , p_rotationStepSize("RotationStepSize", "Initial Step Size Rotation", .5f, .01f, cgt::PIf)
70
        , _lsp()
71
72
        , _referenceReader()
        , _movingReader()
73
        , _rsw(&_canvasSize)
74
75
76
77
        , _sm()
        , _ve(&_canvasSize)
        , _opt(0)
    {
78
        addProcessor(&_lsp);
79
80
        addProcessor(&_referenceReader);
        addProcessor(&_movingReader);
81
        addProcessor(&_rsw);
82
83
84
        addProcessor(&_sm);
        addProcessor(&_ve);

85
86
87
88
89
90
        addProperty(p_optimizer);
        addProperty(p_liveUpdate);
        addProperty(p_performOptimization);
        addProperty(p_forceStop);
        addProperty(p_translationStepSize);
        addProperty(p_rotationStepSize);
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

        p_performOptimization.s_clicked.connect(this, &NloptRegistration::onPerformOptimizationClicked);
        p_forceStop.s_clicked.connect(this, &NloptRegistration::forceStop);

        addEventListenerToBack(&_ve);
    }

    NloptRegistration::~NloptRegistration() {
    }

    void NloptRegistration::init() {
        AutoEvaluationPipeline::init();

        _referenceReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S1median/Volume_2.mhd");
        _referenceReader.p_targetImageID.setValue("Reference Image");
        _referenceReader.p_targetImageID.addSharedProperty(&_sm.p_referenceId);
107
        _referenceReader.p_targetImageID.addSharedProperty(&_rsw.p_sourceImageID);
108
109
110
111

        _movingReader.p_url.setValue("D:/Medical Data/SCR/Data/RegSweeps_phantom_cropped/-1S1median/Volume_3.mhd");
        _movingReader.p_targetImageID.setValue("Moving Image");
        _movingReader.p_targetImageID.addSharedProperty(&_sm.p_movingId);
112
113
114
115
        _movingReader.p_targetImageID.addSharedProperty(&_rsw.p_movingImage);

        _rsw.p_movingTransformationMatrix.setValue("trafoMatrix");
        _rsw.p_targetImageID.setValue("RegistrationSliceView");
116
117
118
119
120

        _sm.p_differenceImageId.addSharedProperty(&_ve.p_inputVolume);
        _sm.p_metric.selectById("NCC");

        _ve.p_outputImage.setValue("volumeexplorer");
121
        _renderTargetID.setValue("RegistrationSliceView");
122

123
124
125
        Geometry1DTransferFunction* dvrTF = new Geometry1DTransferFunction(128, cgt::vec2(-1.f, 1.f));
        dvrTF->addGeometry(TFGeometry1D::createQuad(cgt::vec2(0.f, .5f), cgt::col4(0, 0, 255, 255), cgt::col4(255, 255, 255, 0)));
        dvrTF->addGeometry(TFGeometry1D::createQuad(cgt::vec2(.5f, 1.f), cgt::col4(255, 255, 255, 0), cgt::col4(255, 0, 0, 255)));
126
127
128
        MetaProperty* mp = dynamic_cast<MetaProperty*>(_ve.getProperty("SliceExtractorProperties"));
        dynamic_cast<TransferFunctionProperty*>(mp->getProperty("TransferFunction"))->replaceTF(dvrTF);
        dynamic_cast<TransferFunctionProperty*>(mp->getProperty("TransferFunction"))->setAutoFitWindowToData(false);
129
130
131
132
133
134
135
136
137
138
    }

    void NloptRegistration::deinit() {
        delete _opt;
        _opt = 0;

        AutoEvaluationPipeline::deinit();
    }

    void NloptRegistration::onPerformOptimizationClicked() {
139
        // we want the registration to be performed in a background thread and not in the signal_manager's thread.
140
        std::thread([this] () {
141
142
143
144
145
146
147
148
            // Evaluation of the similarity measure needs an OpenGL context, so we need to acquire the canvas' context.
            // An alternative solution would be to specialize the pipeline and overload the executePipeline() method.
            // Then the registration would be performed in the pipeline's thread, which is probably the mor beautiful
            // solution. 
            cgt::GLContextScopedLock lockGuard(this->_canvas);
            this->performOptimization();
        });

149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
    }

    void NloptRegistration::performOptimization() {
        ImageRepresentationGL::ScopedRepresentation referenceImage(getDataContainer(), _sm.p_referenceId.getValue());
        ImageRepresentationGL::ScopedRepresentation movingImage(getDataContainer(), _sm.p_movingId.getValue());

        if (_opt != 0) {
            LWARNING("Optimization is already running...");
            return;
        }

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

        _opt = new nlopt::opt(p_optimizer.getOptionValue(), 6);
        if (_sm.p_metric.getOptionValue() == "NCC" || _sm.p_metric.getOptionValue() == "SNR") {
            _opt->set_max_objective(&NloptRegistration::optimizerFunc, &mfd);
        }
        else {
            _opt->set_min_objective(&NloptRegistration::optimizerFunc, &mfd);
        }
        _opt->set_xtol_rel(1e-4);

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

        std::vector<double> stepSize(6);
        stepSize[0] = p_translationStepSize.getValue();
        stepSize[1] = p_translationStepSize.getValue();
        stepSize[2] = p_translationStepSize.getValue();
        stepSize[3] = p_rotationStepSize.getValue();
        stepSize[4] = p_rotationStepSize.getValue();
        stepSize[5] = p_rotationStepSize.getValue();
        _opt->set_initial_step(stepSize);

        nlopt::result result = nlopt::SUCCESS;
        try {
190
            double minf;
191
192
193
194
195
196
197
198
            result = _opt->optimize(x, minf);
        }
        catch (std::exception& e) {
            LERROR("Excpetion during optimization: " << e.what());
        }

        if (result >= nlopt::SUCCESS || result <= nlopt::ROUNDOFF_LIMITED) {
            LDEBUG("Optimization successful, took " << mfd._count << " steps.");
199
200
            cgt::vec3 t(x[0], x[1], x[2]);
            cgt::vec3 r(x[3], x[4], x[5]);
201
202
203
204
            _sm.p_translation.setValue(t);
            _sm.p_rotation.setValue(r);

            // compute difference image and render difference volume
205
            _sm.generateDifferenceImage(_dataContainer, referenceImage, movingImage, t, r);
206
207
208
209
210
211
212
213
            _ve.process(getDataContainer());
        }

        delete _opt;
        _opt = 0;
    }

    double NloptRegistration::optimizerFunc(const std::vector<double>& x, std::vector<double>& grad, void* my_func_data) {
214
215
        cgtAssert(x.size() == 6, "Must have 6 values in x.");
        cgtAssert(my_func_data != 0, "my_func_data must not be 0");
216
217
218

        MyFuncData_t* mfd = static_cast<MyFuncData_t*>(my_func_data);
        ++mfd->_count;
219
220
        cgt::vec3 translation(x[0], x[1], x[2]);
        cgt::vec3 rotation(x[3], x[4], x[5]);
221
222
223
224
225
        float similarity = mfd->_object->_sm.computeSimilarity(mfd->_reference, mfd->_moving, translation, rotation);
        LDEBUG(translation << rotation << " : " << similarity);

        // perform interactive update if wished
        if (mfd->_object->p_liveUpdate.getValue()) {
226
227
            cgt::mat4 trafoMatrix = cgt::mat4::createTranslation(translation) * euleranglesToMat4(rotation);
            mfd->_object->getDataContainer().addData("trafoMatrix", new TransformData(trafoMatrix));
228

229
230
            // render slice view
            mfd->_object->_rsw.process(mfd->_object->getDataContainer());
231
232
233
234
235
236
237
238
239
240
241
242

            // update canvas
            mfd->_object->_canvas->getPainter()->paint();
        }

        return similarity;
    }

    void NloptRegistration::forceStop() {
        if (_opt != 0)
            _opt->force_stop();
    }
243
244
    
}
245
246

}