Field3D
SparseFieldIO.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_SparseFieldIO_H_
00047 #define _INCLUDED_Field3D_SparseFieldIO_H_
00048 
00049 //----------------------------------------------------------------------------//
00050 
00051 #include <string>
00052 #include <cmath>
00053 
00054 #include <hdf5.h>
00055 
00056 #include "SparseDataReader.h"
00057 #include "SparseField.h"
00058 #include "SparseFile.h"
00059 #include "FieldIO.h"
00060 #include "Field3DFile.h"
00061 
00062 //----------------------------------------------------------------------------//
00063 
00064 #include "ns.h"
00065 
00066 FIELD3D_NAMESPACE_OPEN
00067 
00068 //----------------------------------------------------------------------------//
00069 // SparseFieldIO
00070 //----------------------------------------------------------------------------//
00071 
00077 //----------------------------------------------------------------------------//
00078 
00079 class SparseFieldIO : public FieldIO 
00080 {
00081 
00082 public:
00083 
00084   // Typedefs ------------------------------------------------------------------
00085   
00086   typedef boost::intrusive_ptr<SparseFieldIO> Ptr;
00087 
00088   // Constructors --------------------------------------------------------------
00089 
00091   SparseFieldIO() 
00092    : FieldIO()
00093   { }
00094 
00096   virtual ~SparseFieldIO() 
00097   { /* Empty */ }
00098 
00099 
00100   static FieldIO::Ptr create()
00101   { return Ptr(new SparseFieldIO); }
00102 
00103   // From FieldIO --------------------------------------------------------------
00104 
00108   virtual FieldBase::Ptr read(hid_t layerGroup, const std::string &filename, 
00109                               const std::string &layerPath,
00110                               DataTypeEnum typeEnum);
00111 
00114   virtual bool write(hid_t layerGroup, FieldBase::Ptr field);
00115 
00117   virtual std::string className() const
00118   { return std::string("SparseField"); }
00119 
00120 private:
00121 
00122   // Internal methods ----------------------------------------------------------
00123 
00125   template <class Data_T>
00126   bool writeInternal(hid_t layerGroup, typename SparseField<Data_T>::Ptr field);
00127 
00129   template <class Data_T>
00130   bool readData(hid_t location, 
00131                 int numBlocks, 
00132                 const std::string &filename, 
00133                 const std::string &layerPath, 
00134                 typename SparseField<Data_T>::Ptr result);
00135 
00136   // Strings -------------------------------------------------------------------
00137 
00138   static const int         k_versionNumber;
00139   static const std::string k_versionAttrName;
00140   static const std::string k_extentsStr;
00141   static const std::string k_dataWindowStr;
00142   static const std::string k_componentsStr;
00143   static const std::string k_blockOrderStr;
00144   static const std::string k_numBlocksStr;
00145   static const std::string k_blockResStr;
00146   static const std::string k_bitsPerComponentStr;
00147   static const std::string k_numOccupiedBlocksStr;
00148   static const std::string k_dataStr;
00149 
00150 };
00151 
00152 //----------------------------------------------------------------------------//
00153 // Template methods
00154 //----------------------------------------------------------------------------//
00155 
00157 template <class Data_T>
00158 bool SparseFieldIO::writeInternal(hid_t layerGroup, 
00159                                   typename SparseField<Data_T>::Ptr field)
00160 {
00161   using namespace std;
00162   using namespace Exc;
00163   using namespace Hdf5Util;
00164   using namespace Sparse;
00165 
00166   Box3i ext(field->extents()), dw(field->dataWindow());
00167 
00168   int components = FieldTraits<Data_T>::dataDims();
00169 
00170   int valuesPerBlock = (1 << (field->m_blockOrder * 3)) * components;
00171 
00172   int size[3];
00173   size[0] = dw.max.x - dw.min.x + 1;
00174   size[1] = dw.max.y - dw.min.y + 1;
00175   size[2] = dw.max.z - dw.min.z + 1;
00176 
00177 
00178   hsize_t totalSize[1];
00179   totalSize[0] = size[0] * size[1] * size[2] * components;
00180 
00181   // Add extents attribute ---
00182 
00183   int extents[6] = 
00184     { ext.min.x, ext.min.y, ext.min.z, ext.max.x, ext.max.y, ext.max.z };
00185 
00186   if (!writeAttribute(layerGroup, k_extentsStr, 6, extents[0])) {
00187     Msg::print(Msg::SevWarning, "Error adding size attribute.");
00188     return false;
00189   }
00190 
00191   // Add data window attribute ---
00192 
00193   int dataWindow[6] = 
00194     { dw.min.x, dw.min.y, dw.min.z, dw.max.x, dw.max.y, dw.max.z };
00195 
00196   if (!writeAttribute(layerGroup, k_dataWindowStr, 6, dataWindow[0])) {
00197     Msg::print(Msg::SevWarning, "Error adding size attribute.");
00198     return false;
00199   }
00200 
00201   // Add components attribute ---
00202 
00203   if (!writeAttribute(layerGroup, k_componentsStr, 1, components)) {
00204     Msg::print(Msg::SevWarning, "Error adding components attribute.");
00205     return false;
00206   }
00207 
00208   // Add block order attribute ---
00209 
00210   int blockOrder = field->m_blockOrder;
00211 
00212   if (!writeAttribute(layerGroup, k_blockOrderStr, 1, blockOrder)) {
00213     Msg::print(Msg::SevWarning, "Error adding block order attribute.");
00214     return false;
00215   }
00216 
00217   // Add number of blocks attribute ---
00218   
00219   V3i &blockRes = field->m_blockRes;
00220   int numBlocks = blockRes.x * blockRes.y * blockRes.z;
00221 
00222   if (!writeAttribute(layerGroup, k_numBlocksStr, 1, numBlocks)) {
00223     Msg::print(Msg::SevWarning, "Error adding number of blocks attribute.");
00224     return false;
00225   }
00226 
00227   // Add block resolution in each dimension ---
00228 
00229   if (!writeAttribute(layerGroup, k_blockResStr, 3, blockRes.x)) {
00230     Msg::print(Msg::SevWarning, "Error adding block res attribute.");
00231     return false;
00232   }
00233 
00234   // Add the bits per component attribute ---
00235 
00236   int bits = DataTypeTraits<Data_T>::h5bits();
00237   if (!writeAttribute(layerGroup, k_bitsPerComponentStr, 1, bits)) {
00238     Msg::print(Msg::SevWarning, "Error adding bits per component attribute.");
00239     return false;    
00240   }
00241 
00242   // Write the block info data sets ---
00243   
00244   // ... Write the isAllocated array
00245   {
00246     vector<char> isAllocated(numBlocks);
00247     vector<char>::iterator i = isAllocated.begin();
00248     typename vector<SparseBlock<Data_T> >::const_iterator b = 
00249       field->m_blocks.begin();
00250     for (; i != isAllocated.end(); ++i, ++b)
00251       *i = static_cast<char>(b->isAllocated);
00252     writeSimpleData<char>(layerGroup, "block_is_allocated_data", isAllocated);
00253   }
00254 
00255   // ... Write the emptyValue array
00256   {
00257     vector<Data_T> emptyValue(numBlocks);
00258     typename vector<Data_T>::iterator i = emptyValue.begin();
00259     typename vector<SparseBlock<Data_T> >::const_iterator b = 
00260       field->m_blocks.begin();
00261     for (; i != emptyValue.end(); ++i, ++b)
00262       *i = static_cast<Data_T>(b->emptyValue);
00263     writeSimpleData<Data_T>(layerGroup, "block_empty_value_data", emptyValue);
00264   }
00265 
00266   // Count the number of occupied blocks ---
00267   int occupiedBlocks = 0;
00268   typename vector<SparseBlock<Data_T> >::iterator b = 
00269     field->m_blocks.begin();
00270   for (; b != field->m_blocks.end(); ++b) {
00271     if (b->isAllocated)
00272       occupiedBlocks++;
00273   }
00274   if (!writeAttribute(layerGroup, k_numOccupiedBlocksStr, 1, occupiedBlocks)) {
00275     throw WriteAttributeException("Couldn't add attribute " + 
00276                                 k_numOccupiedBlocksStr);
00277   }
00278   
00279   if (occupiedBlocks > 0) {
00280 
00281     // Make the memory data space
00282     hsize_t memDims[1];
00283     memDims[0] = valuesPerBlock;
00284     H5ScopedScreate memDataSpace(H5S_SIMPLE);
00285     H5Sset_extent_simple(memDataSpace.id(), 1, memDims, NULL);
00286 
00287     // Make the file data space
00288     hsize_t fileDims[2];
00289     fileDims[0] = occupiedBlocks;
00290     fileDims[1] = valuesPerBlock;
00291     H5ScopedScreate fileDataSpace(H5S_SIMPLE);
00292     H5Sset_extent_simple(fileDataSpace.id(), 2, fileDims, NULL);
00293 
00294     // Set up gzip property list
00295     bool gzipAvailable = checkHdf5Gzip();
00296     hid_t dcpl = H5Pcreate(H5P_DATASET_CREATE);
00297     hsize_t chunkSize[2];
00298     chunkSize[0] = 1;
00299     chunkSize[1] = valuesPerBlock;
00300     if (gzipAvailable) {
00301       herr_t status = H5Pset_deflate(dcpl, 9);
00302       if (status < 0) {
00303         return false;
00304       }
00305       status = H5Pset_chunk(dcpl, 2, chunkSize);
00306       if (status < 0) {
00307         return false;
00308       }    
00309     }
00310 
00311     // Add the data set
00312     H5ScopedDcreate dataSet(layerGroup, k_dataStr, 
00313                             DataTypeTraits<Data_T>::h5type(), 
00314                             fileDataSpace.id(), 
00315                             H5P_DEFAULT, dcpl, H5P_DEFAULT);
00316     if (dataSet.id() < 0)
00317       throw CreateDataSetException("Couldn't create data set in "
00318                                    "SparseFieldIO::writeInternal");
00319 
00320     // For each allocated block ---
00321 
00322     int nextBlockIdx = 0;
00323     hsize_t offset[2];
00324     hsize_t count[2];
00325     herr_t status;
00326 
00327     for (b = field->m_blocks.begin(); b != field->m_blocks.end(); ++b) {
00328       if (b->isAllocated) {
00329         offset[0] = nextBlockIdx;  // Index of next block
00330         offset[1] = 0;             // Index of first data in block. Always 0
00331         count[0] = 1;              // Number of columns to read. Always 1
00332         count[1] = valuesPerBlock; // Number of values in one column
00333         status = H5Sselect_hyperslab(fileDataSpace.id(), H5S_SELECT_SET, 
00334                                      offset, NULL, count, NULL);
00335         if (status < 0) {
00336           throw WriteHyperSlabException(
00337             "Couldn't select slab " + 
00338             boost::lexical_cast<std::string>(nextBlockIdx));
00339         }
00340         Data_T *data = &b->data[0];
00341         status = H5Dwrite(dataSet.id(), DataTypeTraits<Data_T>::h5type(), 
00342                           memDataSpace.id(), 
00343                           fileDataSpace.id(), H5P_DEFAULT, data);
00344         if (status < 0) {
00345           throw WriteHyperSlabException(
00346             "Couldn't write slab " + 
00347             boost::lexical_cast<std::string>(nextBlockIdx));
00348         }
00349         // Increment nextBlockIdx
00350         nextBlockIdx++;
00351       }
00352     }
00353 
00354   } // if occupiedBlocks > 0
00355 
00356   return true; 
00357 
00358 }
00359 
00360 //----------------------------------------------------------------------------//
00361 
00362 template <class Data_T>
00363 bool SparseFieldIO::readData(hid_t location, 
00364                              int numBlocks, 
00365                              const std::string &filename, 
00366                              const std::string &layerPath, 
00367                              typename SparseField<Data_T>::Ptr result)
00368 {
00369   using namespace std;
00370   using namespace Exc;
00371   using namespace Hdf5Util;
00372   using namespace Sparse;
00373 
00374   int occupiedBlocks;
00375 
00376   bool dynamicLoading = SparseFileManager::singleton().doLimitMemUse();
00377 
00378   int components = FieldTraits<Data_T>::dataDims();
00379   int valuesPerBlock = (1 << (result->m_blockOrder * 3)) * components;
00380   
00381   // Read the number of occupied blocks ---
00382 
00383   if (!readAttribute(location, k_numOccupiedBlocksStr, 1, occupiedBlocks)) 
00384     throw MissingAttributeException("Couldn't find attribute: " +
00385                                     k_numOccupiedBlocksStr);
00386 
00387   // Set up the dynamic read info ---
00388 
00389   if (dynamicLoading) {
00390       // Set up the field reference
00391     result->addReference(filename, layerPath,
00392                          valuesPerBlock,
00393                          occupiedBlocks);
00394   }
00395 
00396   // Read the block info data sets ---
00397 
00398   // ... Read the isAllocated array
00399 
00400   {
00401     vector<char> isAllocated(numBlocks);
00402     vector<char>::iterator i = isAllocated.begin();
00403     readSimpleData<char>(location, "block_is_allocated_data", isAllocated);
00404     typename vector<SparseBlock<Data_T> >::iterator b =
00405       result->m_blocks.begin();
00406     typename vector<SparseBlock<Data_T> >::iterator bend = 
00407       result->m_blocks.end();
00408     // We're assuming there are as many blocks in isAllocated as in the field.
00409     for (; b != bend; ++b, ++i) {
00410       b->isAllocated = static_cast<bool>(*i);
00411       if (*i && !dynamicLoading) {
00412         b->data.resize(valuesPerBlock);
00413       }
00414     }
00415   }
00416 
00417   // ... Read the emptyValue array ---
00418 
00419   {
00420     vector<Data_T> emptyValue(numBlocks);
00421     readSimpleData<Data_T>(location, "block_empty_value_data", emptyValue);
00422     typename vector<SparseBlock<Data_T> >::iterator b =
00423       result->m_blocks.begin();
00424     typename vector<SparseBlock<Data_T> >::iterator bend = 
00425       result->m_blocks.end();
00426     typename vector<Data_T>::iterator i = emptyValue.begin();
00427     // We're assuming there are as many blocks in isAllocated as in the field.
00428     for (; b != bend; ++b, ++i) {
00429       b->emptyValue = *i;
00430     }
00431   }
00432 
00433   // Read the data ---
00434 
00435   if (occupiedBlocks > 0) {
00436 
00437     if (dynamicLoading) {
00438 
00439       result->setupReferenceBlocks();
00440 
00441     } else {
00442       
00443       typename vector<SparseBlock<Data_T> >::iterator b =
00444         result->m_blocks.begin();
00445       typename vector<SparseBlock<Data_T> >::iterator bend =
00446         result->m_blocks.end();
00447 
00448       SparseDataReader<Data_T> reader(location, valuesPerBlock, occupiedBlocks);
00449 
00450       // We'll read at most 50meg at a time
00451       static const long maxMemPerPass = 50*1024*1024;
00452 
00453       for (int nextBlockIdx = 0;;) {
00454 
00455         long mem = 0;
00456         std::vector<Data_T*> memoryList;
00457         
00458         for (; b != bend && mem < maxMemPerPass; ++b) {
00459           if (b->isAllocated) {
00460             mem += sizeof(Data_T)*valuesPerBlock;
00461             memoryList.push_back(&b->data[0]);
00462           }
00463         }
00464 
00465         // all done.
00466         if (!memoryList.size()) {
00467           break;
00468         }
00469 
00470         reader.readBlockList(nextBlockIdx, memoryList);
00471         nextBlockIdx += memoryList.size();
00472       }                           
00473 
00474     }
00475 
00476   } // if occupiedBlocks > 0
00477 
00478   return true;
00479   
00480 }
00481 
00482 //----------------------------------------------------------------------------//
00483 
00484 FIELD3D_NAMESPACE_HEADER_CLOSE
00485 
00486 //----------------------------------------------------------------------------//
00487 
00488 #endif