Field3D
FieldInterp.h
Go to the documentation of this file.
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