Field3D
|
00001 //----------------------------------------------------------------------------// 00002 00003 /* 00004 * Copyright (c) 2009 Sony Pictures Imageworks Inc 00005 * 00006 * All rights reserved. 00007 * 00008 * Redistribution and use in source and binary forms, with or without 00009 * modification, are permitted provided that the following conditions 00010 * are met: 00011 * 00012 * Redistributions of source code must retain the above copyright 00013 * notice, this list of conditions and the following disclaimer. 00014 * Redistributions in binary form must reproduce the above copyright 00015 * notice, this list of conditions and the following disclaimer in the 00016 * documentation and/or other materials provided with the 00017 * distribution. Neither the name of Sony Pictures Imageworks nor the 00018 * names of its contributors may be used to endorse or promote 00019 * products derived from this software without specific prior written 00020 * permission. 00021 * 00022 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00023 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 00024 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 00025 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 00026 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 00027 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 00028 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 00029 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 00030 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, 00031 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 00032 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED 00033 * OF THE POSSIBILITY OF SUCH DAMAGE. 00034 */ 00035 00036 //----------------------------------------------------------------------------// 00037 00044 //----------------------------------------------------------------------------// 00045 00046 #ifndef _INCLUDED_Field3D_FieldInterp_H_ 00047 #define _INCLUDED_Field3D_FieldInterp_H_ 00048 00049 #include "Field.h" 00050 #include "DenseField.h" 00051 #include "MACField.h" 00052 #include "ProceduralField.h" 00053 #include "RefCount.h" 00054 //----------------------------------------------------------------------------// 00055 00056 #include "ns.h" 00057 00058 FIELD3D_NAMESPACE_OPEN 00059 00060 //----------------------------------------------------------------------------// 00061 // FieldInterp 00062 //----------------------------------------------------------------------------// 00063 00070 template <class Data_T> 00071 class FieldInterp : public RefBase 00072 { 00073 public: 00074 typedef boost::intrusive_ptr<FieldInterp> Ptr; 00075 virtual ~FieldInterp() 00076 { } 00077 virtual Data_T sample(const Field<Data_T> &data, const V3d &vsP) const = 0; 00078 }; 00079 00080 //----------------------------------------------------------------------------// 00081 // LinearFieldInterp 00082 //----------------------------------------------------------------------------// 00083 00084 /* \class LinearFieldInterp 00085 \ingroup field 00086 \brief Basic linear interpolator using voxel access through Field base class 00087 */ 00088 00089 //----------------------------------------------------------------------------// 00090 00091 template <class Data_T> 00092 class LinearFieldInterp : public FieldInterp<Data_T> 00093 { 00094 public: 00095 typedef boost::intrusive_ptr<LinearFieldInterp> Ptr; 00096 virtual Data_T sample(const Field<Data_T> &data, const V3d &vsP) const; 00097 }; 00098 00099 //----------------------------------------------------------------------------// 00100 // CubicFieldInterp 00101 //----------------------------------------------------------------------------// 00102 00103 /* \class CubicFieldInterp 00104 \ingroup field 00105 \brief Basic cubic interpolator using voxel access through Field base class 00106 */ 00107 00108 //----------------------------------------------------------------------------// 00109 00110 template <class Data_T> 00111 class CubicFieldInterp : public FieldInterp<Data_T> 00112 { 00113 public: 00114 typedef boost::intrusive_ptr<CubicFieldInterp> Ptr; 00115 virtual Data_T sample(const Field<Data_T> &data, const V3d &vsP) const; 00116 }; 00117 00118 //----------------------------------------------------------------------------// 00119 // LinearGenericFieldInterp 00120 //----------------------------------------------------------------------------// 00121 00122 /* \class LinearGenericFieldInterp 00123 \ingroup field 00124 \brief Linear interpolator optimized for fields with a fastValue function 00125 */ 00126 00127 //----------------------------------------------------------------------------// 00128 00129 template <class Field_T> 00130 class LinearGenericFieldInterp : public RefBase 00131 { 00132 public: 00133 typedef boost::intrusive_ptr<LinearGenericFieldInterp> Ptr; 00134 typename Field_T::value_type sample(const Field_T &data, const V3d &vsP) const; 00135 }; 00136 00137 //----------------------------------------------------------------------------// 00138 // LinearMACFieldInterp 00139 //----------------------------------------------------------------------------// 00140 00141 /* \class LinearMACFieldInterp 00142 \ingroup field 00143 \brief Linear interpolator optimized for the MAC fields 00144 */ 00145 00146 //----------------------------------------------------------------------------// 00147 00148 template <class Data_T> 00149 class LinearMACFieldInterp : public RefBase 00150 { 00151 public: 00152 00153 typedef boost::intrusive_ptr<LinearMACFieldInterp> Ptr; 00154 00155 Data_T sample(const MACField<Data_T> &data, const V3d &vsP) const; 00156 }; 00157 00158 //----------------------------------------------------------------------------// 00159 // CubicGenericFieldInterp 00160 //----------------------------------------------------------------------------// 00161 00162 /* \class CubicGenericFieldInterp 00163 \ingroup field 00164 \brief Cubic interpolator optimized for fields with a fastValue function 00165 */ 00166 00167 //----------------------------------------------------------------------------// 00168 00169 template <class Field_T> 00170 class CubicGenericFieldInterp : public RefBase 00171 { 00172 public: 00173 typedef boost::intrusive_ptr<CubicGenericFieldInterp> Ptr; 00174 typename Field_T::value_type sample(const Field_T &data, const V3d &vsP) const; 00175 }; 00176 00177 //----------------------------------------------------------------------------// 00178 // CubicMACFieldInterp 00179 //----------------------------------------------------------------------------// 00180 00181 /* \class CubicMACFieldInterp 00182 \ingroup field 00183 \brief Linear interpolator optimized for MAC fields 00184 */ 00185 00186 //----------------------------------------------------------------------------// 00187 00188 template <class Data_T> 00189 class CubicMACFieldInterp : public RefBase 00190 { 00191 public: 00192 typedef boost::intrusive_ptr<CubicMACFieldInterp> Ptr; 00193 Data_T sample(const MACField<Data_T> &data, const V3d &vsP) const; 00194 }; 00195 00196 //----------------------------------------------------------------------------// 00197 // ProceduralFieldLookup 00198 //----------------------------------------------------------------------------// 00199 00200 /* \class ProceduralFieldLookup 00201 \ingroup field 00202 \brief "Interpolator" for procedural fields - point samples instead of 00203 interpolating. 00204 */ 00205 00206 //----------------------------------------------------------------------------// 00207 00208 template <class Data_T> 00209 class ProceduralFieldLookup : public RefBase 00210 { 00211 public: 00212 typedef boost::intrusive_ptr<ProceduralFieldLookup> Ptr; 00213 Data_T sample(const ProceduralField<Data_T> &data, 00214 const V3d &vsP) const; 00215 }; 00216 00217 //----------------------------------------------------------------------------// 00218 // Interpolation functions 00219 //----------------------------------------------------------------------------// 00220 00222 template <class Data_T> 00223 Data_T wsSample(const typename Field<Data_T>::Ptr f, 00224 const FieldInterp<Data_T> &interp, 00225 const V3d &wsP) 00226 { 00227 V3d vsP; 00228 f->mapping()->worldToVoxel(wsP, vsP); 00229 return interp.sample(*f, vsP); 00230 } 00231 00232 //----------------------------------------------------------------------------// 00233 // Interpolation helper functions 00234 //----------------------------------------------------------------------------// 00235 00237 bool isPointInField(const FieldRes::Ptr f, const V3d &wsP); 00238 00239 //----------------------------------------------------------------------------// 00240 00243 bool isLegalVoxelCoord(const V3d &vsP, const Box3d &vsDataWindow); 00244 00245 //----------------------------------------------------------------------------// 00246 // Math functions 00247 //----------------------------------------------------------------------------// 00248 00251 template <class S, class T> 00252 FIELD3D_VEC3_T<T> operator * (S s, const FIELD3D_VEC3_T<T> vec); 00253 00254 //----------------------------------------------------------------------------// 00255 // Interpolants 00256 //----------------------------------------------------------------------------// 00257 00262 template <class Data_T> 00263 Data_T monotonicCubicInterpolant(const Data_T &f1, const Data_T &f2, 00264 const Data_T &f3, const Data_T &f4, 00265 double t); 00266 00267 //----------------------------------------------------------------------------// 00268 00273 template <class Data_T> 00274 Data_T monotonicCubicInterpolantVec(const Data_T &f1, const Data_T &f2, 00275 const Data_T &f3, const Data_T &f4, 00276 double t); 00277 00278 //----------------------------------------------------------------------------// 00279 // Implementations 00280 //----------------------------------------------------------------------------// 00281 00282 template <class Data_T> 00283 Data_T LinearFieldInterp<Data_T>::sample(const Field<Data_T> &data, 00284 const V3d &vsP) const 00285 { 00286 // Voxel centers are at .5 coordinates 00287 // NOTE: Don't use contToDisc for this, we're looking for sample 00288 // point locations, not coordinate shifts. 00289 FIELD3D_VEC3_T<double> p(vsP - FIELD3D_VEC3_T<double>(0.5)); 00290 00291 // Lower left corner 00292 V3i c1(static_cast<int>(floor(p.x)), 00293 static_cast<int>(floor(p.y)), 00294 static_cast<int>(floor(p.z))); 00295 // Upper right corner 00296 V3i c2(c1 + V3i(1)); 00297 // C1 fractions 00298 FIELD3D_VEC3_T<double> f1(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00299 // C2 fraction 00300 FIELD3D_VEC3_T<double> f2(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00301 00302 // Clamp the indexing coordinates 00303 if (true) { 00304 const Box3i &dataWindow = data.dataWindow(); 00305 c1.x = std::max(dataWindow.min.x, std::min(c1.x, dataWindow.max.x)); 00306 c2.x = std::max(dataWindow.min.x, std::min(c2.x, dataWindow.max.x)); 00307 c1.y = std::max(dataWindow.min.y, std::min(c1.y, dataWindow.max.y)); 00308 c2.y = std::max(dataWindow.min.y, std::min(c2.y, dataWindow.max.y)); 00309 c1.z = std::max(dataWindow.min.z, std::min(c1.z, dataWindow.max.z)); 00310 c2.z = std::max(dataWindow.min.z, std::min(c2.z, dataWindow.max.z)); 00311 } 00312 00313 return static_cast<Data_T> 00314 (f1.x * (f1.y * (f1.z * data.value(c1.x, c1.y, c1.z) + 00315 f2.z * data.value(c1.x, c1.y, c2.z)) + 00316 f2.y * (f1.z * data.value(c1.x, c2.y, c1.z) + 00317 f2.z * data.value(c1.x, c2.y, c2.z))) + 00318 f2.x * (f1.y * (f1.z * data.value(c2.x, c1.y, c1.z) + 00319 f2.z * data.value(c2.x, c1.y, c2.z)) + 00320 f2.y * (f1.z * data.value(c2.x, c2.y, c1.z) + 00321 f2.z * data.value(c2.x, c2.y, c2.z)))); 00322 00323 } 00324 00325 //----------------------------------------------------------------------------// 00326 00327 template <class Data_T> 00328 Data_T CubicFieldInterp<Data_T>::sample(const Field<Data_T> &data, 00329 const V3d &vsP) const 00330 { 00331 // Voxel centers are at .5 coordinates 00332 // NOTE: Don't use contToDisc for this, we're looking for sample 00333 // point locations, not coordinate shifts. 00334 V3d clampedVsP(std::max(0.5, vsP.x), 00335 std::max(0.5, vsP.y), 00336 std::max(0.5, vsP.z)); 00337 FIELD3D_VEC3_T<double> p(clampedVsP - FIELD3D_VEC3_T<double>(0.5)); 00338 00339 // Lower left corner 00340 V3i c(static_cast<int>(floor(p.x)), 00341 static_cast<int>(floor(p.y)), 00342 static_cast<int>(floor(p.z))); 00343 00344 // Fractions 00345 FIELD3D_VEC3_T<double> t(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 00346 00347 const Box3i &dataWindow = data.dataWindow(); 00348 00349 // Clamp the coordinates 00350 int im, jm, km; 00351 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 00352 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 00353 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 00354 int im_1, jm_1, km_1; 00355 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 00356 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 00357 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 00358 int im1, jm1, km1; 00359 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 00360 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 00361 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 00362 int im2, jm2, km2; 00363 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 00364 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 00365 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 00366 00367 // interpolate 16 lines in z: 00368 Data_T z11 = monotonicCubicInterpolant(data.value(im_1, jm_1, km_1), 00369 data.value(im_1, jm_1, km), 00370 data.value(im_1, jm_1, km1), 00371 data.value(im_1, jm_1, km2), t.z); 00372 Data_T z12 = monotonicCubicInterpolant(data.value(im_1, jm, km_1), 00373 data.value(im_1, jm, km), 00374 data.value(im_1, jm, km1), 00375 data.value(im_1, jm, km2), t.z); 00376 Data_T z13 = monotonicCubicInterpolant(data.value(im_1, jm1, km_1), 00377 data.value(im_1, jm1, km), 00378 data.value(im_1, jm1, km1), 00379 data.value(im_1, jm1, km2), t.z); 00380 Data_T z14 = monotonicCubicInterpolant(data.value(im_1, jm2, km_1), 00381 data.value(im_1, jm2, km), 00382 data.value(im_1, jm2, km1), 00383 data.value(im_1, jm2, km2), t.z); 00384 00385 Data_T z21 = monotonicCubicInterpolant(data.value(im, jm_1, km_1), 00386 data.value(im, jm_1, km), 00387 data.value(im, jm_1, km1), 00388 data.value(im, jm_1, km2), t.z); 00389 Data_T z22 = monotonicCubicInterpolant(data.value(im, jm, km_1), 00390 data.value(im, jm, km), 00391 data.value(im, jm, km1), 00392 data.value(im, jm, km2), t.z); 00393 Data_T z23 = monotonicCubicInterpolant(data.value(im, jm1, km_1), 00394 data.value(im, jm1, km), 00395 data.value(im, jm1, km1), 00396 data.value(im, jm1, km2), t.z); 00397 Data_T z24 = monotonicCubicInterpolant(data.value(im, jm2, km_1), 00398 data.value(im, jm2, km), 00399 data.value(im, jm2, km1), 00400 data.value(im, jm2, km2), t.z); 00401 00402 Data_T z31 = monotonicCubicInterpolant(data.value(im1, jm_1, km_1), 00403 data.value(im1, jm_1, km), 00404 data.value(im1, jm_1, km1), 00405 data.value(im1, jm_1, km2), t.z); 00406 Data_T z32 = monotonicCubicInterpolant(data.value(im1, jm, km_1), 00407 data.value(im1, jm, km), 00408 data.value(im1, jm, km1), 00409 data.value(im1, jm, km2), t.z); 00410 Data_T z33 = monotonicCubicInterpolant(data.value(im1, jm1, km_1), 00411 data.value(im1, jm1, km), 00412 data.value(im1, jm1, km1), 00413 data.value(im1, jm1, km2), t.z); 00414 Data_T z34 = monotonicCubicInterpolant(data.value(im1, jm2, km_1), 00415 data.value(im1, jm2, km), 00416 data.value(im1, jm2, km1), 00417 data.value(im1, jm2, km2), t.z); 00418 00419 Data_T z41 = monotonicCubicInterpolant(data.value(im2, jm_1, km_1), 00420 data.value(im2, jm_1, km), 00421 data.value(im2, jm_1, km1), 00422 data.value(im2, jm_1, km2), t.z); 00423 Data_T z42 = monotonicCubicInterpolant(data.value(im2, jm, km_1), 00424 data.value(im2, jm, km), 00425 data.value(im2, jm, km1), 00426 data.value(im2, jm, km2), t.z); 00427 Data_T z43 = monotonicCubicInterpolant(data.value(im2, jm1, km_1), 00428 data.value(im2, jm1, km), 00429 data.value(im2, jm1, km1), 00430 data.value(im2, jm1, km2), t.z); 00431 Data_T z44 = monotonicCubicInterpolant(data.value(im2, jm2, km_1), 00432 data.value(im2, jm2, km), 00433 data.value(im2, jm2, km1), 00434 data.value(im2, jm2, km2), t.z); 00435 00436 Data_T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 00437 Data_T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 00438 Data_T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 00439 Data_T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 00440 00441 Data_T z0 = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 00442 00443 return z0; 00444 } 00445 00446 //----------------------------------------------------------------------------// 00447 00448 template <class Field_T> 00449 typename Field_T::value_type 00450 LinearGenericFieldInterp<Field_T>::sample(const Field_T &data, 00451 const V3d &vsP) const 00452 { 00453 typedef typename Field_T::value_type Data_T; 00454 00455 // Pixel centers are at .5 coordinates 00456 // NOTE: Don't use contToDisc for this, we're looking for sample 00457 // point locations, not coordinate shifts. 00458 FIELD3D_VEC3_T<double> p(vsP - FIELD3D_VEC3_T<double>(0.5)); 00459 00460 // Lower left corner 00461 V3i c1(static_cast<int>(floor(p.x)), 00462 static_cast<int>(floor(p.y)), 00463 static_cast<int>(floor(p.z))); 00464 // Upper right corner 00465 V3i c2(c1 + V3i(1)); 00466 // C1 fractions 00467 FIELD3D_VEC3_T<double> f1(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00468 // C2 fraction 00469 FIELD3D_VEC3_T<double> f2(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00470 00471 const Box3i &dataWindow = data.dataWindow(); 00472 00473 // Clamp the coordinates 00474 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 00475 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 00476 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 00477 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 00478 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 00479 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 00480 00481 return static_cast<Data_T> 00482 (f1.x * (f1.y * (f1.z * data.fastValue(c1.x, c1.y, c1.z) + 00483 f2.z * data.fastValue(c1.x, c1.y, c2.z)) + 00484 f2.y * (f1.z * data.fastValue(c1.x, c2.y, c1.z) + 00485 f2.z * data.fastValue(c1.x, c2.y, c2.z))) + 00486 f2.x * (f1.y * (f1.z * data.fastValue(c2.x, c1.y, c1.z) + 00487 f2.z * data.fastValue(c2.x, c1.y, c2.z)) + 00488 f2.y * (f1.z * data.fastValue(c2.x, c2.y, c1.z) + 00489 f2.z * data.fastValue(c2.x, c2.y, c2.z)))); 00490 } 00491 00492 //----------------------------------------------------------------------------// 00493 00494 template <class Data_T> 00495 Data_T LinearMACFieldInterp<Data_T>::sample(const MACField<Data_T> &data, 00496 const V3d &vsP) const 00497 { 00498 // Pixel centers are at .5 coordinates 00499 // NOTE: Don't use contToDisc for this, we're looking for sample 00500 // point locations, not coordinate shifts. 00501 00502 const Box3i &dataWindow = data.dataWindow(); 00503 00504 Data_T ret; 00505 00506 FIELD3D_VEC3_T<double> p(vsP.x , vsP.y - 0.5, vsP.z - 0.5); 00507 00508 // X component --- 00509 00510 // Lower left corner 00511 V3i c1(static_cast<int>(floor(p.x)), 00512 static_cast<int>(floor(p.y)), 00513 static_cast<int>(floor(p.z))); 00514 00515 // Upper right corner 00516 V3i c2(c1 + V3i(1)); 00517 00518 // C1 fractions 00519 FIELD3D_VEC3_T<double> f1(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00520 // C2 fraction 00521 FIELD3D_VEC3_T<double> f2(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00522 00523 // Clamp the coordinates 00524 c1.x = std::min(dataWindow.max.x + 1, std::max(dataWindow.min.x, c1.x)); 00525 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 00526 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 00527 c2.x = std::min(dataWindow.max.x + 1, std::max(dataWindow.min.x, c2.x)); 00528 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 00529 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 00530 00531 ret.x = (f1.x * (f1.y * (f1.z * data.u(c1.x, c1.y, c1.z) + 00532 f2.z * data.u(c1.x, c1.y, c2.z)) + 00533 f2.y * (f1.z * data.u(c1.x, c2.y, c1.z) + 00534 f2.z * data.u(c1.x, c2.y, c2.z))) + 00535 f2.x * (f1.y * (f1.z * data.u(c2.x, c1.y, c1.z) + 00536 f2.z * data.u(c2.x, c1.y, c2.z)) + 00537 f2.y * (f1.z * data.u(c2.x, c2.y, c1.z) + 00538 f2.z * data.u(c2.x, c2.y, c2.z)))); 00539 00540 // Y component --- 00541 00542 p.setValue(vsP.x - 0.5, vsP.y , vsP.z - 0.5); 00543 00544 // Lower left corner 00545 c1.x = static_cast<int>(floor(p.x )); 00546 c1.y = static_cast<int>(floor(p.y )); 00547 c1.z = static_cast<int>(floor(p.z )); 00548 00549 // Upper right corner 00550 c2.x = c1.x + 1; 00551 c2.y = c1.y + 1; 00552 c2.z = c1.z + 1; 00553 00554 // C1 fractions 00555 f1.setValue(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00556 // C2 fraction 00557 f2.setValue(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00558 00559 // Clamp the coordinates 00560 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 00561 c1.y = std::min(dataWindow.max.y + 1, std::max(dataWindow.min.y, c1.y)); 00562 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 00563 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 00564 c2.y = std::min(dataWindow.max.y + 1, std::max(dataWindow.min.y, c2.y)); 00565 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 00566 00567 ret.y = (f1.x * (f1.y * (f1.z * data.v(c1.x, c1.y, c1.z) + 00568 f2.z * data.v(c1.x, c1.y, c2.z)) + 00569 f2.y * (f1.z * data.v(c1.x, c2.y, c1.z) + 00570 f2.z * data.v(c1.x, c2.y, c2.z))) + 00571 f2.x * (f1.y * (f1.z * data.v(c2.x, c1.y, c1.z) + 00572 f2.z * data.v(c2.x, c1.y, c2.z)) + 00573 f2.y * (f1.z * data.v(c2.x, c2.y, c1.z) + 00574 f2.z * data.v(c2.x, c2.y, c2.z)))); 00575 00576 // Z component --- 00577 00578 p.setValue(vsP.x - 0.5 , vsP.y - 0.5, vsP.z); 00579 00580 // Lower left corner 00581 c1.x = static_cast<int>(floor(p.x )); 00582 c1.y = static_cast<int>(floor(p.y )); 00583 c1.z = static_cast<int>(floor(p.z )); 00584 00585 // Upper right corner 00586 c2.x = c1.x + 1; 00587 c2.y = c1.y + 1; 00588 c2.z = c1.z + 1; 00589 00590 // C1 fractions 00591 f1.setValue(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00592 // C2 fraction 00593 f2.setValue(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00594 00595 // Clamp the coordinates 00596 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 00597 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 00598 c1.z = std::min(dataWindow.max.z + 1, std::max(dataWindow.min.z, c1.z)); 00599 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 00600 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 00601 c2.z = std::min(dataWindow.max.z + 1, std::max(dataWindow.min.z, c2.z)); 00602 00603 ret.z = (f1.x * (f1.y * (f1.z * data.w(c1.x, c1.y, c1.z) + 00604 f2.z * data.w(c1.x, c1.y, c2.z)) + 00605 f2.y * (f1.z * data.w(c1.x, c2.y, c1.z) + 00606 f2.z * data.w(c1.x, c2.y, c2.z))) + 00607 f2.x * (f1.y * (f1.z * data.w(c2.x, c1.y, c1.z) + 00608 f2.z * data.w(c2.x, c1.y, c2.z)) + 00609 f2.y * (f1.z * data.w(c2.x, c2.y, c1.z) + 00610 f2.z * data.w(c2.x, c2.y, c2.z)))); 00611 00612 return ret; 00613 } 00614 00615 //----------------------------------------------------------------------------// 00616 00617 template <class Field_T> 00618 typename Field_T::value_type 00619 CubicGenericFieldInterp<Field_T>::sample(const Field_T &data, 00620 const V3d &vsP) const 00621 { 00622 typedef typename Field_T::value_type Data_T; 00623 00624 // Pixel centers are at .5 coordinates 00625 // NOTE: Don't use contToDisc for this, we're looking for sample 00626 // point locations, not coordinate shifts. 00627 V3d clampedVsP(std::max(0.5, vsP.x), 00628 std::max(0.5, vsP.y), 00629 std::max(0.5, vsP.z)); 00630 FIELD3D_VEC3_T<double> p(clampedVsP - FIELD3D_VEC3_T<double>(0.5)); 00631 00632 const Box3i &dataWindow = data.dataWindow(); 00633 00634 // Lower left corner 00635 V3i c(static_cast<int>(floor(p.x)), 00636 static_cast<int>(floor(p.y)), 00637 static_cast<int>(floor(p.z))); 00638 00639 // Fractions 00640 FIELD3D_VEC3_T<double> t(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 00641 00642 // Clamp the coordinates 00643 int im, jm, km; 00644 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 00645 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 00646 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 00647 int im_1, jm_1, km_1; 00648 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 00649 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 00650 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 00651 int im1, jm1, km1; 00652 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 00653 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 00654 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 00655 int im2, jm2, km2; 00656 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 00657 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 00658 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 00659 00660 Data_T z11 = monotonicCubicInterpolant(data.fastValue(im_1, jm_1, km_1), 00661 data.fastValue(im_1, jm_1, km), 00662 data.fastValue(im_1, jm_1, km1), 00663 data.fastValue(im_1, jm_1, km2), t.z); 00664 Data_T z12 = monotonicCubicInterpolant(data.fastValue(im_1, jm, km_1), 00665 data.fastValue(im_1, jm, km), 00666 data.fastValue(im_1, jm, km1), 00667 data.fastValue(im_1, jm, km2), t.z); 00668 Data_T z13 = monotonicCubicInterpolant(data.fastValue(im_1, jm1, km_1), 00669 data.fastValue(im_1, jm1, km), 00670 data.fastValue(im_1, jm1, km1), 00671 data.fastValue(im_1, jm1, km2), t.z); 00672 Data_T z14 = monotonicCubicInterpolant(data.fastValue(im_1, jm2, km_1), 00673 data.fastValue(im_1, jm2, km), 00674 data.fastValue(im_1, jm2, km1), 00675 data.fastValue(im_1, jm2, km2), t.z); 00676 00677 Data_T z21 = monotonicCubicInterpolant(data.fastValue(im, jm_1, km_1), 00678 data.fastValue(im, jm_1, km), 00679 data.fastValue(im, jm_1, km1), 00680 data.fastValue(im, jm_1, km2), t.z); 00681 Data_T z22 = monotonicCubicInterpolant(data.fastValue(im, jm, km_1), 00682 data.fastValue(im, jm, km), 00683 data.fastValue(im, jm, km1), 00684 data.fastValue(im, jm, km2), t.z); 00685 Data_T z23 = monotonicCubicInterpolant(data.fastValue(im, jm1, km_1), 00686 data.fastValue(im, jm1, km), 00687 data.fastValue(im, jm1, km1), 00688 data.fastValue(im, jm1, km2), t.z); 00689 Data_T z24 = monotonicCubicInterpolant(data.fastValue(im, jm2, km_1), 00690 data.fastValue(im, jm2, km), 00691 data.fastValue(im, jm2, km1), 00692 data.fastValue(im, jm2, km2), t.z); 00693 00694 Data_T z31 = monotonicCubicInterpolant(data.fastValue(im1, jm_1, km_1), 00695 data.fastValue(im1, jm_1, km), 00696 data.fastValue(im1, jm_1, km1), 00697 data.fastValue(im1, jm_1, km2), t.z); 00698 Data_T z32 = monotonicCubicInterpolant(data.fastValue(im1, jm, km_1), 00699 data.fastValue(im1, jm, km), 00700 data.fastValue(im1, jm, km1), 00701 data.fastValue(im1, jm, km2), t.z); 00702 Data_T z33 = monotonicCubicInterpolant(data.fastValue(im1, jm1, km_1), 00703 data.fastValue(im1, jm1, km), 00704 data.fastValue(im1, jm1, km1), 00705 data.fastValue(im1, jm1, km2), t.z); 00706 Data_T z34 = monotonicCubicInterpolant(data.fastValue(im1, jm2, km_1), 00707 data.fastValue(im1, jm2, km), 00708 data.fastValue(im1, jm2, km1), 00709 data.fastValue(im1, jm2, km2), t.z); 00710 00711 Data_T z41 = monotonicCubicInterpolant(data.fastValue(im2, jm_1, km_1), 00712 data.fastValue(im2, jm_1, km), 00713 data.fastValue(im2, jm_1, km1), 00714 data.fastValue(im2, jm_1, km2), t.z); 00715 Data_T z42 = monotonicCubicInterpolant(data.fastValue(im2, jm, km_1), 00716 data.fastValue(im2, jm, km), 00717 data.fastValue(im2, jm, km1), 00718 data.fastValue(im2, jm, km2), t.z); 00719 Data_T z43 = monotonicCubicInterpolant(data.fastValue(im2, jm1, km_1), 00720 data.fastValue(im2, jm1, km), 00721 data.fastValue(im2, jm1, km1), 00722 data.fastValue(im2, jm1, km2), t.z); 00723 Data_T z44 = monotonicCubicInterpolant(data.fastValue(im2, jm2, km_1), 00724 data.fastValue(im2, jm2, km), 00725 data.fastValue(im2, jm2, km1), 00726 data.fastValue(im2, jm2, km2), t.z); 00727 00728 Data_T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 00729 Data_T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 00730 Data_T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 00731 Data_T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 00732 00733 Data_T z0 = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 00734 00735 return z0; 00736 } 00737 00738 //----------------------------------------------------------------------------// 00739 00740 template <class Data_T> 00741 Data_T CubicMACFieldInterp<Data_T>::sample(const MACField<Data_T> &data, 00742 const V3d &vsP) const 00743 { 00744 typedef typename Data_T::BaseType T; 00745 00746 const Box3i &dataWindow = data.dataWindow(); 00747 00748 // Pixel centers are at .5 coordinates 00749 // NOTE: Don't use contToDisc for this, we're looking for sample 00750 // point locations, not coordinate shifts. 00751 00752 Data_T ret; 00753 00754 // X component --- 00755 00756 V3d clampedVsP(std::max(0.5, vsP.x), 00757 std::max(0.5, vsP.y), 00758 std::max(0.5, vsP.z)); 00759 FIELD3D_VEC3_T<double> p(vsP.x, 00760 clampedVsP.y - 0.5, 00761 clampedVsP.z - 0.5); 00762 00763 // Lower left corner 00764 V3i c(static_cast<int>(floor(p.x)), 00765 static_cast<int>(floor(p.y)), 00766 static_cast<int>(floor(p.z))); 00767 00768 FIELD3D_VEC3_T<double> t(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 00769 00770 { 00771 // Clamp the coordinates 00772 int im, jm, km; 00773 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x + 1)); 00774 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 00775 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 00776 int im_1, jm_1, km_1; 00777 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x + 1)); 00778 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 00779 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 00780 int im1, jm1, km1; 00781 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x + 1)); 00782 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 00783 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 00784 int im2, jm2, km2; 00785 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x + 1)); 00786 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 00787 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 00788 00789 T z11 = monotonicCubicInterpolant(data.u(im_1, jm_1, km_1), 00790 data.u(im_1, jm_1, km), 00791 data.u(im_1, jm_1, km1), 00792 data.u(im_1, jm_1, km2), t.z); 00793 T z12 = monotonicCubicInterpolant(data.u(im_1, jm, km_1), 00794 data.u(im_1, jm, km), 00795 data.u(im_1, jm, km1), 00796 data.u(im_1, jm, km2), t.z); 00797 T z13 = monotonicCubicInterpolant(data.u(im_1, jm1, km_1), 00798 data.u(im_1, jm1, km), 00799 data.u(im_1, jm1, km1), 00800 data.u(im_1, jm1, km2), t.z); 00801 T z14 = monotonicCubicInterpolant(data.u(im_1, jm2, km_1), 00802 data.u(im_1, jm2, km), 00803 data.u(im_1, jm2, km1), 00804 data.u(im_1, jm2, km2), t.z); 00805 00806 T z21 = monotonicCubicInterpolant(data.u(im, jm_1, km_1), 00807 data.u(im, jm_1, km), 00808 data.u(im, jm_1, km1), 00809 data.u(im, jm_1, km2), t.z); 00810 T z22 = monotonicCubicInterpolant(data.u(im, jm, km_1), 00811 data.u(im, jm, km), 00812 data.u(im, jm, km1), 00813 data.u(im, jm, km2), t.z); 00814 T z23 = monotonicCubicInterpolant(data.u(im, jm1, km_1), 00815 data.u(im, jm1, km), 00816 data.u(im, jm1, km1), 00817 data.u(im, jm1, km2), t.z); 00818 T z24 = monotonicCubicInterpolant(data.u(im, jm2, km_1), 00819 data.u(im, jm2, km), 00820 data.u(im, jm2, km1), 00821 data.u(im, jm2, km2), t.z); 00822 00823 T z31 = monotonicCubicInterpolant(data.u(im1, jm_1, km_1), 00824 data.u(im1, jm_1, km), 00825 data.u(im1, jm_1, km1), 00826 data.u(im1, jm_1, km2), t.z); 00827 T z32 = monotonicCubicInterpolant(data.u(im1, jm, km_1), 00828 data.u(im1, jm, km), 00829 data.u(im1, jm, km1), 00830 data.u(im1, jm, km2), t.z); 00831 T z33 = monotonicCubicInterpolant(data.u(im1, jm1, km_1), 00832 data.u(im1, jm1, km), 00833 data.u(im1, jm1, km1), 00834 data.u(im1, jm1, km2), t.z); 00835 T z34 = monotonicCubicInterpolant(data.u(im1, jm2, km_1), 00836 data.u(im1, jm2, km), 00837 data.u(im1, jm2, km1), 00838 data.u(im1, jm2, km2), t.z); 00839 00840 T z41 = monotonicCubicInterpolant(data.u(im2, jm_1, km_1), 00841 data.u(im2, jm_1, km), 00842 data.u(im2, jm_1, km1), 00843 data.u(im2, jm_1, km2), t.z); 00844 T z42 = monotonicCubicInterpolant(data.u(im2, jm, km_1), 00845 data.u(im2, jm, km), 00846 data.u(im2, jm, km1), 00847 data.u(im2, jm, km2), t.z); 00848 T z43 = monotonicCubicInterpolant(data.u(im2, jm1, km_1), 00849 data.u(im2, jm1, km), 00850 data.u(im2, jm1, km1), 00851 data.u(im2, jm1, km2), t.z); 00852 T z44 = monotonicCubicInterpolant(data.u(im2, jm2, km_1), 00853 data.u(im2, jm2, km), 00854 data.u(im2, jm2, km1), 00855 data.u(im2, jm2, km2), t.z); 00856 00857 T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 00858 T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 00859 T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 00860 T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 00861 00862 ret.x = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 00863 } 00864 00865 00866 // Y component --- 00867 00868 p.setValue(clampedVsP.x - 0.5, vsP.y , clampedVsP.z - 0.5); 00869 00870 // Lower left corner 00871 c.x = static_cast<int>(floor(p.x)); 00872 c.y = static_cast<int>(floor(p.y)); 00873 c.z = static_cast<int>(floor(p.z)); 00874 00875 t.setValue(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 00876 { 00877 // Clamp the coordinates 00878 int im, jm, km; 00879 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 00880 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y + 1)); 00881 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 00882 int im_1, jm_1, km_1; 00883 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 00884 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y + 1)); 00885 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 00886 int im1, jm1, km1; 00887 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 00888 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y + 1)); 00889 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 00890 int im2, jm2, km2; 00891 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 00892 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y + 1)); 00893 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 00894 00895 T z11 = monotonicCubicInterpolant(data.v(im_1, jm_1, km_1), 00896 data.v(im_1, jm_1, km), 00897 data.v(im_1, jm_1, km1), 00898 data.v(im_1, jm_1, km2), t.z); 00899 T z12 = monotonicCubicInterpolant(data.v(im_1, jm, km_1), 00900 data.v(im_1, jm, km), 00901 data.v(im_1, jm, km1), 00902 data.v(im_1, jm, km2), t.z); 00903 T z13 = monotonicCubicInterpolant(data.v(im_1, jm1, km_1), 00904 data.v(im_1, jm1, km), 00905 data.v(im_1, jm1, km1), 00906 data.v(im_1, jm1, km2), t.z); 00907 T z14 = monotonicCubicInterpolant(data.v(im_1, jm2, km_1), 00908 data.v(im_1, jm2, km), 00909 data.v(im_1, jm2, km1), 00910 data.v(im_1, jm2, km2), t.z); 00911 00912 T z21 = monotonicCubicInterpolant(data.v(im, jm_1, km_1), 00913 data.v(im, jm_1, km), 00914 data.v(im, jm_1, km1), 00915 data.v(im, jm_1, km2), t.z); 00916 T z22 = monotonicCubicInterpolant(data.v(im, jm, km_1), 00917 data.v(im, jm, km), 00918 data.v(im, jm, km1), 00919 data.v(im, jm, km2), t.z); 00920 T z23 = monotonicCubicInterpolant(data.v(im, jm1, km_1), 00921 data.v(im, jm1, km), 00922 data.v(im, jm1, km1), 00923 data.v(im, jm1, km2), t.z); 00924 T z24 = monotonicCubicInterpolant(data.v(im, jm2, km_1), 00925 data.v(im, jm2, km), 00926 data.v(im, jm2, km1), 00927 data.v(im, jm2, km2), t.z); 00928 00929 T z31 = monotonicCubicInterpolant(data.v(im1, jm_1, km_1), 00930 data.v(im1, jm_1, km), 00931 data.v(im1, jm_1, km1), 00932 data.v(im1, jm_1, km2), t.z); 00933 T z32 = monotonicCubicInterpolant(data.v(im1, jm, km_1), 00934 data.v(im1, jm, km), 00935 data.v(im1, jm, km1), 00936 data.v(im1, jm, km2), t.z); 00937 T z33 = monotonicCubicInterpolant(data.v(im1, jm1, km_1), 00938 data.v(im1, jm1, km), 00939 data.v(im1, jm1, km1), 00940 data.v(im1, jm1, km2), t.z); 00941 T z34 = monotonicCubicInterpolant(data.v(im1, jm2, km_1), 00942 data.v(im1, jm2, km), 00943 data.v(im1, jm2, km1), 00944 data.v(im1, jm2, km2), t.z); 00945 00946 T z41 = monotonicCubicInterpolant(data.v(im2, jm_1, km_1), 00947 data.v(im2, jm_1, km), 00948 data.v(im2, jm_1, km1), 00949 data.v(im2, jm_1, km2), t.z); 00950 T z42 = monotonicCubicInterpolant(data.v(im2, jm, km_1), 00951 data.v(im2, jm, km), 00952 data.v(im2, jm, km1), 00953 data.v(im2, jm, km2), t.z); 00954 T z43 = monotonicCubicInterpolant(data.v(im2, jm1, km_1), 00955 data.v(im2, jm1, km), 00956 data.v(im2, jm1, km1), 00957 data.v(im2, jm1, km2), t.z); 00958 T z44 = monotonicCubicInterpolant(data.v(im2, jm2, km_1), 00959 data.v(im2, jm2, km), 00960 data.v(im2, jm2, km1), 00961 data.v(im2, jm2, km2), t.z); 00962 00963 T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 00964 T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 00965 T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 00966 T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 00967 00968 ret.y = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 00969 } 00970 00971 // Z component --- 00972 00973 p.setValue(clampedVsP.x - 0.5 , clampedVsP.y - 0.5, vsP.z); 00974 00975 // Lower left corner 00976 c.x = static_cast<int>(floor(p.x)); 00977 c.y = static_cast<int>(floor(p.y)); 00978 c.z = static_cast<int>(floor(p.z)); 00979 00980 t.setValue(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 00981 { 00982 // Clamp the coordinates 00983 int im, jm, km; 00984 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 00985 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 00986 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z + 1)); 00987 int im_1, jm_1, km_1; 00988 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 00989 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 00990 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z + 1)); 00991 int im1, jm1, km1; 00992 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 00993 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 00994 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z + 1)); 00995 int im2, jm2, km2; 00996 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 00997 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 00998 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z + 1)); 00999 01000 T z11 = monotonicCubicInterpolant(data.w(im_1, jm_1, km_1), 01001 data.w(im_1, jm_1, km), 01002 data.w(im_1, jm_1, km1), 01003 data.w(im_1, jm_1, km2), t.z); 01004 T z12 = monotonicCubicInterpolant(data.w(im_1, jm, km_1), 01005 data.w(im_1, jm, km), 01006 data.w(im_1, jm, km1), 01007 data.w(im_1, jm, km2), t.z); 01008 T z13 = monotonicCubicInterpolant(data.w(im_1, jm1, km_1), 01009 data.w(im_1, jm1, km), 01010 data.w(im_1, jm1, km1), 01011 data.w(im_1, jm1, km2), t.z); 01012 T z14 = monotonicCubicInterpolant(data.w(im_1, jm2, km_1), 01013 data.w(im_1, jm2, km), 01014 data.w(im_1, jm2, km1), 01015 data.w(im_1, jm2, km2), t.z); 01016 01017 T z21 = monotonicCubicInterpolant(data.w(im, jm_1, km_1), 01018 data.w(im, jm_1, km), 01019 data.w(im, jm_1, km1), 01020 data.w(im, jm_1, km2), t.z); 01021 T z22 = monotonicCubicInterpolant(data.w(im, jm, km_1), 01022 data.w(im, jm, km), 01023 data.w(im, jm, km1), 01024 data.w(im, jm, km2), t.z); 01025 T z23 = monotonicCubicInterpolant(data.w(im, jm1, km_1), 01026 data.w(im, jm1, km), 01027 data.w(im, jm1, km1), 01028 data.w(im, jm1, km2), t.z); 01029 T z24 = monotonicCubicInterpolant(data.w(im, jm2, km_1), 01030 data.w(im, jm2, km), 01031 data.w(im, jm2, km1), 01032 data.w(im, jm2, km2), t.z); 01033 01034 T z31 = monotonicCubicInterpolant(data.w(im1, jm_1, km_1), 01035 data.w(im1, jm_1, km), 01036 data.w(im1, jm_1, km1), 01037 data.w(im1, jm_1, km2), t.z); 01038 T z32 = monotonicCubicInterpolant(data.w(im1, jm, km_1), 01039 data.w(im1, jm, km), 01040 data.w(im1, jm, km1), 01041 data.w(im1, jm, km2), t.z); 01042 T z33 = monotonicCubicInterpolant(data.w(im1, jm1, km_1), 01043 data.w(im1, jm1, km), 01044 data.w(im1, jm1, km1), 01045 data.w(im1, jm1, km2), t.z); 01046 T z34 = monotonicCubicInterpolant(data.w(im1, jm2, km_1), 01047 data.w(im1, jm2, km), 01048 data.w(im1, jm2, km1), 01049 data.w(im1, jm2, km2), t.z); 01050 01051 T z41 = monotonicCubicInterpolant(data.w(im2, jm_1, km_1), 01052 data.w(im2, jm_1, km), 01053 data.w(im2, jm_1, km1), 01054 data.w(im2, jm_1, km2), t.z); 01055 T z42 = monotonicCubicInterpolant(data.w(im2, jm, km_1), 01056 data.w(im2, jm, km), 01057 data.w(im2, jm, km1), 01058 data.w(im2, jm, km2), t.z); 01059 T z43 = monotonicCubicInterpolant(data.w(im2, jm1, km_1), 01060 data.w(im2, jm1, km), 01061 data.w(im2, jm1, km1), 01062 data.w(im2, jm1, km2), t.z); 01063 T z44 = monotonicCubicInterpolant(data.w(im2, jm2, km_1), 01064 data.w(im2, jm2, km), 01065 data.w(im2, jm2, km1), 01066 data.w(im2, jm2, km2), t.z); 01067 01068 T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 01069 T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 01070 T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 01071 T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 01072 01073 ret.z = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 01074 } 01075 01076 return ret; 01077 } 01078 01079 //----------------------------------------------------------------------------// 01080 01081 template <class Data_T> 01082 Data_T 01083 ProceduralFieldLookup<Data_T>::sample(const ProceduralField<Data_T> &data, 01084 const V3d &vsP) const 01085 { 01086 V3d voxelScale = V3d(1.0) / data.dataResolution(); 01087 V3d lsP = vsP * voxelScale; 01088 return data.lsSample(lsP); 01089 } 01090 01091 //----------------------------------------------------------------------------// 01092 01093 template <class S, class T> 01094 FIELD3D_VEC3_T<T> operator * (S s, const FIELD3D_VEC3_T<T> vec) 01095 { 01096 return FIELD3D_VEC3_T<T>(vec.x * s, vec.y * s, vec.z * s); 01097 } 01098 01099 //----------------------------------------------------------------------------// 01100 01101 template<class T> 01102 T monotonicCubicInterpolant(const T &f1, const T &f2, const T &f3, const T &f4, 01103 double t) 01104 { 01105 T d_k = T(.5) * (f3 - f1); 01106 T d_k1 = T(.5) * (f4 - f2); 01107 T delta_k = f3 - f2; 01108 01109 if (delta_k == static_cast<T>(0)) { 01110 d_k = static_cast<T>(0); 01111 d_k1 = static_cast<T>(0); 01112 } 01113 01114 T a0 = f2; 01115 T a1 = d_k; 01116 T a2 = (T(3) * delta_k) - (T(2) * d_k) - d_k1; 01117 T a3 = d_k + d_k1 - (T(2) * delta_k); 01118 01119 T t1 = t; 01120 T t2 = t1 * t1; 01121 T t3 = t2 * t1; 01122 01123 return a3 * t3 + a2 * t2 + a1 * t1 + a0; 01124 } 01125 01126 //----------------------------------------------------------------------------// 01127 01129 // References: 01130 // http://en.wikipedia.org/wiki/Monotone_cubic_interpolation 01131 // http://en.wikipedia.org/wiki/Cubic_Hermite_spline 01132 template <class Data_T> 01133 Data_T monotonicCubicInterpolantVec(const Data_T &f1, const Data_T &f2, 01134 const Data_T &f3, const Data_T &f4, 01135 double t) 01136 { 01137 typedef typename Data_T::BaseType T; 01138 01139 Data_T d_k = T(.5) * (f3 - f1); 01140 Data_T d_k1 = T(.5) * (f4 - f2); 01141 Data_T delta_k = f3 - f2; 01142 01143 for (int i = 0; i < 3; i++) { 01144 if (delta_k[i] == static_cast<T>(0)) { 01145 d_k[i] = static_cast<T>(0); 01146 d_k1[i]= static_cast<T>(0); 01147 } 01148 } 01149 01150 Data_T a0 = f2; 01151 Data_T a1 = d_k; 01152 Data_T a2 = (delta_k * T(3)) - (d_k * T(2)) - d_k1; 01153 Data_T a3 = d_k + d_k1 - (delta_k * T(2)); 01154 01155 T t1 = t; 01156 T t2 = t1 * t1; 01157 T t3 = t2 * t1; 01158 01159 return a3 * t3 + a2 * t2 + a1 * t1 + a0; 01160 } 01161 01162 //----------------------------------------------------------------------------// 01163 // Template specializations 01164 //----------------------------------------------------------------------------// 01165 01166 template<> 01167 inline 01168 V3h monotonicCubicInterpolant<V3h>(const V3h &f1, const V3h &f2, 01169 const V3h &f3, const V3h &f4, double t) 01170 { 01171 return monotonicCubicInterpolantVec(f1, f2, f3, f4, t); 01172 } 01173 01174 //----------------------------------------------------------------------------// 01175 01176 template<> 01177 inline 01178 V3f monotonicCubicInterpolant<V3f>(const V3f &f1, const V3f &f2, 01179 const V3f &f3, const V3f &f4, double t) 01180 { 01181 return monotonicCubicInterpolantVec(f1, f2, f3, f4, t); 01182 } 01183 01184 //----------------------------------------------------------------------------// 01185 01186 template<> 01187 inline 01188 V3d monotonicCubicInterpolant<V3d>(const V3d &f1, const V3d &f2, 01189 const V3d &f3, const V3d &f4, double t) 01190 { 01191 return monotonicCubicInterpolantVec(f1, f2, f3, f4, t); 01192 } 01193 01194 //----------------------------------------------------------------------------// 01195 01196 FIELD3D_NAMESPACE_HEADER_CLOSE 01197 01198 //----------------------------------------------------------------------------// 01199 01200 #endif // Include guard