Mercurial > hg > orthanc-stone

comparison Framework/Toolbox/GeometryToolbox.h @ 144:9b83f30fc1c0 wasm

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
bilinear and trilinear interpolation
author Sebastien Jodogne <s.jodogne@gmail.com>
date Mon, 22 Jan 2018 15:55:03 +0100
parents 2115530d3703
children 62670cc2bb50
comparison
equal deleted inserted replaced
143:58c545177c1c 144:9b83f30fc1c0
26 #include <boost/version.hpp> 26 #include <boost/version.hpp>
27 #if BOOST_VERSION >= 106300 // or 64, need to check 27 #if BOOST_VERSION >= 106300 // or 64, need to check
28 # include <boost/serialization/array_wrapper.hpp> 28 # include <boost/serialization/array_wrapper.hpp>
29 #endif 29 #endif
30 30
31 #include <Core/DicomFormat/DicomMap.h>
32
31 #include <boost/numeric/ublas/matrix.hpp> 33 #include <boost/numeric/ublas/matrix.hpp>
32 #include <boost/numeric/ublas/vector.hpp> 34 #include <boost/numeric/ublas/vector.hpp>
33 35 #include <cassert>
34 #include <Core/DicomFormat/DicomMap.h>
35 36
36 namespace OrthancStone 37 namespace OrthancStone
37 { 38 {
38 typedef boost::numeric::ublas::matrix<double> Matrix; 39 typedef boost::numeric::ublas::matrix<double> Matrix;
39 typedef boost::numeric::ublas::vector<double> Vector; 40 typedef boost::numeric::ublas::vector<double> Vector;
123 124
124 Matrix CreateRotationMatrixAlongX(double a); 125 Matrix CreateRotationMatrixAlongX(double a);
125 126
126 Matrix CreateRotationMatrixAlongY(double a); 127 Matrix CreateRotationMatrixAlongY(double a);
127 128
128 Matrix CreateRotationMatrixAlongZ(double a); 129 Matrix CreateRotationMatrixAlongZ(double a);
130
131 inline float ComputeBilinearInterpolationInternal(float x,
132 float y,
133 float f00, // source(x, y)
134 float f01, // source(x + 1, y)
135 float f10, // source(x, y + 1)
136 float f11); // source(x + 1, y + 1)
137
138 inline float ComputeBilinearInterpolation(float x,
139 float y,
140 float f00, // source(x, y)
141 float f01, // source(x + 1, y)
142 float f10, // source(x, y + 1)
143 float f11); // source(x + 1, y + 1)
144
145 inline float ComputeTrilinearInterpolation(float x,
146 float y,
147 float z,
148 float f000, // source(x, y, z)
149 float f001, // source(x + 1, y, z)
150 float f010, // source(x, y + 1, z)
151 float f011, // source(x + 1, y + 1, z)
152 float f100, // source(x, y, z + 1)
153 float f101, // source(x + 1, y, z + 1)
154 float f110, // source(x, y + 1, z + 1)
155 float f111); // source(x + 1, y + 1, z + 1)
129 }; 156 };
130 } 157 }
158
159
160 float OrthancStone::GeometryToolbox::ComputeBilinearInterpolationInternal(float x,
161 float y,
162 float f00,
163 float f01,
164 float f10,
165 float f11)
166 {
167 // This function works on fractional parts
168 assert(x >= 0 && y >= 0 && x < 1 && y < 1);
169
170 // https://en.wikipedia.org/wiki/Bilinear_interpolation#Unit_square
171 return f00 * (1 - x) * (1 - y) + f01 * x * (1 - y) + f10 * (1 - x) * y + f11 * x * y;
172 }
173
174
175 float OrthancStone::GeometryToolbox::ComputeBilinearInterpolation(float x,
176 float y,
177 float f00,
178 float f01,
179 float f10,
180 float f11)
181 {
182 assert(x >= 0 && y >= 0);
183
184 // Compute the fractional part of (x,y) (this is the "floor()"
185 // operation on positive integers)
186 float xx = x - static_cast<float>(static_cast<int>(x));
187 float yy = y - static_cast<float>(static_cast<int>(y));
188
189 return ComputeBilinearInterpolationInternal(xx, yy, f00, f01, f10, f11);
190 }
191
192
193 float OrthancStone::GeometryToolbox::ComputeTrilinearInterpolation(float x,
194 float y,
195 float z,
196 float f000,
197 float f001,
198 float f010,
199 float f011,
200 float f100,
201 float f101,
202 float f110,
203 float f111)
204 {
205 assert(x >= 0 && y >= 0 && z >= 0);
206
207 float xx = x - static_cast<float>(static_cast<int>(x));
208 float yy = y - static_cast<float>(static_cast<int>(y));
209 float zz = z - static_cast<float>(static_cast<int>(z));
210
211 // "In practice, a trilinear interpolation is identical to two
212 // bilinear interpolation combined with a linear interpolation"
213 // https://en.wikipedia.org/wiki/Trilinear_interpolation#Method
214 float a = ComputeBilinearInterpolationInternal(xx, yy, f000, f001, f010, f011);
215 float b = ComputeBilinearInterpolationInternal(xx, yy, f100, f101, f110, f111);
216
217 return (1 - zz) * a + zz * b;
218 }