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