Field3D
SparseField.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 
00042 //----------------------------------------------------------------------------//
00043 
00044 #ifndef _INCLUDED_Field3D_SparseField_H_
00045 #define _INCLUDED_Field3D_SparseField_H_
00046 
00047 //----------------------------------------------------------------------------//
00048 
00049 #include <vector>
00050 
00051 #include <boost/lexical_cast.hpp>
00052 
00053 #include "Field.h"
00054 #include "SparseFile.h"
00055 
00056 #define BLOCK_ORDER 4 // 2^BLOCK_ORDER is the block size along each axis
00057 
00058 //----------------------------------------------------------------------------//
00059 
00060 #include "ns.h"
00061 
00062 FIELD3D_NAMESPACE_OPEN
00063 
00064 //----------------------------------------------------------------------------//
00065 // Forward declarations 
00066 //----------------------------------------------------------------------------//
00067 
00068 template <class Field_T>
00069 class LinearGenericFieldInterp;
00070 template <class Field_T>
00071 class CubicGenericFieldInterp; 
00072 
00073 //----------------------------------------------------------------------------//
00074 // SparseBlock
00075 //----------------------------------------------------------------------------//
00076 
00079 namespace Sparse {
00080 
00084 template <typename Data_T>
00085 struct SparseBlock 
00086 {
00087   // Constructors --------------------------------------------------------------
00088 
00090   SparseBlock()
00091     : isAllocated(false),
00092       emptyValue(static_cast<Data_T>(0))
00093   { /* Empty */ }
00094 
00095   // Main methods --------------------------------------------------------------
00096 
00098   Data_T& value(int i, int j, int k, int blockOrder)
00100   { return data[(k << blockOrder << blockOrder) + (j << blockOrder) + i]; }
00101 
00104   const Data_T& value(int i, int j, int k, int blockOrder) const
00105   { return data[(k << blockOrder << blockOrder) + (j << blockOrder) + i]; }
00106 
00108   void resize(int n)
00109   { this->data.resize(n); }
00110 
00112   void clear()
00113   { std::vector<Data_T>().swap(data); }
00114 
00116   Data_T& dataRef()
00117   { return data[0]; }
00118 
00119   // Data members --------------------------------------------------------------
00120 
00122   bool isAllocated;
00123 
00127   Data_T emptyValue;
00128 
00130   std::vector<Data_T> data;
00131 
00132 private:
00133 
00136 
00137 };
00138 
00139 } // namespace Sparse
00140 
00141 //----------------------------------------------------------------------------//
00142 // SparseField
00143 //----------------------------------------------------------------------------//
00144 
00157 //----------------------------------------------------------------------------//
00158 
00159 template <class Data_T>
00160 class SparseField
00161   : public ResizableField<Data_T>
00162 {
00163 public:
00164 
00165   // Typedefs ------------------------------------------------------------------
00166   
00167   typedef boost::intrusive_ptr<SparseField> Ptr;
00168   typedef std::vector<Ptr> Vec;
00169 
00170   typedef LinearGenericFieldInterp<SparseField<Data_T> > LinearInterp;
00171   typedef CubicGenericFieldInterp<SparseField<Data_T> > CubicInterp;
00172 
00173   // RTTI replacement ----------------------------------------------------------
00174 
00175   typedef SparseField<Data_T> class_type;
00176   DEFINE_FIELD_RTTI_CONCRETE_CLASS;
00177 
00178   static const char *staticClassName()
00179   {
00180     return "SparseField";
00181   }
00182   
00183   static const char *classType()
00184   {
00185     return SparseField<Data_T>::ms_classType.name();
00186   }
00187 
00188   // Constructors --------------------------------------------------------------
00189 
00192 
00194   SparseField();
00195 
00197   SparseField(const SparseField &o);
00198 
00200   ~SparseField();
00201 
00204   SparseField& operator=(const SparseField &o);
00205   
00206   // \}
00207 
00208   // Main methods --------------------------------------------------------------
00209 
00211   virtual void clear(const Data_T &value);
00212 
00215   void setBlockOrder(int order);
00216 
00218   int blockOrder() const;
00219 
00221   int blockSize() const;
00222 
00224   bool voxelIsInAllocatedBlock(int i, int j, int k) const;
00225 
00227   bool blockIsAllocated(int bi, int bj, int bk) const;
00228   
00231   const Data_T getBlockEmptyValue(int bi, int bj, int bk) const;
00232 
00235   void setBlockEmptyValue(int bi, int bj, int bk, const Data_T &val);
00236 
00238   bool blockIndexIsValid(int bi, int bj, int bk) const;
00239 
00241   V3i blockRes() const;
00242 
00248   template <typename Functor_T>
00249   int releaseBlocks(Functor_T func);
00250 
00252   int blockId(int blockI, int blockJ, int blockK) const;
00253 
00257   void getBlockCoord(int i, int j, int k, int &bi, int &bj, int &bk) const;
00258   
00262   void getVoxelInBlock(int i, int j, int k, int &vi, int &vj, int &vk) const;
00263 
00265   void applyDataWindowOffset(int &i, int &j, int &k) const
00266   {
00267     i -= base::m_dataWindow.min.x;
00268     j -= base::m_dataWindow.min.y;
00269     k -= base::m_dataWindow.min.z;
00270   }
00271   
00272   // From Field base class -----------------------------------------------------
00273 
00276   virtual Data_T value(int i, int j, int k) const;
00277   virtual long long int memSize() const;
00279   
00280   // From WritableField base class ---------------------------------------------
00281 
00284   virtual Data_T& lvalue(int i, int j, int k);
00286 
00287   // Concrete voxel access -----------------------------------------------------
00288 
00290   Data_T fastValue(int i, int j, int k) const;
00292   Data_T& fastLValue(int i, int j, int k);
00293 
00294   // From FieldBase ------------------------------------------------------------
00295 
00298   virtual std::string className() const
00299   { return staticClassName(); }
00300 
00301   virtual FieldBase::Ptr clone() const
00302   { return Ptr(new SparseField(*this)); }
00303 
00305   
00306   // Iterators -----------------------------------------------------------------
00307 
00310 
00312   class const_iterator;
00313 
00315   const_iterator cbegin() const;
00317   const_iterator cbegin(const Box3i &subset) const;
00319   const_iterator cend() const;
00322   const_iterator cend(const Box3i &subset) const;
00323 
00327   class iterator;
00328 
00330   iterator begin();
00332   iterator begin(const Box3i &subset);
00334   iterator end();
00337   iterator end(const Box3i &subset);
00338 
00342   class block_iterator;
00343 
00344   block_iterator blockBegin() const;
00346   block_iterator blockEnd() const;
00347 
00349 
00350   // Internal utility functions ------------------------------------------------
00351 
00354   void addReference(const std::string &filename, const std::string &layerPath,
00355                     int valuesPerBlock, int occupiedBlocks);
00358   void setupReferenceBlocks();
00359   
00360  protected:
00361 
00362   friend class SparseFieldIO;
00363 
00364   // Typedefs ------------------------------------------------------------------
00365 
00366   typedef ResizableField<Data_T> base;
00367   typedef Sparse::SparseBlock<Data_T> Block;
00368 
00369   // From ResizableField class -------------------------------------------------
00370 
00371   virtual void sizeChanged()
00372   { 
00373     // Call base class
00374     base::sizeChanged();
00375     setupBlocks(); 
00376   }
00377 
00378   // Convenience methods -------------------------------------------------------
00379 
00382 
00384   void setupBlocks();
00385 
00387   void deallocBlock(Block &block, const Data_T &emptyValue);
00388 
00390 
00391   // Data members --------------------------------------------------------------
00392 
00394   int m_blockOrder;
00396   V3i m_blockRes;
00398   int m_blockXYSize;
00400   std::vector<Block> m_blocks;
00401 
00404   SparseFileManager *m_fileManager;
00406   int m_fileId;
00407 
00409   Data_T m_dummy;
00410 
00411 private:
00412 
00413   // Static data members -------------------------------------------------------
00414 
00415   static TemplatedFieldType<SparseField<Data_T> > ms_classType;
00416 
00417   // Utility methods -----------------------------------------------------------
00418 
00421   void copySparseField(const SparseField &o);
00422 
00425   void copyBlockStates(const SparseField<Data_T> &o);
00426 
00427 };
00428 
00429 //----------------------------------------------------------------------------//
00430 // Static member instantiations
00431 //----------------------------------------------------------------------------//
00432 
00433 FIELD3D_CLASSTYPE_TEMPL_INSTANTIATION(SparseField);
00434 
00435 //----------------------------------------------------------------------------//
00436 // Typedefs
00437 //----------------------------------------------------------------------------//
00438 
00439 typedef SparseField<half>   SparseFieldh;
00440 typedef SparseField<float>  SparseFieldf;
00441 typedef SparseField<double> SparseFieldd;
00442 typedef SparseField<V3h>    SparseField3h;
00443 typedef SparseField<V3f>    SparseField3f;
00444 typedef SparseField<V3d>    SparseField3d;
00445 
00446 //----------------------------------------------------------------------------//
00447 // Helper functors
00448 //----------------------------------------------------------------------------//
00449 
00450 namespace Sparse {
00451 
00454 template <typename Data_T>
00455 struct CheckAllEqual
00456 {
00465   bool check(const SparseBlock<Data_T> &block, Data_T &retEmptyValue,
00466              const V3i &validSize, const V3i &blockSize)
00467   {
00468     // Store first value
00469     Data_T first = block.data[0];
00470     // Iterate over rest
00471     bool match = true;
00472     if (validSize == blockSize) {
00473       // interior block so look at all voxels
00474       for (typename std::vector<Data_T>::const_iterator i = block.data.begin();
00475            i != block.data.end(); ++i) {
00476         if (*i != first) {
00477           match = false;
00478           break;
00479         }
00480       }
00481     } else {
00482       // only look at valid voxels
00483       int x=0, y=0, z=0;
00484       for (typename std::vector<Data_T>::const_iterator i = block.data.begin();
00485            i != block.data.end(); ++i, ++x) {
00486         if (x >= blockSize.x) {
00487           x = 0;
00488           ++y;
00489           if (y >= blockSize.y) {
00490             y = 0;
00491             ++z;
00492           }
00493         }
00494         if (x >= validSize.x || y >= validSize.y || z >= validSize.z) {
00495           continue;
00496         }
00497 
00498         if (*i != first) {
00499           match = false;
00500           break;
00501         }
00502       }
00503     } // end of interior block test
00504 
00505     if (match) {
00506       retEmptyValue = first;
00507       return true;
00508     } else {
00509       return false;
00510     }
00511   }
00512 };
00513 
00514 //----------------------------------------------------------------------------//
00515 
00516 template <typename Data_T>
00517 inline bool isAnyLess(const Data_T &left, const Data_T &right)
00518 {
00519   return (std::abs(left) < right);
00520 }
00521 
00522 //----------------------------------------------------------------------------//
00523 
00524 template <>
00525 inline bool isAnyLess(const V3h &left, const V3h &right)
00526 {
00527   return (std::abs(left.x) < right.x || 
00528           std::abs(left.y) < right.y || 
00529           std::abs(left.z) < right.z );
00530 }
00531 
00532 //----------------------------------------------------------------------------//
00533 
00534 template <>
00535 inline bool isAnyLess(const V3f &left, const V3f &right)
00536 {
00537   return (std::abs(left.x) < right.x || 
00538           std::abs(left.y) < right.y || 
00539           std::abs(left.z) < right.z );
00540 }
00541 
00542 //----------------------------------------------------------------------------//
00543 
00544 template <>
00545 inline bool isAnyLess(const V3d &left, const V3d &right)
00546 {
00547   return (std::abs(left.x) < right.x || 
00548           std::abs(left.y) < right.y || 
00549           std::abs(left.z) < right.z );
00550 }
00551 
00552 //----------------------------------------------------------------------------//
00553 
00557 template <typename Data_T>
00558 struct CheckMaxAbs
00559 {
00561   CheckMaxAbs(Data_T maxValue)
00562     : m_maxValue(maxValue)
00563   { }
00572   bool check(const SparseBlock<Data_T> &block, Data_T &retEmptyValue,
00573              const V3i &validSize, const V3i &blockSize)
00574   {
00575     // Store first value
00576     Data_T first = block.data[0];
00577     // Iterate over rest
00578     bool allGreater = true;
00579     if (validSize == blockSize) {
00580       // interior block so look at all voxels
00581       for (typename std::vector<Data_T>::const_iterator i = block.data.begin();
00582            i != block.data.end(); ++i) {
00583         if (isAnyLess<Data_T>(*i, m_maxValue)) {
00584           allGreater = false;
00585           break;
00586         }
00587       }
00588     } else {
00589       // only look at valid voxels
00590       int x=0, y=0, z=0;
00591       for (typename std::vector<Data_T>::const_iterator i = block.data.begin();
00592            i != block.data.end(); ++i, ++x) {
00593         if (x >= blockSize.x) {
00594           x = 0;
00595           ++y;
00596           if (y >= blockSize.y) {
00597             y = 0;
00598             ++z;
00599           }
00600         }
00601         if (x >= validSize.x || y >= validSize.y || z >= validSize.z) {
00602           continue;
00603         }
00604         if (isAnyLess<Data_T>(*i, m_maxValue)) {
00605           allGreater = false;
00606           break;
00607         }
00608       }
00609     } // end of interior block test
00610 
00611     if (allGreater) {
00612       retEmptyValue = first;
00613       return true;
00614     } else {
00615       return false;
00616     }
00617   }
00618 private:
00619   Data_T m_maxValue;
00620 };
00621 
00622 } // namespace Sparse
00623 
00624 //----------------------------------------------------------------------------//
00625 // SparseField::const_iterator
00626 //----------------------------------------------------------------------------//
00627 
00629 template <class Data_T>
00630 class SparseField<Data_T>::const_iterator
00631 {
00632  public:
00633   typedef SparseField<Data_T> class_type;
00634   const_iterator(const class_type &field,
00635                  const Box3i &window,
00636                  const V3i &currentPos, int blockOrder)
00637     : x(currentPos.x), y(currentPos.y), z(currentPos.z), 
00638       m_p(NULL), m_blockIsActivated(false), 
00639       m_blockStepsTicker(0), m_blockOrder(blockOrder), 
00640       m_blockId(-1), m_window(window), m_field(&field)
00641   {
00642     m_manager = m_field->m_fileManager;
00643     setupNextBlock(x, y, z);
00644   }
00645   ~const_iterator() {
00646     if (m_manager && m_blockId >= 0 && 
00647         m_blockId < static_cast<int>(m_field->m_blocks.size())) {
00648       if (m_field->m_blocks[m_blockId].isAllocated)
00649         m_manager->decBlockRef<Data_T>(m_field->m_fileId, m_blockId);
00650     }
00651   }
00652   const const_iterator& operator ++ ()
00653   {
00654     bool resetPtr = false;
00655     // Check against end of data window
00656     if (x == m_window.max.x) {
00657       if (y == m_window.max.y) {
00658         x = m_window.min.x;
00659         y = m_window.min.y;
00660         ++z;
00661         resetPtr = true;
00662       } else {
00663         x = m_window.min.x;
00664         ++y;
00665         resetPtr = true;
00666       }
00667     } else {
00668       ++x;
00669     }
00670     // These can both safely be incremented here
00671     ++m_blockStepsTicker;
00672     // ... but only step forward if we're in a non-empty block
00673     if (!m_isEmptyBlock && (!m_manager || m_blockIsActivated))
00674       ++m_p;   
00675     // Check if we've reached the end of this block
00676     if (m_blockStepsTicker == (1 << m_blockOrder))
00677       resetPtr = true;
00678     if (resetPtr) {
00679       // If we have, we need to reset the current block, etc.
00680       m_blockStepsTicker = 0;
00681       setupNextBlock(x, y, z);
00682     }
00683     return *this;
00684   }
00685   template <class Iter_T>
00686   inline bool operator == (const Iter_T &rhs) const
00687   {
00688     return x == rhs.x && y == rhs.y && z == rhs.z;
00689   }
00690   template <class Iter_T>
00691   inline bool operator != (const Iter_T &rhs) const
00692   {
00693     return x != rhs.x || y != rhs.y || z != rhs.z;
00694   }
00695   inline const Data_T& operator * () const
00696   {
00697     if (!m_isEmptyBlock && m_manager && !m_blockIsActivated) {
00698       m_manager->activateBlock<Data_T>(m_field->m_fileId, m_blockId);
00699       m_blockIsActivated = true;
00700       const Block &block = m_field->m_blocks[m_blockId];
00701       int vi, vj, vk;
00702       m_field->getVoxelInBlock(x, y, z, vi, vj, vk);
00703       m_p = &block.value(vi, vj, vk, m_blockOrder);
00704     }
00705     return *m_p;
00706   }
00707   inline const Data_T* operator -> () const
00708   {
00709     if (!m_isEmptyBlock && m_manager && !m_blockIsActivated) {
00710       SparseFileManager *manager = m_field->m_fileManager;
00711       manager->activateBlock<Data_T>(m_field->m_fileId, m_blockId);
00712       m_blockIsActivated = true;
00713       const Block &block = m_field->m_blocks[m_blockId];
00714       int vi, vj, vk;
00715       m_field->getVoxelInBlock(x, y, z, vi, vj, vk);
00716       m_p = &block.value(vi, vj, vk, m_blockOrder);
00717     }
00718     return m_p;
00719   }
00720 
00721   // Public data members -------------------------------------------------------
00722 
00724   int x, y, z;
00725 
00726 private:
00727 
00728   // Typedefs ------------------------------------------------------------------
00729 
00730   typedef Sparse::SparseBlock<Data_T> Block;
00731 
00732   // Convenience methods -------------------------------------------------------
00733 
00734   void setupNextBlock(int i, int j, int k) 
00735   {
00736     m_field->applyDataWindowOffset(i, j, k);
00737     m_field->getBlockCoord(i, j, k, m_blockI, m_blockJ, m_blockK);
00738     int oldBlockId = m_blockId;
00739     m_blockId = m_field->blockId(m_blockI, m_blockJ, m_blockK);
00740     if (m_manager && oldBlockId != m_blockId &&
00741         oldBlockId >= 0 && 
00742         oldBlockId < static_cast<int>(m_field->m_blocks.size()) &&
00743         m_field->m_blocks[oldBlockId].isAllocated) {
00744       m_manager->decBlockRef<Data_T>(m_field->m_fileId, oldBlockId);
00745     }
00746     if (m_blockId >= m_field->m_blockXYSize * m_field->m_blockRes.z) {
00747       m_isEmptyBlock = true;
00748       return;
00749     }
00750 
00751     const Block &block = m_field->m_blocks[m_blockId];
00752     int vi, vj, vk;
00753     m_field->getVoxelInBlock(i, j, k, vi, vj, vk);      
00754     m_blockStepsTicker = vi;
00755     if (block.isAllocated) {
00756       if (m_manager && oldBlockId != m_blockId && m_blockId >= 0) {
00757         m_manager->incBlockRef<Data_T>(m_field->m_fileId, m_blockId);
00758         // this is a managed field, so the block may not be loaded
00759         // yet, so don't bother setting m_p yet (it'll get set in the
00760         // * and -> operators when the block is activated)
00761       } else {
00762         // only set m_p to the voxel's address if this is not a
00763         // managed field, i.e., if the data is already in memory.
00764         m_p = &block.value(vi, vj, vk, m_blockOrder);
00765       }
00766       m_isEmptyBlock = false;
00767     } else {
00768       m_p = &block.emptyValue;
00769       m_isEmptyBlock = true;
00770     }
00771     if (m_field->m_fileManager) {
00772       m_blockIsActivated = false;
00773     }
00774   }
00775 
00777   mutable const Data_T *m_p;
00779   bool m_isEmptyBlock;
00782   mutable bool m_blockIsActivated;
00784   int m_blockStepsTicker;
00786   int m_blockOrder;
00788   int m_blockI, m_blockJ, m_blockK, m_blockId;
00790   Box3i m_window;
00792   const class_type *m_field;
00794   SparseFileManager *m_manager;
00795 };
00796 
00797 //----------------------------------------------------------------------------//
00798 // SparseField::iterator
00799 //----------------------------------------------------------------------------/
00800 
00802 template <class Data_T>
00803 class SparseField<Data_T>::iterator
00804 {
00805  public:
00806   typedef SparseField<Data_T> class_type;
00807   iterator(class_type &field,
00808            const Box3i &window,
00809            const V3i &currentPos, int blockOrder)
00810     : x(currentPos.x), y(currentPos.y), z(currentPos.z),
00811       m_p(NULL), m_blockStepsTicker(0), m_blockOrder(blockOrder), 
00812       m_blockId(-1), m_window(window), m_field(&field)
00813   {
00814     setupNextBlock(x, y, z);
00815   }
00816   const iterator& operator ++ ()
00817   {
00818     bool resetPtr = false;
00819     // Check against end of data window
00820     if (x == m_window.max.x) {
00821       if (y == m_window.max.y) {
00822         x = m_window.min.x;
00823         y = m_window.min.y;
00824         ++z;
00825         resetPtr = true;
00826       } else {
00827         x = m_window.min.x;
00828         ++y;
00829         resetPtr = true;
00830       }
00831     } else {
00832       ++x;
00833     }
00834     // These can both safely be incremented here
00835     ++m_blockStepsTicker;
00836     // ... but only step forward if we're in a non-empty block
00837     if (!m_isEmptyBlock)
00838       ++m_p;   
00839     // Check if we've reached the end of this block
00840     if (m_blockStepsTicker == (1 << m_blockOrder))
00841       resetPtr = true;
00842     if (resetPtr) {
00843       // If we have, we need to reset the current block, etc.
00844       m_blockStepsTicker = 0;
00845       setupNextBlock(x, y, z);
00846     }
00847     return *this;
00848   }
00849   inline bool operator == (const iterator &rhs) const
00850   {
00851     return x == rhs.x && y == rhs.y && z == rhs.z;
00852   }
00853   inline bool operator != (const iterator &rhs) const
00854   {
00855     return x != rhs.x || y != rhs.y || z != rhs.z;
00856   }
00857   inline Data_T& operator * ()
00858   {
00859     if (m_field->m_fileManager) {
00860       assert(false && "Dereferencing iterator on a dynamic-read sparse field");
00861       Msg::print(Msg::SevWarning, "Dereferencing iterator on a dynamic-read "
00862                 "sparse field");
00863       return *m_p;      
00864     }
00865     // If the block is currently empty, we must allocate it
00866     if (m_isEmptyBlock) {
00867       // Touch the voxel to allocate the block
00868       m_field->lvalue(x, y, z);
00869       // Set up the block again
00870       setupNextBlock(x, y, z);
00871     }
00872     return *m_p;
00873   }
00874   inline Data_T* operator -> ()
00875   {
00876     if (m_field->m_fileManager) {
00877       assert(false && "Dereferencing iterator on a dynamic-read sparse field");
00878       Msg::print(Msg::SevWarning, "Dereferencing iterator on a dynamic-read "
00879                 "sparse field");
00880       return m_p;      
00881     }
00882     // If the block is currently empty, we must allocate it
00883     if (m_isEmptyBlock) {
00884       // Touch the voxel to allocate the block
00885       m_field->lvalue(x, y, z);
00886       // Set up the block again
00887       setupNextBlock(x, y, z);
00888     }
00889     return m_p;
00890   }
00891   // Public data members
00892   int x, y, z;
00893 private:
00894   typedef Sparse::SparseBlock<Data_T> Block;
00896   void setupNextBlock(int i, int j, int k) 
00897   {
00898     m_field->applyDataWindowOffset(i, j, k);
00899     m_field->getBlockCoord(i, j, k, m_blockI, m_blockJ, m_blockK);
00900     m_blockId = m_field->blockId(m_blockI, m_blockJ, m_blockK);
00901     if (m_blockId >= m_field->m_blockXYSize * m_field->m_blockRes.z) {
00902       m_isEmptyBlock = true;
00903       return;
00904     }
00905     Block &block = m_field->m_blocks[m_blockId];
00906     int vi, vj, vk;
00907     m_field->getVoxelInBlock(i, j, k, vi, vj, vk);      
00908     m_blockStepsTicker = vi;
00909     if (block.isAllocated) {
00910       m_p = &block.value(vi, vj, vk, m_blockOrder);
00911       m_isEmptyBlock = false;
00912     } else {
00913       m_p = &block.emptyValue;
00914       m_isEmptyBlock = true;
00915     }
00916   }
00918   Data_T *m_p;
00920   bool m_isEmptyBlock;
00922   int m_blockStepsTicker;
00924   int m_blockOrder;
00926   int m_blockI, m_blockJ, m_blockK, m_blockId;
00928   Box3i m_window;
00930   class_type *m_field;
00931 };
00932 
00933 //----------------------------------------------------------------------------//
00934 // SparseField::block_iterator
00935 //----------------------------------------------------------------------------/
00936 
00939 template <class Data_T>
00940 class SparseField<Data_T>::block_iterator
00941 {
00942  public:
00944   typedef SparseField<Data_T> class_type;
00946   block_iterator(const class_type &field, const Box3i &window,
00947                  const V3i &currentPos)
00948     : x(currentPos.x), y(currentPos.y), z(currentPos.z), 
00949       m_window(window), m_field(field)
00950   {
00951     recomputeBlockBoundingBox();
00952   }
00954   const block_iterator& operator ++ ()
00955   {
00956     if (x == m_window.max.x) {
00957       if (y == m_window.max.y) {
00958         x = m_window.min.x;
00959         y = m_window.min.y;
00960         ++z;
00961       } else {
00962         x = m_window.min.x; 
00963         ++y;
00964       }
00965     } else {
00966       ++x;
00967     }
00968     recomputeBlockBoundingBox();
00969     return *this;
00970   }
00972   inline bool operator == (const block_iterator &rhs) const
00973   {
00974     return x == rhs.x && y == rhs.y && z == rhs.z;
00975   }
00977   inline bool operator != (const block_iterator &rhs) const
00978   {
00979     return x != rhs.x || y != rhs.y || z != rhs.z;
00980   }
00982   const Box3i& blockBoundingBox()
00983   {
00984     return m_currentBlockWindow;
00985   }
00987   int x, y, z;
00988 private:
00989   void recomputeBlockBoundingBox()
00990   {
00991     Box3i box;
00992     int blockSize = m_field.blockSize();
00993     box.min = V3i(x * blockSize, y * blockSize, z * blockSize);
00994     box.max = box.min + V3i(blockSize - 1, blockSize - 1, blockSize - 1);
00995     // Clamp the box 
00996     box.min = FIELD3D_CLIP(box.min, m_field.dataWindow());
00997     box.max = FIELD3D_CLIP(box.max, m_field.dataWindow());
00998     // Set the member variable
00999     m_currentBlockWindow = box;
01000   }
01002   Box3i m_window;
01004   const class_type& m_field;
01006   Box3i m_currentBlockWindow;
01007 };
01008 
01009 //----------------------------------------------------------------------------//
01010 // SparseField implementations
01011 //----------------------------------------------------------------------------//
01012 
01013 template <class Data_T>
01014 SparseField<Data_T>::SparseField()
01015   : base(),
01016     m_blockOrder(BLOCK_ORDER),
01017     m_fileManager(NULL)
01018 { 
01019   setupBlocks();
01020 }
01021 
01022 //----------------------------------------------------------------------------//
01023 
01024 template <class Data_T>
01025 SparseField<Data_T>::SparseField(const SparseField<Data_T> &o) :
01026   base(o)
01027 {
01028   copySparseField(o);
01029 }
01030 
01031 //----------------------------------------------------------------------------//
01032 
01033 template <class Data_T>
01034 SparseField<Data_T>::~SparseField()
01035 {
01036   if (m_fileManager) {
01037     // this file is dynamically managed, so we need to ensure the
01038     // cache doesn't point to this field's blocks because they are
01039     // about to be deleted
01040     m_fileManager->removeFieldFromCache<Data_T>(m_fileId);
01041   }
01042 }
01043 
01044 //----------------------------------------------------------------------------//
01045 
01046 template <class Data_T>
01047 SparseField<Data_T> &
01048 SparseField<Data_T>::operator=(const SparseField<Data_T> &o)
01049 {
01050   if (this != &o) {
01051     this->base::operator=(o);
01052     copySparseField(o);
01053   }
01054   return *this;
01055 }
01056 
01057 //----------------------------------------------------------------------------//
01058 
01059 template <class Data_T>
01060 void
01061 SparseField<Data_T>::copySparseField(const SparseField<Data_T> &o)
01062 {
01063   m_blockOrder = o.m_blockOrder;
01064   if (o.m_fileManager) {
01065     // allocate m_blocks, sets m_blockRes, m_blockXYSize, m_blocks
01066     setupBlocks();
01067 
01068     m_fileManager = o.m_fileManager;
01069     SparseFile::Reference<Data_T> &oldReference = 
01070       m_fileManager->reference<Data_T>(o.m_fileId);
01071     addReference(oldReference.filename, oldReference.layerPath,
01072                  oldReference.valuesPerBlock,
01073                  oldReference.occupiedBlocks);
01074     copyBlockStates(o);
01075     setupReferenceBlocks();
01076   } else {
01077     // directly copy all values and blocks from the source, no extra setup
01078     m_blockRes = o.m_blockRes;
01079     m_blockXYSize = o.m_blockXYSize;
01080     m_blocks = o.m_blocks;
01081     m_fileId = -1;
01082     m_fileManager = NULL;
01083   }
01084 }
01085 
01086 //----------------------------------------------------------------------------//
01087 
01088 template <class Data_T>
01089 void SparseField<Data_T>::addReference(const std::string &filename,
01090                                        const std::string &layerPath,
01091                                        int valuesPerBlock,
01092                                        int occupiedBlocks)
01093 {
01094   m_fileManager = &SparseFileManager::singleton();
01095   m_fileId = m_fileManager->getNextId<Data_T>(filename, layerPath);
01096   // Set up the manager data
01097   SparseFile::Reference<Data_T> &reference = 
01098     m_fileManager->reference<Data_T>(m_fileId);
01099   reference.valuesPerBlock = valuesPerBlock;
01100   reference.occupiedBlocks = occupiedBlocks;
01101   reference.setNumBlocks(m_blocks.size());
01102 }
01103 
01104 //----------------------------------------------------------------------------//
01105 
01106 template <class Data_T>
01107 void SparseField<Data_T>::copyBlockStates(const SparseField<Data_T> &o)
01108 {
01109   if (m_blocks.size() != o.m_blocks.size()) return;
01110 
01111   typename std::vector<Sparse::SparseBlock<Data_T> >::iterator b =
01112     m_blocks.begin();
01113   typename std::vector<Sparse::SparseBlock<Data_T> >::iterator bend = 
01114     m_blocks.end();
01115   typename std::vector<Sparse::SparseBlock<Data_T> >::const_iterator ob =
01116     o.m_blocks.begin();
01117 
01118   for (; b != bend; ++b, ++ob) {
01119     b->isAllocated = ob->isAllocated;
01120     b->emptyValue = ob->emptyValue;
01121     b->clear();
01122   }
01123 }
01124 
01125 //----------------------------------------------------------------------------//
01126 
01127 template <class Data_T>
01128 void SparseField<Data_T>::setupReferenceBlocks()
01129 {
01130   if (!m_fileManager || m_fileId < 0) return;
01131 
01132   SparseFile::Reference<Data_T> &reference = 
01133     m_fileManager->reference<Data_T>(m_fileId);
01134 
01135   std::vector<int>::iterator fb = reference.fileBlockIndices.begin();
01136   typename SparseFile::Reference<Data_T>::BlockPtrs::iterator bp = 
01137     reference.blocks.begin();
01138   typename std::vector<Sparse::SparseBlock<Data_T> >::iterator b =
01139     m_blocks.begin();
01140   typename std::vector<Sparse::SparseBlock<Data_T> >::iterator bend =
01141     m_blocks.end();
01142   int nextBlockIdx = 0;
01143 
01144   for (; b != bend; ++b, ++fb, ++bp) {
01145     if (b->isAllocated) {
01146       *fb = nextBlockIdx;
01147       *bp = &(*b);
01148       nextBlockIdx++;
01149     } else {
01150       *fb = -1;
01151     }
01152   }
01153 }
01154 
01155 //----------------------------------------------------------------------------//
01156 
01157 template <class Data_T>
01158 void SparseField<Data_T>::clear(const Data_T &value)
01159 {
01160   // If we're clearing, we can get rid of all current blocks
01161   setupBlocks();
01162   // Then just fill in the default values
01163   typename std::vector<Block>::iterator i;
01164   typename std::vector<Block>::iterator end;
01165   for (i = m_blocks.begin(), end = m_blocks.end(); i != end; ++i) {
01166     i->emptyValue = value;
01167   }
01168 }
01169 
01170 //----------------------------------------------------------------------------//
01171 
01172 template <class Data_T>
01173 void SparseField<Data_T>::setBlockOrder(int order)
01174 {
01175   m_blockOrder = order;
01176   setupBlocks();
01177 }
01178 
01179 //----------------------------------------------------------------------------//
01180 
01181 template <class Data_T>
01182 int SparseField<Data_T>::blockOrder() const 
01183 {
01184   return m_blockOrder;
01185 }
01186 
01187 //----------------------------------------------------------------------------//
01188 
01189 template <class Data_T>
01190 int SparseField<Data_T>::blockSize() const
01191 { 
01192   return 1 << m_blockOrder;
01193 }
01194 
01195 //----------------------------------------------------------------------------//
01196 
01197 template <class Data_T>
01198 bool SparseField<Data_T>::voxelIsInAllocatedBlock(int i, int j, int k) const
01199 {
01200   int bi, bj, bk;
01201   applyDataWindowOffset(i, j, k);
01202   getBlockCoord(i, j, k, bi, bj, bk);
01203   return blockIsAllocated(bi, bj, bk);
01204 }
01205 
01206 //----------------------------------------------------------------------------//
01207 
01208 template <class Data_T>
01209 bool SparseField<Data_T>::blockIsAllocated(int bi, int bj, int bk) const
01210 {
01211   const Block &block = m_blocks[blockId(bi, bj, bk)];
01212   return block.isAllocated;
01213 }
01214   
01215 //----------------------------------------------------------------------------//
01216 
01217 template <class Data_T>
01218 const Data_T SparseField<Data_T>::getBlockEmptyValue(int bi, int bj, int bk) const
01219 {
01220   return m_blocks[blockId(bi, bj, bk)].emptyValue;
01221 }
01222 
01223 //----------------------------------------------------------------------------//
01224 
01225 template <class Data_T>
01226 void SparseField<Data_T>::setBlockEmptyValue(int bi, int bj, int bk, 
01227                                              const Data_T &val)
01228 {
01229   Block &block = m_blocks[blockId(bi, bj, bk)];
01230   if (block.isAllocated) {
01231     deallocBlock(block, val);
01232   } else {
01233     block.emptyValue = val;
01234   }
01235 }
01236 
01237 //----------------------------------------------------------------------------//
01238 
01239 template <class Data_T>
01240 bool SparseField<Data_T>::blockIndexIsValid(int bi, int bj, int bk) const
01241 {
01242   return bi >= 0 && bj >= 0 && bk >= 0 && 
01243     bi < m_blockRes.x && bj < m_blockRes.y && bk < m_blockRes.z;
01244 }
01245 
01246 //----------------------------------------------------------------------------//
01247 
01248 template <class Data_T>
01249 V3i SparseField<Data_T>::blockRes() const
01250 {
01251   return m_blockRes;
01252 }
01253 
01254 //----------------------------------------------------------------------------//
01255 
01256 template <class Data_T>
01257 template <typename Functor_T>
01258 int SparseField<Data_T>::releaseBlocks(Functor_T func)
01259 {
01260   Data_T emptyValue;
01261   int numDeallocs = 0;
01262   typename std::vector<Block>::iterator i;
01263 
01264   // If the block is on the edge of the field, it may have unused
01265   // voxels, with undefined values.  We need to pass the range of
01266   // valid voxels into the check function, so it only looks at valid
01267   // voxels.
01268   V3i dataRes = FieldRes::dataResolution();
01269   V3i validSize;
01270   V3i blockAllocSize(blockSize());
01271   int bx, by, bz;
01272 
01273   for (i = m_blocks.begin(), bx=0, by=0, bz=0; i != m_blocks.end(); ++i, ++bx) {
01274     if (bx >= m_blockRes.x) {
01275       bx = 0;
01276       ++by;
01277       if (by >= m_blockRes.y) {
01278         by = 0;
01279         ++bz;
01280       }
01281     }
01282     validSize = blockAllocSize;
01283     if (bx == m_blockRes.x-1) {
01284       validSize.x = dataRes.x - bx * blockAllocSize.x;
01285     }
01286     if (by == m_blockRes.y-1) {
01287       validSize.y = dataRes.y - by * blockAllocSize.y;
01288     }
01289     if (bz == m_blockRes.z-1) {
01290       validSize.z = dataRes.z - bz * blockAllocSize.z;
01291     }
01292 
01293     if (i->isAllocated) {
01294       if (func.check(*i, emptyValue, validSize, blockAllocSize)) {
01295         deallocBlock(*i, emptyValue);
01296         numDeallocs++;
01297       }
01298     }
01299   }
01300   return numDeallocs;
01301 }
01302 
01303 //----------------------------------------------------------------------------//
01304 
01305 template <class Data_T>
01306 Data_T SparseField<Data_T>::value(int i, int j, int k) const
01307 {
01308   return fastValue(i, j, k);
01309 }
01310 
01311 //----------------------------------------------------------------------------//
01312 
01313 template <class Data_T>
01314 Data_T& SparseField<Data_T>::lvalue(int i, int j, int k)
01315 {
01316   return fastLValue(i, j, k);
01317 }
01318 
01319 //----------------------------------------------------------------------------//
01320 
01321 template <class Data_T>
01322 Data_T SparseField<Data_T>::fastValue(int i, int j, int k) const
01323 {
01324   assert (i >= base::m_dataWindow.min.x);
01325   assert (i <= base::m_dataWindow.max.x);
01326   assert (j >= base::m_dataWindow.min.y);
01327   assert (j <= base::m_dataWindow.max.y);
01328   assert (k >= base::m_dataWindow.min.z);
01329   assert (k <= base::m_dataWindow.max.z);
01330   // Add crop window offset
01331   applyDataWindowOffset(i, j, k);
01332   // Find block coord
01333   int bi, bj, bk;
01334   getBlockCoord(i, j, k, bi, bj, bk);
01335   // Find coord in block
01336   int vi, vj, vk;
01337   getVoxelInBlock(i, j, k, vi, vj, vk);
01338   // Get the actual block
01339   int id = blockId(bi, bj, bk);
01340 
01341   const Block &block = m_blocks[id];
01342   // Check if block data is allocated
01343   if (block.isAllocated) {
01344     if (m_fileManager) {
01345       m_fileManager->incBlockRef<Data_T>(m_fileId, id);
01346       m_fileManager->activateBlock<Data_T>(m_fileId, id);
01347       Data_T tmpValue = block.value(vi, vj, vk, m_blockOrder);
01348       m_fileManager->decBlockRef<Data_T>(m_fileId, id);
01349       return tmpValue;
01350     } else {
01351       return block.value(vi, vj, vk, m_blockOrder);
01352     }
01353   } else {
01354     return block.emptyValue;
01355   }
01356 }
01357 
01358 //----------------------------------------------------------------------------//
01359 
01361 template <class Data_T>
01362 Data_T& SparseField<Data_T>::fastLValue(int i, int j, int k)
01363 {
01364   assert (i >= base::m_dataWindow.min.x);
01365   assert (i <= base::m_dataWindow.max.x);
01366   assert (j >= base::m_dataWindow.min.y);
01367   assert (j <= base::m_dataWindow.max.y);
01368   assert (k >= base::m_dataWindow.min.z);
01369   assert (k <= base::m_dataWindow.max.z);
01370 
01371   if (m_fileManager) {
01372     assert(false && "Called fastLValue() on a dynamic-read sparse field");
01373     Msg::print(Msg::SevWarning, "Called fastLValue() on a dynamic-read "
01374               "sparse field");
01375     return m_dummy;
01376   }
01377 
01378   // Add crop window offset
01379   applyDataWindowOffset(i, j, k);
01380   // Find block coord
01381   int bi, bj, bk;
01382   getBlockCoord(i, j, k, bi, bj, bk);
01383   // Find coord in block
01384   int vi, vj, vk;
01385   getVoxelInBlock(i, j, k, vi, vj, vk);
01386   // Get the actual block
01387   int id = blockId(bi, bj, bk);
01388   Block &block = m_blocks[id];
01389   // If block is allocated, return a reference to the data
01390   if (block.isAllocated) {
01391     return block.value(vi, vj, vk, m_blockOrder);
01392   } else {
01393     // ... Otherwise, allocate block
01394     block.isAllocated = true;
01395     block.data.resize(1 << m_blockOrder << m_blockOrder << m_blockOrder);
01396     std::fill(block.data.begin(), block.data.end(), block.emptyValue);
01397     return block.value(vi, vj, vk, m_blockOrder);
01398   }
01399 }
01400 
01401 //----------------------------------------------------------------------------//
01402 
01403 template <class Data_T>
01404 long long int SparseField<Data_T>::memSize() const
01405 {
01406   long long int blockSize = m_blocks.capacity() * sizeof(Block);
01407   long long int dataSize = 0;
01408   typename std::vector<Block>::const_iterator i;
01409   for (i = m_blocks.begin(); i != m_blocks.end(); ++i) {
01410     if (i->isAllocated) {
01411       dataSize += i->data.capacity() * sizeof(Data_T);
01412     }
01413   }
01414   return sizeof(*this) + dataSize + blockSize;
01415 }
01416 
01417 //----------------------------------------------------------------------------//
01418 
01419 template <class Data_T>
01420 typename SparseField<Data_T>::const_iterator 
01421 SparseField<Data_T>::cbegin() const
01422 { 
01423   if (FieldRes::dataResolution() == V3i(0))
01424     return cend();
01425   return const_iterator(*this, base::m_dataWindow, base::m_dataWindow.min,
01426                         m_blockOrder); 
01427 }
01428 
01429 //----------------------------------------------------------------------------//
01430 
01431 template <class Data_T>
01432 typename SparseField<Data_T>::const_iterator 
01433 SparseField<Data_T>::cbegin(const Box3i &subset) const
01434 { 
01435   if (subset.isEmpty())
01436     return cend(subset);
01437   return const_iterator(*this, subset, subset.min, m_blockOrder); 
01438 }
01439 
01440 //----------------------------------------------------------------------------//
01441 
01442 template <class Data_T>
01443 typename SparseField<Data_T>::const_iterator 
01444 SparseField<Data_T>::cend() const
01445 { 
01446   return const_iterator(*this, base::m_dataWindow, 
01447                         V3i(base::m_dataWindow.min.x, 
01448                                    base::m_dataWindow.min.y,
01449                                    base::m_dataWindow.max.z + 1),
01450                         m_blockOrder);
01451 }
01452 
01453 //----------------------------------------------------------------------------//
01454 
01455 template <class Data_T>
01456 typename SparseField<Data_T>::const_iterator 
01457 SparseField<Data_T>::cend(const Box3i &subset) const
01458 { 
01459   return const_iterator(*this, subset, 
01460                         V3i(subset.min.x, 
01461                                    subset.min.y,
01462                                    subset.max.z + 1), m_blockOrder);
01463 }
01464 
01465 //----------------------------------------------------------------------------//
01466 
01467 template <class Data_T>
01468 typename SparseField<Data_T>::iterator 
01469 SparseField<Data_T>::begin()
01470 { 
01471   if (FieldRes::dataResolution() == V3i(0))
01472     return end();
01473   return iterator(*this, base::m_dataWindow, 
01474                   base::m_dataWindow.min, m_blockOrder); }
01475 
01476 //----------------------------------------------------------------------------//
01477 
01478 template <class Data_T>
01479 typename SparseField<Data_T>::iterator 
01480 SparseField<Data_T>::begin(const Box3i &subset)
01481 { 
01482   if (subset.isEmpty())
01483     return end(subset);
01484   return iterator(*this, subset, subset.min, m_blockOrder); 
01485 }
01486 
01487 //----------------------------------------------------------------------------//
01488 
01489 template <class Data_T>
01490 typename SparseField<Data_T>::iterator 
01491 SparseField<Data_T>::end()
01492 { 
01493   return iterator(*this, base::m_dataWindow, 
01494                   V3i(base::m_dataWindow.min.x, 
01495                              base::m_dataWindow.min.y,
01496                              base::m_dataWindow.max.z + 1), m_blockOrder);
01497 }
01498 
01499 //----------------------------------------------------------------------------//
01500 
01501 template <class Data_T>
01502 typename SparseField<Data_T>::iterator 
01503 SparseField<Data_T>::end(const Box3i &subset)
01504 { 
01505   return iterator(*this, subset, 
01506                   V3i(subset.min.x, subset.min.y, subset.max.z + 1), 
01507                   m_blockOrder);
01508 }
01509 
01510 //----------------------------------------------------------------------------//
01511 
01512 template <class Data_T>
01513 typename SparseField<Data_T>::block_iterator 
01514 SparseField<Data_T>::blockBegin() const
01515 { 
01516   if (FieldRes::dataResolution() == V3i(0))
01517     return blockEnd();
01518   return block_iterator(*this, Box3i(V3i(0), m_blockRes - V3i(1)), 
01519                         V3i(0)); 
01520 }
01521 
01522 //----------------------------------------------------------------------------//
01523 
01524 template <class Data_T>
01525 typename SparseField<Data_T>::block_iterator 
01526 SparseField<Data_T>::blockEnd() const
01527 { 
01528   return block_iterator(*this, Box3i(V3i(0), m_blockRes - V3i(1)), 
01529                         V3i(0, 0, m_blockRes.z));
01530 }
01531 
01532 //----------------------------------------------------------------------------//
01533 
01534 template <class Data_T>
01535 void SparseField<Data_T>::setupBlocks()
01536 {
01537   // Do calculation in floating point so we can round up later
01538   V3f res(base::m_dataWindow.size() + V3i(1));
01539   V3f blockRes(res / (1 << m_blockOrder));
01540   blockRes.x = ceil(blockRes.x);
01541   blockRes.y = ceil(blockRes.y);
01542   blockRes.z = ceil(blockRes.z);
01543   V3i intBlockRes(static_cast<int>(blockRes.x),
01544                          static_cast<int>(blockRes.y),
01545                          static_cast<int>(blockRes.z));
01546   m_blockRes = intBlockRes;
01547   m_blockXYSize = m_blockRes.x * m_blockRes.y;
01548   // clear() won't deallocate data. Do the swap trick.
01549   //m_blocks.clear();
01550   std::vector<Block>().swap(m_blocks);
01551   
01552   m_blocks.resize(intBlockRes.x * intBlockRes.y * intBlockRes.z);
01553 }
01554 
01555 //----------------------------------------------------------------------------//
01556 
01557 template <class Data_T>
01558 int SparseField<Data_T>::blockId(int blockI, int blockJ, int blockK) const
01559 {
01560   return blockK * m_blockXYSize + blockJ * m_blockRes.x + blockI;
01561 }
01562 
01563 //----------------------------------------------------------------------------//
01564 
01566 template <class Data_T>
01567 void SparseField<Data_T>::getBlockCoord(int i, int j, int k, 
01568                                         int &bi, int &bj, int &bk) const 
01569 {
01570   assert(i >= 0);
01571   assert(j >= 0);
01572   assert(k >= 0);
01573   bi = i >> m_blockOrder;
01574   bj = j >> m_blockOrder;
01575   bk = k >> m_blockOrder;
01576 }
01577 
01578 //----------------------------------------------------------------------------//
01579 
01581 template <class Data_T>
01582 void SparseField<Data_T>::getVoxelInBlock(int i, int j, int k, 
01583                                           int &vi, int &vj, int &vk) const
01584 {
01585   assert(i >= 0);
01586   assert(j >= 0);
01587   assert(k >= 0);
01588   vi = i & ((1 << m_blockOrder) - 1);
01589   vj = j & ((1 << m_blockOrder) - 1);
01590   vk = k & ((1 << m_blockOrder) - 1);
01591 }
01592 
01593 //----------------------------------------------------------------------------//
01594 
01595 template <class Data_T>
01596 void SparseField<Data_T>::deallocBlock(Block &block, const Data_T &emptyValue)
01597 {
01598   block.isAllocated = false;
01600   block.clear();
01601   block.emptyValue = emptyValue;
01602 }
01603 
01604 //----------------------------------------------------------------------------//
01605 
01606 FIELD3D_NAMESPACE_HEADER_CLOSE
01607 
01608 //----------------------------------------------------------------------------//
01609 
01610 #endif // Include guard