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 00057 #include "ns.h" 00058 00059 FIELD3D_NAMESPACE_OPEN 00060 00061 //----------------------------------------------------------------------------// 00062 // FieldInterp 00063 //----------------------------------------------------------------------------// 00064 00071 template <class Data_T> 00072 class FieldInterp : public RefBase 00073 { 00074 public: 00075 00076 // Typedefs ------------------------------------------------------------------ 00077 00078 typedef Data_T value_type; 00079 typedef boost::intrusive_ptr<FieldInterp> Ptr; 00080 00081 // RTTI replacement ---------------------------------------------------------- 00082 00083 typedef FieldInterp class_type; 00084 DEFINE_FIELD_RTTI_ABSTRACT_CLASS; 00085 00086 static const char *staticClassName() 00087 { 00088 return "FieldInterp"; 00089 } 00090 00091 static const char* classType() 00092 { 00093 return ms_classType.name(); 00094 } 00095 00096 // Ctor, dtor ---------------------------------------------------------------- 00097 00098 virtual ~FieldInterp() 00099 { } 00100 00101 // Main methods -------------------------------------------------------------- 00102 00103 virtual Data_T sample(const Field<Data_T> &data, const V3d &vsP) const = 0; 00104 00105 private: 00106 00107 // Static data members ------------------------------------------------------- 00108 00109 static TemplatedFieldType<FieldInterp<Data_T> > ms_classType; 00110 00111 // Typedefs ------------------------------------------------------------------ 00112 00114 typedef RefBase base; 00115 00116 }; 00117 00118 //----------------------------------------------------------------------------// 00119 // Static member instantiation 00120 //----------------------------------------------------------------------------// 00121 00122 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(FieldInterp); 00123 00124 //----------------------------------------------------------------------------// 00125 // LinearFieldInterp 00126 //----------------------------------------------------------------------------// 00127 00128 /* \class LinearFieldInterp 00129 \ingroup field 00130 \brief Basic linear interpolator using voxel access through Field base class 00131 */ 00132 00133 //----------------------------------------------------------------------------// 00134 00135 template <class Data_T> 00136 class LinearFieldInterp : public FieldInterp<Data_T> 00137 { 00138 public: 00139 00140 // Typedefs ------------------------------------------------------------------ 00141 00142 typedef Data_T value_type; 00143 typedef boost::intrusive_ptr<LinearFieldInterp> Ptr; 00144 00145 // RTTI replacement ---------------------------------------------------------- 00146 00147 typedef LinearFieldInterp class_type; 00148 DEFINE_FIELD_RTTI_CONCRETE_CLASS; 00149 00150 static const char *staticClassName() 00151 { 00152 return "LinearFieldInterp"; 00153 } 00154 00155 static const char* classType() 00156 { 00157 return ms_classType.name(); 00158 } 00159 00160 // From FieldInterp ---------------------------------------------------------- 00161 00162 virtual Data_T sample(const Field<Data_T> &data, const V3d &vsP) const; 00163 00164 private: 00165 00166 // Static data members ------------------------------------------------------- 00167 00168 static TemplatedFieldType<LinearFieldInterp<Data_T> > ms_classType; 00169 00170 // Typedefs ------------------------------------------------------------------ 00171 00173 typedef FieldInterp<Data_T> base; 00174 00175 }; 00176 00177 //----------------------------------------------------------------------------// 00178 // Static data member instantiation 00179 //----------------------------------------------------------------------------// 00180 00181 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(LinearFieldInterp); 00182 00183 //----------------------------------------------------------------------------// 00184 // CubicFieldInterp 00185 //----------------------------------------------------------------------------// 00186 00187 /* \class CubicFieldInterp 00188 \ingroup field 00189 \brief Basic cubic interpolator using voxel access through Field base class 00190 */ 00191 00192 //----------------------------------------------------------------------------// 00193 00194 template <class Data_T> 00195 class CubicFieldInterp : public FieldInterp<Data_T> 00196 { 00197 public: 00198 00199 // Typedefs ------------------------------------------------------------------ 00200 00201 typedef Data_T value_type; 00202 typedef boost::intrusive_ptr<CubicFieldInterp> Ptr; 00203 00204 // RTTI replacement ---------------------------------------------------------- 00205 00206 typedef CubicFieldInterp class_type; 00207 DEFINE_FIELD_RTTI_CONCRETE_CLASS; 00208 00209 static const char *staticClassName() 00210 { 00211 return "CubicFieldInterp"; 00212 } 00213 00214 static const char* classType() 00215 { 00216 return ms_classType.name(); 00217 } 00218 00219 // From FieldInterp ---------------------------------------------------------- 00220 00221 virtual Data_T sample(const Field<Data_T> &data, const V3d &vsP) const; 00222 00223 private: 00224 00225 // Static data members ------------------------------------------------------- 00226 00227 static TemplatedFieldType<CubicFieldInterp<Data_T> > ms_classType; 00228 00229 // Typedefs ------------------------------------------------------------------ 00230 00232 typedef FieldInterp<Data_T> base; 00233 }; 00234 00235 //----------------------------------------------------------------------------// 00236 // Static data member instantiation 00237 //----------------------------------------------------------------------------// 00238 00239 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(CubicFieldInterp); 00240 00241 //----------------------------------------------------------------------------// 00242 // LinearGenericFieldInterp 00243 //----------------------------------------------------------------------------// 00244 00245 /* \class LinearGenericFieldInterp 00246 \ingroup field 00247 \brief Linear interpolator optimized for fields with a fastValue function 00248 */ 00249 00250 //----------------------------------------------------------------------------// 00251 00252 template <class Field_T> 00253 class LinearGenericFieldInterp : public RefBase 00254 { 00255 public: 00256 00257 // Typedefs ------------------------------------------------------------------ 00258 00259 typedef typename Field_T::value_type value_type; 00260 typedef boost::intrusive_ptr<LinearGenericFieldInterp> Ptr; 00261 00262 // RTTI replacement ---------------------------------------------------------- 00263 00264 typedef LinearGenericFieldInterp class_type; 00265 DEFINE_FIELD_RTTI_CONCRETE_CLASS; 00266 00267 static const char *staticClassName() 00268 { 00269 return "LinearGenericFieldInterp"; 00270 } 00271 00272 static const char* classType() 00273 { 00274 return ms_classType.name(); 00275 } 00276 00277 // Main methods -------------------------------------------------------------- 00278 00279 value_type sample(const Field_T &data, const V3d &vsP) const; 00280 00281 private: 00282 00283 // Static data members ------------------------------------------------------- 00284 00285 static TemplatedFieldType<LinearGenericFieldInterp<Field_T> > ms_classType; 00286 00287 // Typedefs ------------------------------------------------------------------ 00288 00290 typedef RefBase base; 00291 00292 }; 00293 00294 //----------------------------------------------------------------------------// 00295 // Static data member instantiation 00296 //----------------------------------------------------------------------------// 00297 00298 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(LinearGenericFieldInterp); 00299 00300 //----------------------------------------------------------------------------// 00301 // LinearMACFieldInterp 00302 //----------------------------------------------------------------------------// 00303 00304 /* \class LinearMACFieldInterp 00305 \ingroup field 00306 \brief Linear interpolator optimized for the MAC fields 00307 */ 00308 00309 //----------------------------------------------------------------------------// 00310 00311 template <class Data_T> 00312 class LinearMACFieldInterp : public RefBase 00313 { 00314 public: 00315 00316 // Typedefs ------------------------------------------------------------------ 00317 00318 typedef Data_T value_type; 00319 typedef boost::intrusive_ptr<LinearMACFieldInterp> Ptr; 00320 00321 // RTTI replacement ---------------------------------------------------------- 00322 00323 typedef LinearMACFieldInterp class_type; 00324 DEFINE_FIELD_RTTI_CONCRETE_CLASS; 00325 00326 static const char *staticClassName() 00327 { 00328 return "LinearMACFieldInterp"; 00329 } 00330 00332 static const char* classType() 00333 { 00334 return ms_classType.name(); 00335 } 00336 00337 // Main methods -------------------------------------------------------------- 00338 00339 Data_T sample(const MACField<Data_T> &data, const V3d &vsP) const; 00340 00341 double sample(const MACField<Data_T> &data, 00342 const MACComponent &comp, 00343 const V3d &vsP) const; 00344 00345 private: 00346 00347 // Static data members ------------------------------------------------------- 00348 00349 static TemplatedFieldType<LinearMACFieldInterp<Data_T> > ms_classType; 00350 00351 // Typedefs ------------------------------------------------------------------ 00352 00354 typedef RefBase base; 00355 }; 00356 00357 //----------------------------------------------------------------------------// 00358 // Static data member instantiation 00359 //----------------------------------------------------------------------------// 00360 00361 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(LinearMACFieldInterp); 00362 00363 //----------------------------------------------------------------------------// 00364 // CubicGenericFieldInterp 00365 //----------------------------------------------------------------------------// 00366 00367 /* \class CubicGenericFieldInterp 00368 \ingroup field 00369 \brief Cubic interpolator optimized for fields with a fastValue function 00370 */ 00371 00372 //----------------------------------------------------------------------------// 00373 00374 template <class Field_T> 00375 class CubicGenericFieldInterp : public RefBase 00376 { 00377 public: 00378 00379 // Typedefs ------------------------------------------------------------------ 00380 00381 typedef typename Field_T::value_type value_type; 00382 typedef boost::intrusive_ptr<CubicGenericFieldInterp> Ptr; 00383 00384 // RTTI replacement ---------------------------------------------------------- 00385 00386 typedef CubicGenericFieldInterp class_type; 00387 DEFINE_FIELD_RTTI_CONCRETE_CLASS; 00388 00389 static const char *staticClassName() 00390 { 00391 return "CubicGenericFieldInterp"; 00392 } 00393 00394 static const char* classType() 00395 { 00396 return ms_classType.name(); 00397 } 00398 00399 // Main methods -------------------------------------------------------------- 00400 00401 value_type sample(const Field_T &data, const V3d &vsP) const; 00402 00403 00404 private: 00405 00406 // Static data members ------------------------------------------------------- 00407 00408 static TemplatedFieldType<CubicGenericFieldInterp<Field_T> > ms_classType; 00409 00410 // Typedefs ------------------------------------------------------------------ 00411 00413 typedef RefBase base; 00414 }; 00415 00416 //----------------------------------------------------------------------------// 00417 // Static data member instantiation 00418 //----------------------------------------------------------------------------// 00419 00420 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(CubicGenericFieldInterp); 00421 00422 //----------------------------------------------------------------------------// 00423 // CubicMACFieldInterp 00424 //----------------------------------------------------------------------------// 00425 00426 /* \class CubicMACFieldInterp 00427 \ingroup field 00428 \brief Linear interpolator optimized for MAC fields 00429 */ 00430 00431 //----------------------------------------------------------------------------// 00432 00433 template <class Data_T> 00434 class CubicMACFieldInterp : public RefBase 00435 { 00436 public: 00437 00438 // Typedefs ------------------------------------------------------------------ 00439 00440 typedef Data_T value_type; 00441 typedef boost::intrusive_ptr<CubicMACFieldInterp> Ptr; 00442 00443 // RTTI replacement ---------------------------------------------------------- 00444 00445 typedef CubicMACFieldInterp class_type; 00446 DEFINE_FIELD_RTTI_CONCRETE_CLASS; 00447 00448 static const char *staticClassName() 00449 { 00450 return "CubicMACFieldInterp"; 00451 } 00452 00453 static const char* classType() 00454 { 00455 return CubicMACFieldInterp<Data_T>::ms_classType.name(); 00456 } 00457 00458 // Main methods -------------------------------------------------------------- 00459 00460 Data_T sample(const MACField<Data_T> &data, const V3d &vsP) const; 00461 00462 private: 00463 00464 // Static data members ------------------------------------------------------- 00465 00466 static TemplatedFieldType<CubicMACFieldInterp<Data_T> > ms_classType; 00467 00468 // Typedefs ------------------------------------------------------------------ 00469 00471 typedef RefBase base; 00472 }; 00473 00474 //----------------------------------------------------------------------------// 00475 // Static data member instantiation 00476 //----------------------------------------------------------------------------// 00477 00478 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(CubicMACFieldInterp); 00479 00480 //----------------------------------------------------------------------------// 00481 // ProceduralFieldLookup 00482 //----------------------------------------------------------------------------// 00483 00484 /* \class ProceduralFieldLookup 00485 \ingroup field 00486 \brief "Interpolator" for procedural fields - point samples instead of 00487 interpolating. 00488 */ 00489 00490 //----------------------------------------------------------------------------// 00491 00492 template <class Data_T> 00493 class ProceduralFieldLookup : public RefBase 00494 { 00495 public: 00496 00497 // Typedefs ------------------------------------------------------------------ 00498 00499 typedef Data_T value_type; 00500 typedef boost::intrusive_ptr<ProceduralFieldLookup> Ptr; 00501 00502 // RTTI replacement ---------------------------------------------------------- 00503 00504 typedef ProceduralFieldLookup class_type; 00505 DEFINE_FIELD_RTTI_CONCRETE_CLASS; 00506 00507 static const char *staticClassName() 00508 { 00509 return "ProceduralFieldLookup"; 00510 } 00511 00512 static const char* classType() 00513 { 00514 return ProceduralFieldLookup<Data_T>::ms_classType.name(); 00515 } 00516 00517 // Main methods -------------------------------------------------------------- 00518 00519 Data_T sample(const ProceduralField<Data_T> &data, const V3d &vsP) const; 00520 00521 private: 00522 00523 // Static data members ------------------------------------------------------- 00524 00525 static TemplatedFieldType<ProceduralFieldLookup<Data_T> > ms_classType; 00526 00527 // Typedefs ------------------------------------------------------------------ 00528 00530 typedef RefBase base; 00531 00532 }; 00533 00534 //----------------------------------------------------------------------------// 00535 // Static data member instantiation 00536 //----------------------------------------------------------------------------// 00537 00538 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(ProceduralFieldLookup); 00539 00540 //----------------------------------------------------------------------------// 00541 // Interpolation functions 00542 //----------------------------------------------------------------------------// 00543 00545 template <class Data_T> 00546 Data_T wsSample(const typename Field<Data_T>::Ptr f, 00547 const FieldInterp<Data_T> &interp, 00548 const V3d &wsP) 00549 { 00550 V3d vsP; 00551 f->mapping()->worldToVoxel(wsP, vsP); 00552 return interp.sample(*f, vsP); 00553 } 00554 00555 //----------------------------------------------------------------------------// 00556 // Interpolation helper functions 00557 //----------------------------------------------------------------------------// 00558 00560 bool isPointInField(const FieldRes::Ptr f, const V3d &wsP); 00561 00562 //----------------------------------------------------------------------------// 00563 00566 bool isLegalVoxelCoord(const V3d &vsP, const Box3d &vsDataWindow); 00567 00568 //----------------------------------------------------------------------------// 00569 // Math functions 00570 //----------------------------------------------------------------------------// 00571 00574 template <class S, class T> 00575 FIELD3D_VEC3_T<T> operator * (S s, const FIELD3D_VEC3_T<T> vec); 00576 00577 //----------------------------------------------------------------------------// 00578 // Interpolants 00579 //----------------------------------------------------------------------------// 00580 00585 template <class Data_T> 00586 Data_T monotonicCubicInterpolant(const Data_T &f1, const Data_T &f2, 00587 const Data_T &f3, const Data_T &f4, 00588 double t); 00589 00590 //----------------------------------------------------------------------------// 00591 00596 template <class Data_T> 00597 Data_T monotonicCubicInterpolantVec(const Data_T &f1, const Data_T &f2, 00598 const Data_T &f3, const Data_T &f4, 00599 double t); 00600 00601 //----------------------------------------------------------------------------// 00602 // Implementations 00603 //----------------------------------------------------------------------------// 00604 00605 template <class Data_T> 00606 Data_T LinearFieldInterp<Data_T>::sample(const Field<Data_T> &data, 00607 const V3d &vsP) const 00608 { 00609 // Voxel centers are at .5 coordinates 00610 // NOTE: Don't use contToDisc for this, we're looking for sample 00611 // point locations, not coordinate shifts. 00612 FIELD3D_VEC3_T<double> p(vsP - FIELD3D_VEC3_T<double>(0.5)); 00613 00614 // Lower left corner 00615 V3i c1(static_cast<int>(floor(p.x)), 00616 static_cast<int>(floor(p.y)), 00617 static_cast<int>(floor(p.z))); 00618 // Upper right corner 00619 V3i c2(c1 + V3i(1)); 00620 // C1 fractions 00621 FIELD3D_VEC3_T<double> f1(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00622 // C2 fraction 00623 FIELD3D_VEC3_T<double> f2(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00624 00625 // Clamp the indexing coordinates 00626 if (true) { 00627 const Box3i &dataWindow = data.dataWindow(); 00628 c1.x = std::max(dataWindow.min.x, std::min(c1.x, dataWindow.max.x)); 00629 c2.x = std::max(dataWindow.min.x, std::min(c2.x, dataWindow.max.x)); 00630 c1.y = std::max(dataWindow.min.y, std::min(c1.y, dataWindow.max.y)); 00631 c2.y = std::max(dataWindow.min.y, std::min(c2.y, dataWindow.max.y)); 00632 c1.z = std::max(dataWindow.min.z, std::min(c1.z, dataWindow.max.z)); 00633 c2.z = std::max(dataWindow.min.z, std::min(c2.z, dataWindow.max.z)); 00634 } 00635 00636 return static_cast<Data_T> 00637 (f1.x * (f1.y * (f1.z * data.value(c1.x, c1.y, c1.z) + 00638 f2.z * data.value(c1.x, c1.y, c2.z)) + 00639 f2.y * (f1.z * data.value(c1.x, c2.y, c1.z) + 00640 f2.z * data.value(c1.x, c2.y, c2.z))) + 00641 f2.x * (f1.y * (f1.z * data.value(c2.x, c1.y, c1.z) + 00642 f2.z * data.value(c2.x, c1.y, c2.z)) + 00643 f2.y * (f1.z * data.value(c2.x, c2.y, c1.z) + 00644 f2.z * data.value(c2.x, c2.y, c2.z)))); 00645 00646 } 00647 00648 //----------------------------------------------------------------------------// 00649 00650 template <class Data_T> 00651 Data_T CubicFieldInterp<Data_T>::sample(const Field<Data_T> &data, 00652 const V3d &vsP) const 00653 { 00654 // Voxel centers are at .5 coordinates 00655 // NOTE: Don't use contToDisc for this, we're looking for sample 00656 // point locations, not coordinate shifts. 00657 V3d clampedVsP(std::max(0.5, vsP.x), 00658 std::max(0.5, vsP.y), 00659 std::max(0.5, vsP.z)); 00660 FIELD3D_VEC3_T<double> p(clampedVsP - FIELD3D_VEC3_T<double>(0.5)); 00661 00662 // Lower left corner 00663 V3i c(static_cast<int>(floor(p.x)), 00664 static_cast<int>(floor(p.y)), 00665 static_cast<int>(floor(p.z))); 00666 00667 // Fractions 00668 FIELD3D_VEC3_T<double> t(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 00669 00670 const Box3i &dataWindow = data.dataWindow(); 00671 00672 // Clamp the coordinates 00673 int im, jm, km; 00674 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 00675 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 00676 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 00677 int im_1, jm_1, km_1; 00678 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 00679 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 00680 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 00681 int im1, jm1, km1; 00682 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 00683 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 00684 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 00685 int im2, jm2, km2; 00686 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 00687 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 00688 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 00689 00690 // interpolate 16 lines in z: 00691 Data_T z11 = monotonicCubicInterpolant(data.value(im_1, jm_1, km_1), 00692 data.value(im_1, jm_1, km), 00693 data.value(im_1, jm_1, km1), 00694 data.value(im_1, jm_1, km2), t.z); 00695 Data_T z12 = monotonicCubicInterpolant(data.value(im_1, jm, km_1), 00696 data.value(im_1, jm, km), 00697 data.value(im_1, jm, km1), 00698 data.value(im_1, jm, km2), t.z); 00699 Data_T z13 = monotonicCubicInterpolant(data.value(im_1, jm1, km_1), 00700 data.value(im_1, jm1, km), 00701 data.value(im_1, jm1, km1), 00702 data.value(im_1, jm1, km2), t.z); 00703 Data_T z14 = monotonicCubicInterpolant(data.value(im_1, jm2, km_1), 00704 data.value(im_1, jm2, km), 00705 data.value(im_1, jm2, km1), 00706 data.value(im_1, jm2, km2), t.z); 00707 00708 Data_T z21 = monotonicCubicInterpolant(data.value(im, jm_1, km_1), 00709 data.value(im, jm_1, km), 00710 data.value(im, jm_1, km1), 00711 data.value(im, jm_1, km2), t.z); 00712 Data_T z22 = monotonicCubicInterpolant(data.value(im, jm, km_1), 00713 data.value(im, jm, km), 00714 data.value(im, jm, km1), 00715 data.value(im, jm, km2), t.z); 00716 Data_T z23 = monotonicCubicInterpolant(data.value(im, jm1, km_1), 00717 data.value(im, jm1, km), 00718 data.value(im, jm1, km1), 00719 data.value(im, jm1, km2), t.z); 00720 Data_T z24 = monotonicCubicInterpolant(data.value(im, jm2, km_1), 00721 data.value(im, jm2, km), 00722 data.value(im, jm2, km1), 00723 data.value(im, jm2, km2), t.z); 00724 00725 Data_T z31 = monotonicCubicInterpolant(data.value(im1, jm_1, km_1), 00726 data.value(im1, jm_1, km), 00727 data.value(im1, jm_1, km1), 00728 data.value(im1, jm_1, km2), t.z); 00729 Data_T z32 = monotonicCubicInterpolant(data.value(im1, jm, km_1), 00730 data.value(im1, jm, km), 00731 data.value(im1, jm, km1), 00732 data.value(im1, jm, km2), t.z); 00733 Data_T z33 = monotonicCubicInterpolant(data.value(im1, jm1, km_1), 00734 data.value(im1, jm1, km), 00735 data.value(im1, jm1, km1), 00736 data.value(im1, jm1, km2), t.z); 00737 Data_T z34 = monotonicCubicInterpolant(data.value(im1, jm2, km_1), 00738 data.value(im1, jm2, km), 00739 data.value(im1, jm2, km1), 00740 data.value(im1, jm2, km2), t.z); 00741 00742 Data_T z41 = monotonicCubicInterpolant(data.value(im2, jm_1, km_1), 00743 data.value(im2, jm_1, km), 00744 data.value(im2, jm_1, km1), 00745 data.value(im2, jm_1, km2), t.z); 00746 Data_T z42 = monotonicCubicInterpolant(data.value(im2, jm, km_1), 00747 data.value(im2, jm, km), 00748 data.value(im2, jm, km1), 00749 data.value(im2, jm, km2), t.z); 00750 Data_T z43 = monotonicCubicInterpolant(data.value(im2, jm1, km_1), 00751 data.value(im2, jm1, km), 00752 data.value(im2, jm1, km1), 00753 data.value(im2, jm1, km2), t.z); 00754 Data_T z44 = monotonicCubicInterpolant(data.value(im2, jm2, km_1), 00755 data.value(im2, jm2, km), 00756 data.value(im2, jm2, km1), 00757 data.value(im2, jm2, km2), t.z); 00758 00759 Data_T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 00760 Data_T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 00761 Data_T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 00762 Data_T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 00763 00764 Data_T z0 = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 00765 00766 return z0; 00767 } 00768 00769 //----------------------------------------------------------------------------// 00770 00771 template <class Field_T> 00772 typename Field_T::value_type 00773 LinearGenericFieldInterp<Field_T>::sample(const Field_T &data, 00774 const V3d &vsP) const 00775 { 00776 typedef typename Field_T::value_type Data_T; 00777 00778 // Pixel centers are at .5 coordinates 00779 // NOTE: Don't use contToDisc for this, we're looking for sample 00780 // point locations, not coordinate shifts. 00781 FIELD3D_VEC3_T<double> p(vsP - FIELD3D_VEC3_T<double>(0.5)); 00782 00783 // Lower left corner 00784 V3i c1(static_cast<int>(floor(p.x)), 00785 static_cast<int>(floor(p.y)), 00786 static_cast<int>(floor(p.z))); 00787 // Upper right corner 00788 V3i c2(c1 + V3i(1)); 00789 // C1 fractions 00790 FIELD3D_VEC3_T<double> f1(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00791 // C2 fraction 00792 FIELD3D_VEC3_T<double> f2(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00793 00794 const Box3i &dataWindow = data.dataWindow(); 00795 00796 // Clamp the coordinates 00797 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 00798 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 00799 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 00800 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 00801 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 00802 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 00803 00804 return static_cast<Data_T> 00805 (f1.x * (f1.y * (f1.z * data.fastValue(c1.x, c1.y, c1.z) + 00806 f2.z * data.fastValue(c1.x, c1.y, c2.z)) + 00807 f2.y * (f1.z * data.fastValue(c1.x, c2.y, c1.z) + 00808 f2.z * data.fastValue(c1.x, c2.y, c2.z))) + 00809 f2.x * (f1.y * (f1.z * data.fastValue(c2.x, c1.y, c1.z) + 00810 f2.z * data.fastValue(c2.x, c1.y, c2.z)) + 00811 f2.y * (f1.z * data.fastValue(c2.x, c2.y, c1.z) + 00812 f2.z * data.fastValue(c2.x, c2.y, c2.z)))); 00813 } 00814 00815 //----------------------------------------------------------------------------// 00816 00817 template <class Data_T> 00818 Data_T LinearMACFieldInterp<Data_T>::sample(const MACField<Data_T> &data, 00819 const V3d &vsP) const 00820 { 00821 // Pixel centers are at .5 coordinates 00822 // NOTE: Don't use contToDisc for this, we're looking for sample 00823 // point locations, not coordinate shifts. 00824 00825 const Box3i &dataWindow = data.dataWindow(); 00826 00827 Data_T ret; 00828 00829 FIELD3D_VEC3_T<double> p(vsP.x , vsP.y - 0.5, vsP.z - 0.5); 00830 00831 // X component --- 00832 00833 // Lower left corner 00834 V3i c1(static_cast<int>(floor(p.x)), 00835 static_cast<int>(floor(p.y)), 00836 static_cast<int>(floor(p.z))); 00837 00838 // Upper right corner 00839 V3i c2(c1 + V3i(1)); 00840 00841 // C1 fractions 00842 FIELD3D_VEC3_T<double> f1(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00843 // C2 fraction 00844 FIELD3D_VEC3_T<double> f2(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00845 00846 // Clamp the coordinates 00847 c1.x = std::min(dataWindow.max.x + 1, std::max(dataWindow.min.x, c1.x)); 00848 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 00849 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 00850 c2.x = std::min(dataWindow.max.x + 1, std::max(dataWindow.min.x, c2.x)); 00851 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 00852 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 00853 00854 ret.x = (f1.x * (f1.y * (f1.z * data.u(c1.x, c1.y, c1.z) + 00855 f2.z * data.u(c1.x, c1.y, c2.z)) + 00856 f2.y * (f1.z * data.u(c1.x, c2.y, c1.z) + 00857 f2.z * data.u(c1.x, c2.y, c2.z))) + 00858 f2.x * (f1.y * (f1.z * data.u(c2.x, c1.y, c1.z) + 00859 f2.z * data.u(c2.x, c1.y, c2.z)) + 00860 f2.y * (f1.z * data.u(c2.x, c2.y, c1.z) + 00861 f2.z * data.u(c2.x, c2.y, c2.z)))); 00862 00863 // Y component --- 00864 00865 p.setValue(vsP.x - 0.5, vsP.y , vsP.z - 0.5); 00866 00867 // Lower left corner 00868 c1.x = static_cast<int>(floor(p.x )); 00869 c1.y = static_cast<int>(floor(p.y )); 00870 c1.z = static_cast<int>(floor(p.z )); 00871 00872 // Upper right corner 00873 c2.x = c1.x + 1; 00874 c2.y = c1.y + 1; 00875 c2.z = c1.z + 1; 00876 00877 // C1 fractions 00878 f1.setValue(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00879 // C2 fraction 00880 f2.setValue(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00881 00882 // Clamp the coordinates 00883 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 00884 c1.y = std::min(dataWindow.max.y + 1, std::max(dataWindow.min.y, c1.y)); 00885 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 00886 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 00887 c2.y = std::min(dataWindow.max.y + 1, std::max(dataWindow.min.y, c2.y)); 00888 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 00889 00890 ret.y = (f1.x * (f1.y * (f1.z * data.v(c1.x, c1.y, c1.z) + 00891 f2.z * data.v(c1.x, c1.y, c2.z)) + 00892 f2.y * (f1.z * data.v(c1.x, c2.y, c1.z) + 00893 f2.z * data.v(c1.x, c2.y, c2.z))) + 00894 f2.x * (f1.y * (f1.z * data.v(c2.x, c1.y, c1.z) + 00895 f2.z * data.v(c2.x, c1.y, c2.z)) + 00896 f2.y * (f1.z * data.v(c2.x, c2.y, c1.z) + 00897 f2.z * data.v(c2.x, c2.y, c2.z)))); 00898 00899 // Z component --- 00900 00901 p.setValue(vsP.x - 0.5 , vsP.y - 0.5, vsP.z); 00902 00903 // Lower left corner 00904 c1.x = static_cast<int>(floor(p.x )); 00905 c1.y = static_cast<int>(floor(p.y )); 00906 c1.z = static_cast<int>(floor(p.z )); 00907 00908 // Upper right corner 00909 c2.x = c1.x + 1; 00910 c2.y = c1.y + 1; 00911 c2.z = c1.z + 1; 00912 00913 // C1 fractions 00914 f1.setValue(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00915 // C2 fraction 00916 f2.setValue(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00917 00918 // Clamp the coordinates 00919 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 00920 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 00921 c1.z = std::min(dataWindow.max.z + 1, std::max(dataWindow.min.z, c1.z)); 00922 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 00923 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 00924 c2.z = std::min(dataWindow.max.z + 1, std::max(dataWindow.min.z, c2.z)); 00925 00926 ret.z = (f1.x * (f1.y * (f1.z * data.w(c1.x, c1.y, c1.z) + 00927 f2.z * data.w(c1.x, c1.y, c2.z)) + 00928 f2.y * (f1.z * data.w(c1.x, c2.y, c1.z) + 00929 f2.z * data.w(c1.x, c2.y, c2.z))) + 00930 f2.x * (f1.y * (f1.z * data.w(c2.x, c1.y, c1.z) + 00931 f2.z * data.w(c2.x, c1.y, c2.z)) + 00932 f2.y * (f1.z * data.w(c2.x, c2.y, c1.z) + 00933 f2.z * data.w(c2.x, c2.y, c2.z)))); 00934 00935 return ret; 00936 } 00937 00938 //----------------------------------------------------------------------------// 00939 00940 template <class Data_T> 00941 double LinearMACFieldInterp<Data_T>::sample(const MACField<Data_T> &data, 00942 const MACComponent &comp, 00943 const V3d &vsP) const 00944 { 00945 // Pixel centers are at .5 coordinates 00946 // NOTE: Don't use contToDisc for this, we're looking for sample 00947 // point locations, not coordinate shifts. 00948 00949 const Box3i &dataWindow = data.dataWindow(); 00950 00951 double ret = 0.0; 00952 FIELD3D_VEC3_T<double> p; 00953 V3i c1, c2; 00954 FIELD3D_VEC3_T<double> f1; 00955 FIELD3D_VEC3_T<double> f2; 00956 00957 switch(comp) { 00958 // U component --- 00959 case MACCompU: 00960 { 00961 p.setValue<>(vsP.x, vsP.y-0.5, vsP.z-0.5); 00962 00963 // Lower left corner 00964 c1.x = static_cast<int>(floor(p.x)); 00965 c1.y = static_cast<int>(floor(p.y)); 00966 c1.z = static_cast<int>(floor(p.z)); 00967 00968 // Upper right corner 00969 c2.x = c1.x + 1; 00970 c2.y = c1.y + 1; 00971 c2.z = c1.z + 1; 00972 00973 // C1 fractions 00974 f1.setValue(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 00975 // C2 fraction 00976 f2.setValue(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 00977 00978 // Clamp the coordinates 00979 c1.x = std::min(dataWindow.max.x + 1, std::max(dataWindow.min.x, c1.x)); 00980 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 00981 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 00982 c2.x = std::min(dataWindow.max.x + 1, std::max(dataWindow.min.x, c2.x)); 00983 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 00984 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 00985 00986 ret = (f1.x * (f1.y * (f1.z * data.u(c1.x, c1.y, c1.z) + 00987 f2.z * data.u(c1.x, c1.y, c2.z)) + 00988 f2.y * (f1.z * data.u(c1.x, c2.y, c1.z) + 00989 f2.z * data.u(c1.x, c2.y, c2.z))) + 00990 f2.x * (f1.y * (f1.z * data.u(c2.x, c1.y, c1.z) + 00991 f2.z * data.u(c2.x, c1.y, c2.z)) + 00992 f2.y * (f1.z * data.u(c2.x, c2.y, c1.z) + 00993 f2.z * data.u(c2.x, c2.y, c2.z)))); 00994 break; 00995 } 00996 // Y component --- 00997 case MACCompV: 00998 { 00999 p.setValue(vsP.x-0.5, vsP.y, vsP.z-0.5); 01000 01001 // Lower left corner 01002 c1.x = static_cast<int>(floor(p.x )); 01003 c1.y = static_cast<int>(floor(p.y )); 01004 c1.z = static_cast<int>(floor(p.z )); 01005 01006 // Upper right corner 01007 c2.x = c1.x + 1; 01008 c2.y = c1.y + 1; 01009 c2.z = c1.z + 1; 01010 01011 // C1 fractions 01012 f1.setValue(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 01013 // C2 fraction 01014 f2.setValue(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 01015 01016 // Clamp the coordinates 01017 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 01018 c1.y = std::min(dataWindow.max.y + 1, std::max(dataWindow.min.y, c1.y)); 01019 c1.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c1.z)); 01020 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 01021 c2.y = std::min(dataWindow.max.y + 1, std::max(dataWindow.min.y, c2.y)); 01022 c2.z = std::min(dataWindow.max.z, std::max(dataWindow.min.z, c2.z)); 01023 01024 ret = (f1.x * (f1.y * (f1.z * data.v(c1.x, c1.y, c1.z) + 01025 f2.z * data.v(c1.x, c1.y, c2.z)) + 01026 f2.y * (f1.z * data.v(c1.x, c2.y, c1.z) + 01027 f2.z * data.v(c1.x, c2.y, c2.z))) + 01028 f2.x * (f1.y * (f1.z * data.v(c2.x, c1.y, c1.z) + 01029 f2.z * data.v(c2.x, c1.y, c2.z)) + 01030 f2.y * (f1.z * data.v(c2.x, c2.y, c1.z) + 01031 f2.z * data.v(c2.x, c2.y, c2.z)))); 01032 break; 01033 } 01034 // W component --- 01035 case MACCompW: 01036 { 01037 p.setValue(vsP.x-0.5, vsP.y-0.5, vsP.z); 01038 01039 // Lower left corner 01040 c1.x = static_cast<int>(floor(p.x )); 01041 c1.y = static_cast<int>(floor(p.y )); 01042 c1.z = static_cast<int>(floor(p.z )); 01043 01044 // Upper right corner 01045 c2.x = c1.x + 1; 01046 c2.y = c1.y + 1; 01047 c2.z = c1.z + 1; 01048 01049 // C1 fractions 01050 f1.setValue(static_cast<FIELD3D_VEC3_T<double> >(c2) - p); 01051 // C2 fraction 01052 f2.setValue(static_cast<FIELD3D_VEC3_T<double> >(1.0) - f1); 01053 01054 // Clamp the coordinates 01055 c1.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c1.x)); 01056 c1.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c1.y)); 01057 c1.z = std::min(dataWindow.max.z + 1, std::max(dataWindow.min.z, c1.z)); 01058 c2.x = std::min(dataWindow.max.x, std::max(dataWindow.min.x, c2.x)); 01059 c2.y = std::min(dataWindow.max.y, std::max(dataWindow.min.y, c2.y)); 01060 c2.z = std::min(dataWindow.max.z + 1, std::max(dataWindow.min.z, c2.z)); 01061 01062 ret = (f1.x * (f1.y * (f1.z * data.w(c1.x, c1.y, c1.z) + 01063 f2.z * data.w(c1.x, c1.y, c2.z)) + 01064 f2.y * (f1.z * data.w(c1.x, c2.y, c1.z) + 01065 f2.z * data.w(c1.x, c2.y, c2.z))) + 01066 f2.x * (f1.y * (f1.z * data.w(c2.x, c1.y, c1.z) + 01067 f2.z * data.w(c2.x, c1.y, c2.z)) + 01068 f2.y * (f1.z * data.w(c2.x, c2.y, c1.z) + 01069 f2.z * data.w(c2.x, c2.y, c2.z)))); 01070 break; 01071 } 01072 default: 01073 break; 01074 } 01075 01076 return ret; 01077 } 01078 01079 //----------------------------------------------------------------------------// 01080 01081 template <class Field_T> 01082 typename Field_T::value_type 01083 CubicGenericFieldInterp<Field_T>::sample(const Field_T &data, 01084 const V3d &vsP) const 01085 { 01086 typedef typename Field_T::value_type Data_T; 01087 01088 // Pixel centers are at .5 coordinates 01089 // NOTE: Don't use contToDisc for this, we're looking for sample 01090 // point locations, not coordinate shifts. 01091 V3d clampedVsP(std::max(0.5, vsP.x), 01092 std::max(0.5, vsP.y), 01093 std::max(0.5, vsP.z)); 01094 FIELD3D_VEC3_T<double> p(clampedVsP - FIELD3D_VEC3_T<double>(0.5)); 01095 01096 const Box3i &dataWindow = data.dataWindow(); 01097 01098 // Lower left corner 01099 V3i c(static_cast<int>(floor(p.x)), 01100 static_cast<int>(floor(p.y)), 01101 static_cast<int>(floor(p.z))); 01102 01103 // Fractions 01104 FIELD3D_VEC3_T<double> t(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 01105 01106 // Clamp the coordinates 01107 int im, jm, km; 01108 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 01109 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 01110 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 01111 int im_1, jm_1, km_1; 01112 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 01113 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 01114 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 01115 int im1, jm1, km1; 01116 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 01117 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 01118 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 01119 int im2, jm2, km2; 01120 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 01121 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 01122 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 01123 01124 Data_T z11 = monotonicCubicInterpolant(data.fastValue(im_1, jm_1, km_1), 01125 data.fastValue(im_1, jm_1, km), 01126 data.fastValue(im_1, jm_1, km1), 01127 data.fastValue(im_1, jm_1, km2), t.z); 01128 Data_T z12 = monotonicCubicInterpolant(data.fastValue(im_1, jm, km_1), 01129 data.fastValue(im_1, jm, km), 01130 data.fastValue(im_1, jm, km1), 01131 data.fastValue(im_1, jm, km2), t.z); 01132 Data_T z13 = monotonicCubicInterpolant(data.fastValue(im_1, jm1, km_1), 01133 data.fastValue(im_1, jm1, km), 01134 data.fastValue(im_1, jm1, km1), 01135 data.fastValue(im_1, jm1, km2), t.z); 01136 Data_T z14 = monotonicCubicInterpolant(data.fastValue(im_1, jm2, km_1), 01137 data.fastValue(im_1, jm2, km), 01138 data.fastValue(im_1, jm2, km1), 01139 data.fastValue(im_1, jm2, km2), t.z); 01140 01141 Data_T z21 = monotonicCubicInterpolant(data.fastValue(im, jm_1, km_1), 01142 data.fastValue(im, jm_1, km), 01143 data.fastValue(im, jm_1, km1), 01144 data.fastValue(im, jm_1, km2), t.z); 01145 Data_T z22 = monotonicCubicInterpolant(data.fastValue(im, jm, km_1), 01146 data.fastValue(im, jm, km), 01147 data.fastValue(im, jm, km1), 01148 data.fastValue(im, jm, km2), t.z); 01149 Data_T z23 = monotonicCubicInterpolant(data.fastValue(im, jm1, km_1), 01150 data.fastValue(im, jm1, km), 01151 data.fastValue(im, jm1, km1), 01152 data.fastValue(im, jm1, km2), t.z); 01153 Data_T z24 = monotonicCubicInterpolant(data.fastValue(im, jm2, km_1), 01154 data.fastValue(im, jm2, km), 01155 data.fastValue(im, jm2, km1), 01156 data.fastValue(im, jm2, km2), t.z); 01157 01158 Data_T z31 = monotonicCubicInterpolant(data.fastValue(im1, jm_1, km_1), 01159 data.fastValue(im1, jm_1, km), 01160 data.fastValue(im1, jm_1, km1), 01161 data.fastValue(im1, jm_1, km2), t.z); 01162 Data_T z32 = monotonicCubicInterpolant(data.fastValue(im1, jm, km_1), 01163 data.fastValue(im1, jm, km), 01164 data.fastValue(im1, jm, km1), 01165 data.fastValue(im1, jm, km2), t.z); 01166 Data_T z33 = monotonicCubicInterpolant(data.fastValue(im1, jm1, km_1), 01167 data.fastValue(im1, jm1, km), 01168 data.fastValue(im1, jm1, km1), 01169 data.fastValue(im1, jm1, km2), t.z); 01170 Data_T z34 = monotonicCubicInterpolant(data.fastValue(im1, jm2, km_1), 01171 data.fastValue(im1, jm2, km), 01172 data.fastValue(im1, jm2, km1), 01173 data.fastValue(im1, jm2, km2), t.z); 01174 01175 Data_T z41 = monotonicCubicInterpolant(data.fastValue(im2, jm_1, km_1), 01176 data.fastValue(im2, jm_1, km), 01177 data.fastValue(im2, jm_1, km1), 01178 data.fastValue(im2, jm_1, km2), t.z); 01179 Data_T z42 = monotonicCubicInterpolant(data.fastValue(im2, jm, km_1), 01180 data.fastValue(im2, jm, km), 01181 data.fastValue(im2, jm, km1), 01182 data.fastValue(im2, jm, km2), t.z); 01183 Data_T z43 = monotonicCubicInterpolant(data.fastValue(im2, jm1, km_1), 01184 data.fastValue(im2, jm1, km), 01185 data.fastValue(im2, jm1, km1), 01186 data.fastValue(im2, jm1, km2), t.z); 01187 Data_T z44 = monotonicCubicInterpolant(data.fastValue(im2, jm2, km_1), 01188 data.fastValue(im2, jm2, km), 01189 data.fastValue(im2, jm2, km1), 01190 data.fastValue(im2, jm2, km2), t.z); 01191 01192 Data_T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 01193 Data_T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 01194 Data_T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 01195 Data_T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 01196 01197 Data_T z0 = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 01198 01199 return z0; 01200 } 01201 01202 //----------------------------------------------------------------------------// 01203 01204 template <class Data_T> 01205 Data_T CubicMACFieldInterp<Data_T>::sample(const MACField<Data_T> &data, 01206 const V3d &vsP) const 01207 { 01208 typedef typename Data_T::BaseType T; 01209 01210 const Box3i &dataWindow = data.dataWindow(); 01211 01212 // Pixel centers are at .5 coordinates 01213 // NOTE: Don't use contToDisc for this, we're looking for sample 01214 // point locations, not coordinate shifts. 01215 01216 Data_T ret; 01217 01218 // X component --- 01219 01220 V3d clampedVsP(std::max(0.5, vsP.x), 01221 std::max(0.5, vsP.y), 01222 std::max(0.5, vsP.z)); 01223 FIELD3D_VEC3_T<double> p(vsP.x, 01224 clampedVsP.y - 0.5, 01225 clampedVsP.z - 0.5); 01226 01227 // Lower left corner 01228 V3i c(static_cast<int>(floor(p.x)), 01229 static_cast<int>(floor(p.y)), 01230 static_cast<int>(floor(p.z))); 01231 01232 FIELD3D_VEC3_T<double> t(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 01233 01234 { 01235 // Clamp the coordinates 01236 int im, jm, km; 01237 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x + 1)); 01238 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 01239 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 01240 int im_1, jm_1, km_1; 01241 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x + 1)); 01242 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 01243 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 01244 int im1, jm1, km1; 01245 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x + 1)); 01246 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 01247 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 01248 int im2, jm2, km2; 01249 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x + 1)); 01250 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 01251 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 01252 01253 T z11 = monotonicCubicInterpolant(data.u(im_1, jm_1, km_1), 01254 data.u(im_1, jm_1, km), 01255 data.u(im_1, jm_1, km1), 01256 data.u(im_1, jm_1, km2), t.z); 01257 T z12 = monotonicCubicInterpolant(data.u(im_1, jm, km_1), 01258 data.u(im_1, jm, km), 01259 data.u(im_1, jm, km1), 01260 data.u(im_1, jm, km2), t.z); 01261 T z13 = monotonicCubicInterpolant(data.u(im_1, jm1, km_1), 01262 data.u(im_1, jm1, km), 01263 data.u(im_1, jm1, km1), 01264 data.u(im_1, jm1, km2), t.z); 01265 T z14 = monotonicCubicInterpolant(data.u(im_1, jm2, km_1), 01266 data.u(im_1, jm2, km), 01267 data.u(im_1, jm2, km1), 01268 data.u(im_1, jm2, km2), t.z); 01269 01270 T z21 = monotonicCubicInterpolant(data.u(im, jm_1, km_1), 01271 data.u(im, jm_1, km), 01272 data.u(im, jm_1, km1), 01273 data.u(im, jm_1, km2), t.z); 01274 T z22 = monotonicCubicInterpolant(data.u(im, jm, km_1), 01275 data.u(im, jm, km), 01276 data.u(im, jm, km1), 01277 data.u(im, jm, km2), t.z); 01278 T z23 = monotonicCubicInterpolant(data.u(im, jm1, km_1), 01279 data.u(im, jm1, km), 01280 data.u(im, jm1, km1), 01281 data.u(im, jm1, km2), t.z); 01282 T z24 = monotonicCubicInterpolant(data.u(im, jm2, km_1), 01283 data.u(im, jm2, km), 01284 data.u(im, jm2, km1), 01285 data.u(im, jm2, km2), t.z); 01286 01287 T z31 = monotonicCubicInterpolant(data.u(im1, jm_1, km_1), 01288 data.u(im1, jm_1, km), 01289 data.u(im1, jm_1, km1), 01290 data.u(im1, jm_1, km2), t.z); 01291 T z32 = monotonicCubicInterpolant(data.u(im1, jm, km_1), 01292 data.u(im1, jm, km), 01293 data.u(im1, jm, km1), 01294 data.u(im1, jm, km2), t.z); 01295 T z33 = monotonicCubicInterpolant(data.u(im1, jm1, km_1), 01296 data.u(im1, jm1, km), 01297 data.u(im1, jm1, km1), 01298 data.u(im1, jm1, km2), t.z); 01299 T z34 = monotonicCubicInterpolant(data.u(im1, jm2, km_1), 01300 data.u(im1, jm2, km), 01301 data.u(im1, jm2, km1), 01302 data.u(im1, jm2, km2), t.z); 01303 01304 T z41 = monotonicCubicInterpolant(data.u(im2, jm_1, km_1), 01305 data.u(im2, jm_1, km), 01306 data.u(im2, jm_1, km1), 01307 data.u(im2, jm_1, km2), t.z); 01308 T z42 = monotonicCubicInterpolant(data.u(im2, jm, km_1), 01309 data.u(im2, jm, km), 01310 data.u(im2, jm, km1), 01311 data.u(im2, jm, km2), t.z); 01312 T z43 = monotonicCubicInterpolant(data.u(im2, jm1, km_1), 01313 data.u(im2, jm1, km), 01314 data.u(im2, jm1, km1), 01315 data.u(im2, jm1, km2), t.z); 01316 T z44 = monotonicCubicInterpolant(data.u(im2, jm2, km_1), 01317 data.u(im2, jm2, km), 01318 data.u(im2, jm2, km1), 01319 data.u(im2, jm2, km2), t.z); 01320 01321 T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 01322 T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 01323 T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 01324 T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 01325 01326 ret.x = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 01327 } 01328 01329 01330 // Y component --- 01331 01332 p.setValue(clampedVsP.x - 0.5, vsP.y , clampedVsP.z - 0.5); 01333 01334 // Lower left corner 01335 c.x = static_cast<int>(floor(p.x)); 01336 c.y = static_cast<int>(floor(p.y)); 01337 c.z = static_cast<int>(floor(p.z)); 01338 01339 t.setValue(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 01340 { 01341 // Clamp the coordinates 01342 int im, jm, km; 01343 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 01344 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y + 1)); 01345 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z)); 01346 int im_1, jm_1, km_1; 01347 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 01348 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y + 1)); 01349 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z)); 01350 int im1, jm1, km1; 01351 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 01352 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y + 1)); 01353 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z)); 01354 int im2, jm2, km2; 01355 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 01356 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y + 1)); 01357 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z)); 01358 01359 T z11 = monotonicCubicInterpolant(data.v(im_1, jm_1, km_1), 01360 data.v(im_1, jm_1, km), 01361 data.v(im_1, jm_1, km1), 01362 data.v(im_1, jm_1, km2), t.z); 01363 T z12 = monotonicCubicInterpolant(data.v(im_1, jm, km_1), 01364 data.v(im_1, jm, km), 01365 data.v(im_1, jm, km1), 01366 data.v(im_1, jm, km2), t.z); 01367 T z13 = monotonicCubicInterpolant(data.v(im_1, jm1, km_1), 01368 data.v(im_1, jm1, km), 01369 data.v(im_1, jm1, km1), 01370 data.v(im_1, jm1, km2), t.z); 01371 T z14 = monotonicCubicInterpolant(data.v(im_1, jm2, km_1), 01372 data.v(im_1, jm2, km), 01373 data.v(im_1, jm2, km1), 01374 data.v(im_1, jm2, km2), t.z); 01375 01376 T z21 = monotonicCubicInterpolant(data.v(im, jm_1, km_1), 01377 data.v(im, jm_1, km), 01378 data.v(im, jm_1, km1), 01379 data.v(im, jm_1, km2), t.z); 01380 T z22 = monotonicCubicInterpolant(data.v(im, jm, km_1), 01381 data.v(im, jm, km), 01382 data.v(im, jm, km1), 01383 data.v(im, jm, km2), t.z); 01384 T z23 = monotonicCubicInterpolant(data.v(im, jm1, km_1), 01385 data.v(im, jm1, km), 01386 data.v(im, jm1, km1), 01387 data.v(im, jm1, km2), t.z); 01388 T z24 = monotonicCubicInterpolant(data.v(im, jm2, km_1), 01389 data.v(im, jm2, km), 01390 data.v(im, jm2, km1), 01391 data.v(im, jm2, km2), t.z); 01392 01393 T z31 = monotonicCubicInterpolant(data.v(im1, jm_1, km_1), 01394 data.v(im1, jm_1, km), 01395 data.v(im1, jm_1, km1), 01396 data.v(im1, jm_1, km2), t.z); 01397 T z32 = monotonicCubicInterpolant(data.v(im1, jm, km_1), 01398 data.v(im1, jm, km), 01399 data.v(im1, jm, km1), 01400 data.v(im1, jm, km2), t.z); 01401 T z33 = monotonicCubicInterpolant(data.v(im1, jm1, km_1), 01402 data.v(im1, jm1, km), 01403 data.v(im1, jm1, km1), 01404 data.v(im1, jm1, km2), t.z); 01405 T z34 = monotonicCubicInterpolant(data.v(im1, jm2, km_1), 01406 data.v(im1, jm2, km), 01407 data.v(im1, jm2, km1), 01408 data.v(im1, jm2, km2), t.z); 01409 01410 T z41 = monotonicCubicInterpolant(data.v(im2, jm_1, km_1), 01411 data.v(im2, jm_1, km), 01412 data.v(im2, jm_1, km1), 01413 data.v(im2, jm_1, km2), t.z); 01414 T z42 = monotonicCubicInterpolant(data.v(im2, jm, km_1), 01415 data.v(im2, jm, km), 01416 data.v(im2, jm, km1), 01417 data.v(im2, jm, km2), t.z); 01418 T z43 = monotonicCubicInterpolant(data.v(im2, jm1, km_1), 01419 data.v(im2, jm1, km), 01420 data.v(im2, jm1, km1), 01421 data.v(im2, jm1, km2), t.z); 01422 T z44 = monotonicCubicInterpolant(data.v(im2, jm2, km_1), 01423 data.v(im2, jm2, km), 01424 data.v(im2, jm2, km1), 01425 data.v(im2, jm2, km2), t.z); 01426 01427 T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 01428 T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 01429 T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 01430 T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 01431 01432 ret.y = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 01433 } 01434 01435 // Z component --- 01436 01437 p.setValue(clampedVsP.x - 0.5 , clampedVsP.y - 0.5, vsP.z); 01438 01439 // Lower left corner 01440 c.x = static_cast<int>(floor(p.x)); 01441 c.y = static_cast<int>(floor(p.y)); 01442 c.z = static_cast<int>(floor(p.z)); 01443 01444 t.setValue(p - static_cast<FIELD3D_VEC3_T<double> >(c)); 01445 { 01446 // Clamp the coordinates 01447 int im, jm, km; 01448 im = std::max(dataWindow.min.x, std::min(c.x, dataWindow.max.x)); 01449 jm = std::max(dataWindow.min.y, std::min(c.y, dataWindow.max.y)); 01450 km = std::max(dataWindow.min.z, std::min(c.z, dataWindow.max.z + 1)); 01451 int im_1, jm_1, km_1; 01452 im_1 = std::max(dataWindow.min.x, std::min(im - 1, dataWindow.max.x)); 01453 jm_1 = std::max(dataWindow.min.y, std::min(jm - 1, dataWindow.max.y)); 01454 km_1 = std::max(dataWindow.min.z, std::min(km - 1, dataWindow.max.z + 1)); 01455 int im1, jm1, km1; 01456 im1 = std::max(dataWindow.min.x, std::min(im + 1, dataWindow.max.x)); 01457 jm1 = std::max(dataWindow.min.y, std::min(jm + 1, dataWindow.max.y)); 01458 km1 = std::max(dataWindow.min.z, std::min(km + 1, dataWindow.max.z + 1)); 01459 int im2, jm2, km2; 01460 im2 = std::max(dataWindow.min.x, std::min(im + 2, dataWindow.max.x)); 01461 jm2 = std::max(dataWindow.min.y, std::min(jm + 2, dataWindow.max.y)); 01462 km2 = std::max(dataWindow.min.z, std::min(km + 2, dataWindow.max.z + 1)); 01463 01464 T z11 = monotonicCubicInterpolant(data.w(im_1, jm_1, km_1), 01465 data.w(im_1, jm_1, km), 01466 data.w(im_1, jm_1, km1), 01467 data.w(im_1, jm_1, km2), t.z); 01468 T z12 = monotonicCubicInterpolant(data.w(im_1, jm, km_1), 01469 data.w(im_1, jm, km), 01470 data.w(im_1, jm, km1), 01471 data.w(im_1, jm, km2), t.z); 01472 T z13 = monotonicCubicInterpolant(data.w(im_1, jm1, km_1), 01473 data.w(im_1, jm1, km), 01474 data.w(im_1, jm1, km1), 01475 data.w(im_1, jm1, km2), t.z); 01476 T z14 = monotonicCubicInterpolant(data.w(im_1, jm2, km_1), 01477 data.w(im_1, jm2, km), 01478 data.w(im_1, jm2, km1), 01479 data.w(im_1, jm2, km2), t.z); 01480 01481 T z21 = monotonicCubicInterpolant(data.w(im, jm_1, km_1), 01482 data.w(im, jm_1, km), 01483 data.w(im, jm_1, km1), 01484 data.w(im, jm_1, km2), t.z); 01485 T z22 = monotonicCubicInterpolant(data.w(im, jm, km_1), 01486 data.w(im, jm, km), 01487 data.w(im, jm, km1), 01488 data.w(im, jm, km2), t.z); 01489 T z23 = monotonicCubicInterpolant(data.w(im, jm1, km_1), 01490 data.w(im, jm1, km), 01491 data.w(im, jm1, km1), 01492 data.w(im, jm1, km2), t.z); 01493 T z24 = monotonicCubicInterpolant(data.w(im, jm2, km_1), 01494 data.w(im, jm2, km), 01495 data.w(im, jm2, km1), 01496 data.w(im, jm2, km2), t.z); 01497 01498 T z31 = monotonicCubicInterpolant(data.w(im1, jm_1, km_1), 01499 data.w(im1, jm_1, km), 01500 data.w(im1, jm_1, km1), 01501 data.w(im1, jm_1, km2), t.z); 01502 T z32 = monotonicCubicInterpolant(data.w(im1, jm, km_1), 01503 data.w(im1, jm, km), 01504 data.w(im1, jm, km1), 01505 data.w(im1, jm, km2), t.z); 01506 T z33 = monotonicCubicInterpolant(data.w(im1, jm1, km_1), 01507 data.w(im1, jm1, km), 01508 data.w(im1, jm1, km1), 01509 data.w(im1, jm1, km2), t.z); 01510 T z34 = monotonicCubicInterpolant(data.w(im1, jm2, km_1), 01511 data.w(im1, jm2, km), 01512 data.w(im1, jm2, km1), 01513 data.w(im1, jm2, km2), t.z); 01514 01515 T z41 = monotonicCubicInterpolant(data.w(im2, jm_1, km_1), 01516 data.w(im2, jm_1, km), 01517 data.w(im2, jm_1, km1), 01518 data.w(im2, jm_1, km2), t.z); 01519 T z42 = monotonicCubicInterpolant(data.w(im2, jm, km_1), 01520 data.w(im2, jm, km), 01521 data.w(im2, jm, km1), 01522 data.w(im2, jm, km2), t.z); 01523 T z43 = monotonicCubicInterpolant(data.w(im2, jm1, km_1), 01524 data.w(im2, jm1, km), 01525 data.w(im2, jm1, km1), 01526 data.w(im2, jm1, km2), t.z); 01527 T z44 = monotonicCubicInterpolant(data.w(im2, jm2, km_1), 01528 data.w(im2, jm2, km), 01529 data.w(im2, jm2, km1), 01530 data.w(im2, jm2, km2), t.z); 01531 01532 T y1 = monotonicCubicInterpolant(z11, z12, z13, z14, t.y); 01533 T y2 = monotonicCubicInterpolant(z21, z22, z23, z24, t.y); 01534 T y3 = monotonicCubicInterpolant(z31, z32, z33, z34, t.y); 01535 T y4 = monotonicCubicInterpolant(z41, z42, z43, z44, t.y); 01536 01537 ret.z = monotonicCubicInterpolant(y1, y2, y3, y4, t.x); 01538 } 01539 01540 return ret; 01541 } 01542 01543 //----------------------------------------------------------------------------// 01544 01545 template <class Data_T> 01546 Data_T 01547 ProceduralFieldLookup<Data_T>::sample(const ProceduralField<Data_T> &data, 01548 const V3d &vsP) const 01549 { 01550 V3d voxelScale = V3d(1.0) / data.dataResolution(); 01551 V3d lsP = vsP * voxelScale; 01552 return data.lsSample(lsP); 01553 } 01554 01555 //----------------------------------------------------------------------------// 01556 01557 template <class S, class T> 01558 FIELD3D_VEC3_T<T> operator * (S s, const FIELD3D_VEC3_T<T> vec) 01559 { 01560 return FIELD3D_VEC3_T<T>(vec.x * s, vec.y * s, vec.z * s); 01561 } 01562 01563 //----------------------------------------------------------------------------// 01564 01565 template<class T> 01566 T monotonicCubicInterpolant(const T &f1, const T &f2, const T &f3, const T &f4, 01567 double t) 01568 { 01569 T d_k = T(.5) * (f3 - f1); 01570 T d_k1 = T(.5) * (f4 - f2); 01571 T delta_k = f3 - f2; 01572 01573 if (delta_k == static_cast<T>(0)) { 01574 d_k = static_cast<T>(0); 01575 d_k1 = static_cast<T>(0); 01576 } 01577 01578 T a0 = f2; 01579 T a1 = d_k; 01580 T a2 = (T(3) * delta_k) - (T(2) * d_k) - d_k1; 01581 T a3 = d_k + d_k1 - (T(2) * delta_k); 01582 01583 T t1 = t; 01584 T t2 = t1 * t1; 01585 T t3 = t2 * t1; 01586 01587 return a3 * t3 + a2 * t2 + a1 * t1 + a0; 01588 } 01589 01590 //----------------------------------------------------------------------------// 01591 01593 // References: 01594 // http://en.wikipedia.org/wiki/Monotone_cubic_interpolation 01595 // http://en.wikipedia.org/wiki/Cubic_Hermite_spline 01596 template <class Data_T> 01597 Data_T monotonicCubicInterpolantVec(const Data_T &f1, const Data_T &f2, 01598 const Data_T &f3, const Data_T &f4, 01599 double t) 01600 { 01601 typedef typename Data_T::BaseType T; 01602 01603 Data_T d_k = T(.5) * (f3 - f1); 01604 Data_T d_k1 = T(.5) * (f4 - f2); 01605 Data_T delta_k = f3 - f2; 01606 01607 for (int i = 0; i < 3; i++) { 01608 if (delta_k[i] == static_cast<T>(0)) { 01609 d_k[i] = static_cast<T>(0); 01610 d_k1[i]= static_cast<T>(0); 01611 } 01612 } 01613 01614 Data_T a0 = f2; 01615 Data_T a1 = d_k; 01616 Data_T a2 = (delta_k * T(3)) - (d_k * T(2)) - d_k1; 01617 Data_T a3 = d_k + d_k1 - (delta_k * T(2)); 01618 01619 T t1 = t; 01620 T t2 = t1 * t1; 01621 T t3 = t2 * t1; 01622 01623 return a3 * t3 + a2 * t2 + a1 * t1 + a0; 01624 } 01625 01626 //----------------------------------------------------------------------------// 01627 // Template specializations 01628 //----------------------------------------------------------------------------// 01629 01630 template<> 01631 inline 01632 V3h monotonicCubicInterpolant<V3h>(const V3h &f1, const V3h &f2, 01633 const V3h &f3, const V3h &f4, double t) 01634 { 01635 return monotonicCubicInterpolantVec(f1, f2, f3, f4, t); 01636 } 01637 01638 //----------------------------------------------------------------------------// 01639 01640 template<> 01641 inline 01642 V3f monotonicCubicInterpolant<V3f>(const V3f &f1, const V3f &f2, 01643 const V3f &f3, const V3f &f4, double t) 01644 { 01645 return monotonicCubicInterpolantVec(f1, f2, f3, f4, t); 01646 } 01647 01648 //----------------------------------------------------------------------------// 01649 01650 template<> 01651 inline 01652 V3d monotonicCubicInterpolant<V3d>(const V3d &f1, const V3d &f2, 01653 const V3d &f3, const V3d &f4, double t) 01654 { 01655 return monotonicCubicInterpolantVec(f1, f2, f3, f4, t); 01656 } 01657 01658 //----------------------------------------------------------------------------// 01659 01660 FIELD3D_NAMESPACE_HEADER_CLOSE 01661 01662 //----------------------------------------------------------------------------// 01663 01664 #endif // Include guard