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