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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
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
139
140
141
142
143
144
145
146
147
148
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
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
|
/*
//
// Copyright 2016 Google, Inc.
//
// Copyright 1997-2009 Torsten Rohlfing
//
// Copyright 2004-2011 SRI International
//
// This file is part of the Computational Morphometry Toolkit.
//
// http://www.nitrc.org/projects/cmtk/
//
// The Computational Morphometry Toolkit is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// The Computational Morphometry Toolkit is distributed in the hope that it
// will be useful, but WITHOUT ANY WARRANTY; without even the implied
// warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with the Computational Morphometry Toolkit. If not, see
// <http://www.gnu.org/licenses/>.
//
// $Revision: 5436 $
//
// $LastChangedDate: 2018-12-10 19:01:20 -0800 (Mon, 10 Dec 2018) $
//
// $LastChangedBy: torstenrohlfing $
//
*/
#ifndef __cmtkVolumeInjectionReconstruction_h_included_
#define __cmtkVolumeInjectionReconstruction_h_included_
#include <cmtkconfig.h>
#include <Registration/cmtkAffineRegistration.h>
#include <Base/cmtkAffineXform.h>
#include <Base/cmtkXform.h>
#include <Base/cmtkUniformVolume.h>
#include "Numerics/ap.h"
#include <vector>
namespace
cmtk
{
/** \addtogroup Recon */
//@{
/** Class for volume reconstruction using volume injection.
*
* This class implements some side effects that are not strictly part of the volume
* injection algorithm, but which supply functionality to the derived
* igsInverseInterpolationVolumeReconstructionBase class. These include computation
* of regional minimum/maximum intensity ranges, between-pass registration, and
* KLD metric computation.
*
*\author Torsten Rohlfing
*/
class VolumeInjectionReconstruction
{
public:
/// This class.
typedef VolumeInjectionReconstruction Self;
/** Constructor from single interleaved image.
* Take original image. Set interleaved image count and stacking axis. Construct separate 3D image stacks for interleaved passes. Allocate corrected image.
*\param originalImage Smart pointer to the original image with motion artifacts.
*\param numberOfPasses The number of interleaved passes, i.e., the number of pass images that comprise the final image.
*\param interleaveAxis Between-slice axis of the interleaved acquisition.
*/
VolumeInjectionReconstruction( const UniformVolume* originalImage, const Types::GridIndexType numberOfPasses, const int interleaveAxis );
/** Constructor for general volume reconstruction from multiple acquired images.
*/
VolumeInjectionReconstruction( const UniformVolume* reconstructionGrid, std::vector<UniformVolume::SmartPtr>& images );
/// Virtual destructor stub.
virtual ~VolumeInjectionReconstruction() {}
/** Static helper function: guess interleaved axis.
* Basically, we assume that images are acquired as interleaved stacks of
* square 2D images with pixel size different from inter-slice spacing.
* So we guess that the interleaved axis is the one that does not match the
* other two axes' image dimensions, and if all three are the same, the
* one that doesn't match their spacing (delta).
*/
static int GuessInterleaveAxis( const UniformVolume* image /*!< The interleaved image.*/,
const int defaultAxis = 2 /*!< In case all guessing fails, this is the default axis we return.*/ );
/** Compute transformations between the reference image grid and the original pass images.
* The resulting transformations are stored in the m_TransformationsToPassImages vector.
* If a high-resolution reference image is set in m_ReferenceImage, then all subimages are registered to it.
* Otherwise, the 0th subimage is used as the reference and the remaining subimages are registered to it.
* In this case, the transformation for the 0th subimage is set as the identity transformation.
*\param registrationMetric Similarity metric for registration of the passes to the reference image.
*/
void ComputeTransformationsToPassImages( const int registrationMetric = 0 );
/// Set transformations to pass images externally (e.g., imported from disk).
void SetTransformationsToPassImages( std::vector<Xform::SmartPtr>& transformations )
{
this->m_TransformationsToPassImages = transformations;
}
/// Get transformation to one pass image.
Xform::SmartPtr& GetTransformationToPassImage( const size_t passIdx )
{
if ( passIdx < this->m_TransformationsToPassImages.size() )
return this->m_TransformationsToPassImages[passIdx];
else
return Xform::SmartPtr::Null();
}
/// Get transformation to one pass image.
std::vector<Xform::SmartPtr>& GetTransformationsToPassImages()
{
return this->m_TransformationsToPassImages;
}
/// Create initial approximation using isotropic volume injection.
void VolumeInjectionIsotropic( const Types::Coordinate kernelSigma /*!< Gaussian kernel sigma (standard deviation) parameter*/,
const Types::Coordinate kernelRadius /*!< Gaussian kernel cut-off radius.*/ );
/// Create initial approximation using anisotropic volume injection.
void VolumeInjectionAnisotropic( const Types::Coordinate kernelSigmaFactor /*!< Gaussian kernel sigma (standard deviation) factor (multiple of per-dimension pass image spacing)*/,
const Types::Coordinate kernelRadiusFactor /*!< Gaussian kernel cut-off radius factor (multiple of per-dimension pass image spacing)*/ );
/// Returns the corrected image.
UniformVolume::SmartPtr& GetCorrectedImage();
/// Set optional separate reference image for motion parameter estimation.
void SetReferenceImage( UniformVolume::SmartPtr& referenceImage );
/** Set pass weight.
* Each pass weight should be between 0 and 1. If the weight for a pass is zero, then that
* pass is effectively excluded from the reconstruction. This can be useful if one of the
* passes shows severe within-pass motion artifacts that would otherwise disturb the
* across-pass correction.
*
* By default, all pass weights are set to 1, i.e., all passes contribute equally.
*/
void SetPassWeight( const size_t pass, const Types::Coordinate weight )
{
this->m_PassWeights[pass] = weight;
}
/// Get Kullback-Leibler Divergence between intensity distributions in original and corrected image.
ap::real_value_type GetOriginalToCorrectedImageKLD( const ap::real_1d_array& x );
protected:
/// Number of interleaved passes.
Types::GridIndexType m_NumberOfPasses;
/// Relative weights of the passes in the correction; can be used to underweight or even exclude passes.
std::vector<Types::Coordinate> m_PassWeights;
/// Original volume pixel intensity range.
Types::DataItemRange m_OriginalImageRange;
/// Original pass images.
std::vector<UniformVolume::SmartPtr> m_OriginalPassImages;
/// Histogram type.
typedef Histogram<double> HistogramType;
/// Original image histogram.
HistogramType::SmartPtr m_OriginalImageHistogram;
/// Corrected image histogram.
HistogramType::SmartPtr m_CorrectedImageHistogram;
/// Original image intensity noise kernel.
std::vector<HistogramType::BinType> m_OriginalImageIntensityNoiseKernel;
/// Optional high-resolution non-interleaved reference image.
UniformVolume::SmartPtr m_ReferenceImage;
/// Affine transformations that map FROM the corrected image TO each of the subimages.
std::vector<Xform::SmartPtr> m_TransformationsToPassImages;
/// Developing corrected image.
UniformVolume::SmartPtr m_CorrectedImage;
/// Corrected image Laplacian.
std::vector<ap::real_value_type> m_CorrectedImageLaplacians;
/** Compute norm of the corrected image Laplacian.
* Side effect: this function first computes the Laplacian image, which is stored in
* m_CorrectedImageLaplacian for use in the AddLaplacianGradientImage function.
*/
ap::real_value_type ComputeCorrectedImageLaplacianNorm( const ap::real_1d_array& correctedImagePixels /*!< Current vector of corrected image pixels.*/ );
/// Add weighted gradient image of Laplacian to already computed cost function gradient.
void AddLaplacianGradientImage(
ap::real_1d_array& g,
const ap::real_1d_array& correctedImagePixels,
const ap::real_value_type weight ) const;
/// Maximum neighborhood pixel values in the corrected image.
ap::real_1d_array m_NeighorhoodMaxPixelValues;
/// Minimum neighborhood pixel values in the corrected image.
ap::real_1d_array m_NeighorhoodMinPixelValues;
private:
/// Number of histogram bins for image entropy estimation.
static const unsigned int NumberOfHistogramBins = 64;
/// Setup kernels and histograms for image entropy estimation.
void SetupHistogramKernels( const TypedArray* originalData );
};
//@}
} // namespace cmtk
#endif // #ifndef __cmtkVolumeInjectionReconstruction_h_included_
|
