polaroid-pp

Schlieren and contour plot tool
git clone https://git.0xfab.ch/polaroid-pp.git
Log | Files | Refs | Submodules | README | LICENSE

Polaroid.cpp (12644B)

      1 // File       : Polaroid.cpp
      2 // Date       : Tue Apr 26 22:14:08 2016
      3 // Author     : Fabian Wermelinger
      4 // Description: Polaroid Implementation
      5 // Copyright 2016 ETH Zurich. All Rights Reserved.
      6 #include <cassert>
      7 #include <cmath>
      8 #include <cstring>
      9 #include <algorithm>
     10 #include <string>
     11 #include <sstream>
     12 #include <cstdlib>
     13 #ifdef _USE_HDF_
     14 #include <hdf5.h>
     15 #endif /* _USE_HDF_ */
     16 
     17 #include "Polaroid.h"
     18 #ifdef _USE_CUBISMZ_
     19 #include "Reader_WaveletCompression.h"
     20 #endif /* _USE_CUBISMZ_ */
     21 
     22 using namespace std;
     23 
     24 #define ID3(ix, iy, iz, NX, NY) ((ix) + (NX) * ((iy) + (NY) * (iz)))
     25 
     26 template <typename T>
     27 inline void _streamSwap(const size_t n, unsigned char * const stream)
     28 {
     29     constexpr size_t _e = sizeof(T);
     30     unsigned char _buf[_e];
     31     for (size_t i = 0; i < n*_e; i += _e)
     32     {
     33         const unsigned char * const upstream = stream + i + (_e - 1);
     34 #pragma unroll(_e)
     35         for (int j = 0; j < _e; ++j)
     36             _buf[j] = *(upstream - j);
     37         std::memcpy(stream+i, &_buf[0], _e);
     38     }
     39 }
     40 
     41 void Polaroid::load_hdf5(const char* filename, ArgumentParser& parser)
     42 {
     43 #ifdef _USE_HDF_
     44     const double sf      = parser("-sf").asDouble(0.5); // slice fraction (default)
     45     const int si         = parser("-si").asInt(-1); // slice index (need to specify)
     46     const char direction = parser("-sd").asString("x")[0]; // slice normal direction
     47     const int channel    = parser("-channel").asInt(0); // data channel
     48     const bool magnitude = parser("-magnitude").asBool(false); // vector magnitude (only if NCH == 3)
     49     const bool swap      = parser("-swap").asBool(false); // swap bytes?
     50 
     51     /* open data */
     52     hid_t file_id, dataset_id, dataspace_id, file_dataspace_id;
     53     hsize_t* dims;
     54     hssize_t num_elem;
     55     int rank, ndims, NCH;
     56     int maxDim[3];
     57     Real* data;
     58 
     59     file_id           = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
     60     dataset_id        = H5Dopen(file_id, "/data", H5P_DEFAULT);
     61     file_dataspace_id = H5Dget_space(dataset_id);
     62     rank              = H5Sget_simple_extent_ndims(file_dataspace_id);
     63     dims              = new hsize_t[rank];
     64     ndims             = H5Sget_simple_extent_dims(file_dataspace_id, dims, NULL);
     65     num_elem          = H5Sget_simple_extent_npoints(file_dataspace_id);
     66     data              = new Real[num_elem];
     67     maxDim[2]         = dims[0];
     68     maxDim[1]         = dims[1];
     69     maxDim[0]         = dims[2];
     70     NCH               = dims[3];
     71     dataspace_id      = H5Screate_simple(rank, dims, NULL);
     72     int status        = H5Dread(dataset_id, HDF_PRECISION, dataspace_id, file_dataspace_id, H5P_DEFAULT, data);
     73 
     74     const int Nmax = std::max(maxDim[0], std::max(maxDim[1], maxDim[2]));
     75     m_data.set_max3Ddim(Nmax);
     76 
     77     /* release stuff */
     78     delete [] dims;
     79     status = H5Dclose(dataset_id);
     80     status = H5Sclose(dataspace_id);
     81     status = H5Sclose(file_dataspace_id);
     82     status = H5Fclose(file_id);
     83 
     84     assert(channel < NCH);
     85 
     86     /* extract plane */
     87     m_data.set_sliceDir(direction);
     88     if ('x' == direction) // y-z plane
     89     {
     90         const int fixed = (si >= 0) ? si : static_cast<int>(maxDim[0]*sf);
     91         m_data.set_sliceID(fixed);
     92         m_data.resize(maxDim[1], maxDim[2]);
     93         if (magnitude && NCH == 3)
     94         {
     95             for (int h=0; h < m_data.height(); ++h)
     96                 for (int w=0; w < m_data.width(); ++w)
     97                 {
     98                     const Real a = data[0 + NCH * ID3(fixed,w,h,maxDim[0], maxDim[1])];
     99                     const Real b = data[1 + NCH * ID3(fixed,w,h,maxDim[0], maxDim[1])];
    100                     const Real c = data[2 + NCH * ID3(fixed,w,h,maxDim[0], maxDim[1])];
    101                     m_data(w,h) = std::sqrt(a*a + b*b + c*c);
    102                 }
    103         }
    104         else
    105         {
    106             for (int h=0; h < m_data.height(); ++h)
    107                 for (int w=0; w < m_data.width(); ++w)
    108                     m_data(w,h) = data[channel + NCH * ID3(fixed,w,h,maxDim[0], maxDim[1])];
    109         }
    110     }
    111     else if ('y' == direction) // x-z plane
    112     {
    113         const int fixed = (si >= 0) ? si : static_cast<int>(maxDim[1]*sf);
    114         m_data.set_sliceID(fixed);
    115         m_data.resize(maxDim[0], maxDim[2]);
    116         if (magnitude && NCH == 3)
    117         {
    118             for (int h=0; h < m_data.height(); ++h)
    119                 for (int w=0; w < m_data.width(); ++w)
    120                 {
    121                     const Real a = data[0 + NCH * ID3(w,fixed,h,maxDim[0], maxDim[1])];
    122                     const Real b = data[1 + NCH * ID3(w,fixed,h,maxDim[0], maxDim[1])];
    123                     const Real c = data[2 + NCH * ID3(w,fixed,h,maxDim[0], maxDim[1])];
    124                     m_data(w,h) = std::sqrt(a*a + b*b + c*c);
    125                 }
    126         }
    127         else
    128         {
    129             for (int h=0; h < m_data.height(); ++h)
    130                 for (int w=0; w < m_data.width(); ++w)
    131                     m_data(w,h) = data[channel + NCH * ID3(w,fixed,h,maxDim[0], maxDim[1])];
    132         }
    133     }
    134     else if ('z' == direction) // x-y plane
    135     {
    136         const int fixed = (si >= 0) ? si : static_cast<int>(maxDim[2]*sf);
    137         m_data.set_sliceID(fixed);
    138         m_data.resize(maxDim[0], maxDim[1]);
    139         if (magnitude && NCH == 3)
    140         {
    141             for (int h=0; h < m_data.height(); ++h)
    142                 for (int w=0; w < m_data.width(); ++w)
    143                 {
    144                     const Real a = data[0 + NCH * ID3(w,h,fixed,maxDim[0], maxDim[1])];
    145                     const Real b = data[1 + NCH * ID3(w,h,fixed,maxDim[0], maxDim[1])];
    146                     const Real c = data[2 + NCH * ID3(w,h,fixed,maxDim[0], maxDim[1])];
    147                     m_data(w,h) = std::sqrt(a*a + b*b + c*c);
    148                 }
    149         }
    150         else
    151         {
    152             for (int h=0; h < m_data.height(); ++h)
    153                 for (int w=0; w < m_data.width(); ++w)
    154                     m_data(w,h) = data[channel + NCH * ID3(w,h,fixed,maxDim[0], maxDim[1])];
    155         }
    156     }
    157     else
    158     {
    159         fprintf(stderr, "ERROR: Direction '%c' unsupported.", direction);
    160         abort();
    161     }
    162 
    163     delete [] data;
    164     m_dataLoaded = true;
    165 
    166 #else
    167     fprintf(stderr, "WARNING: Executable was compiled without HDF support...\n");
    168 #endif /* _USE_HDF_ */
    169 }
    170 
    171 void Polaroid::load_hdf5_slice(const char* filename, ArgumentParser& parser)
    172 {
    173 #ifdef _USE_HDF_
    174     const int channel    = parser("-channel").asInt(0); // data channel
    175     const bool magnitude = parser("-magnitude").asBool(false); // vector magnitude (only if NCH == 3)
    176     const bool swap      = parser("-swap").asBool(false); // swap bytes?
    177 
    178     /* open data */
    179     hid_t file_id, dataset_id, dataspace_id, file_dataspace_id;
    180     hsize_t* dims;
    181     hssize_t num_elem;
    182     int rank, ndims, NCH;
    183     int maxDim[2];
    184     Real* data;
    185 
    186     file_id           = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
    187     dataset_id        = H5Dopen(file_id, "/data", H5P_DEFAULT);
    188     file_dataspace_id = H5Dget_space(dataset_id);
    189     rank              = H5Sget_simple_extent_ndims(file_dataspace_id);
    190     dims              = new hsize_t[rank];
    191     ndims             = H5Sget_simple_extent_dims(file_dataspace_id, dims, NULL);
    192     num_elem          = H5Sget_simple_extent_npoints(file_dataspace_id);
    193     data              = new Real[num_elem];
    194     maxDim[1]         = dims[0];
    195     maxDim[0]         = dims[1];
    196     NCH               = dims[2];
    197     dataspace_id      = H5Screate_simple(rank, dims, NULL);
    198     int status        = H5Dread(dataset_id, HDF_PRECISION, dataspace_id, file_dataspace_id, H5P_DEFAULT, data);
    199 
    200     const int Nmax = std::max(maxDim[0], maxDim[1]);
    201     m_data.set_max3Ddim(Nmax);
    202 
    203     /* release stuff */
    204     delete [] dims;
    205     status = H5Dclose(dataset_id);
    206     status = H5Sclose(dataspace_id);
    207     status = H5Sclose(file_dataspace_id);
    208     status = H5Fclose(file_id);
    209 
    210     assert(channel < NCH);
    211 
    212     /* extract plane */
    213     m_data.set_sliceDir('s'); // from 2D slice
    214     m_data.set_sliceID(-1);
    215     m_data.resize(maxDim[0], maxDim[1]);
    216     if (magnitude && NCH == 3)
    217     {
    218         for (int h=0; h < m_data.height(); ++h)
    219             for (int w=0; w < m_data.width(); ++w)
    220             {
    221                 const Real a = data[0 + NCH*(w + maxDim[0]*h)];
    222                 const Real b = data[1 + NCH*(w + maxDim[0]*h)];
    223                 const Real c = data[2 + NCH*(w + maxDim[0]*h)];
    224                 m_data(w,h) = std::sqrt(a*a + b*b + c*c);
    225             }
    226     }
    227     else
    228     {
    229         for (int h=0; h < m_data.height(); ++h)
    230             for (int w=0; w < m_data.width(); ++w)
    231                 m_data(w,h) = data[channel + NCH*(w + maxDim[0]*h)];
    232     }
    233     delete [] data;
    234     m_dataLoaded = true;
    235 
    236 #else
    237     fprintf(stderr, "WARNING: Executable was compiled without HDF support...\n");
    238 #endif /* _USE_HDF_ */
    239 }
    240 
    241 void Polaroid::load_wavelet(const char* filename, ArgumentParser& parser)
    242 {
    243 #ifdef _USE_CUBISMZ_
    244     const double sf        = parser("-sf").asDouble(0.5); // slice fraction (default)
    245     const int si           = parser("-si").asInt(-1); // slice index (need to specify)
    246     const char direction   = parser("-sd").asString("x")[0]; // slice normal direction
    247     const bool byte_swap   = parser("-swap").asBool(false);
    248     const int wavelet_type = parser("-wtype").asInt(1);
    249 
    250     const string fname(filename);
    251     Reader_WaveletCompression myreader(fname, byte_swap, wavelet_type);
    252     myreader.load_file();
    253     const int NBX = myreader.xblocks();
    254     const int NBY = myreader.yblocks();
    255     const int NBZ = myreader.zblocks();
    256     const int maxDim[3] = {NBX*_BLOCKSIZE_, NBY*_BLOCKSIZE_, NBZ*_BLOCKSIZE_};
    257     const int Nmax = std::max(maxDim[0], std::max(maxDim[1], maxDim[2]));
    258     m_data.set_max3Ddim(Nmax);
    259 
    260     Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_];
    261 
    262     /* extract plane */
    263     m_data.set_sliceDir(direction);
    264     if ('x' == direction) // y-z plane
    265     {
    266         const int fixed = (si >= 0) ? si : static_cast<int>(maxDim[0]*sf);
    267         m_data.set_sliceID(fixed);
    268         m_data.resize(maxDim[1], maxDim[2]);
    269 
    270         const int fixedBID = fixed/_BLOCKSIZE_;
    271         const int BlocalID = fixed % _BLOCKSIZE_;
    272         for (int iz=0; iz < NBZ; ++iz)
    273             for (int iy=0; iy < NBY; ++iy)
    274             {
    275                 const double zratio = myreader.load_block2(fixedBID, iy, iz, blockdata);
    276                 for (int z=0; z < _BLOCKSIZE_; ++z)
    277                     for (int y=0; y < _BLOCKSIZE_; ++y)
    278                     {
    279                         assert(iy*_BLOCKSIZE_+y < m_data.width());
    280                         assert(iz*_BLOCKSIZE_+z < m_data.height());
    281                         m_data(iy*_BLOCKSIZE_+y, iz*_BLOCKSIZE_+z) = blockdata[z][y][BlocalID];
    282                     }
    283             }
    284     }
    285     else if ('y' == direction) // x-z plane
    286     {
    287         const int fixed = (si >= 0) ? si : static_cast<int>(maxDim[1]*sf);
    288         m_data.set_sliceID(fixed);
    289         m_data.resize(maxDim[0], maxDim[2]);
    290 
    291         const int fixedBID = fixed/_BLOCKSIZE_;
    292         const int BlocalID = fixed % _BLOCKSIZE_;
    293         for (int iz=0; iz < NBZ; ++iz)
    294             for (int ix=0; ix < NBX; ++ix)
    295             {
    296                 const double zratio = myreader.load_block2(ix, fixedBID, iz, blockdata);
    297                 for (int z=0; z < _BLOCKSIZE_; ++z)
    298                     for (int x=0; x < _BLOCKSIZE_; ++x)
    299                     {
    300                         assert(ix*_BLOCKSIZE_+x < m_data.width());
    301                         assert(iz*_BLOCKSIZE_+z < m_data.height());
    302                         m_data(ix*_BLOCKSIZE_+x, iz*_BLOCKSIZE_+z) = blockdata[z][BlocalID][x];
    303                     }
    304             }
    305     }
    306     else if ('z' == direction) // x-y plane
    307     {
    308         const int fixed = (si >= 0) ? si : static_cast<int>(maxDim[2]*sf);
    309         m_data.set_sliceID(fixed);
    310         m_data.resize(maxDim[0], maxDim[1]);
    311 
    312         const int fixedBID = fixed/_BLOCKSIZE_;
    313         const int BlocalID = fixed % _BLOCKSIZE_;
    314         for (int iy=0; iy < NBY; ++iy)
    315             for (int ix=0; ix < NBX; ++ix)
    316             {
    317                 const double zratio = myreader.load_block2(ix, iy, fixedBID, blockdata);
    318                 for (int y=0; y < _BLOCKSIZE_; ++y)
    319                     for (int x=0; x < _BLOCKSIZE_; ++x)
    320                     {
    321                         assert(ix*_BLOCKSIZE_+x < m_data.width());
    322                         assert(iy*_BLOCKSIZE_+y < m_data.height());
    323                         m_data(ix*_BLOCKSIZE_+x, iy*_BLOCKSIZE_+y) = blockdata[BlocalID][y][x];
    324                     }
    325             }
    326     }
    327     else
    328     {
    329         fprintf(stderr, "ERROR: Direction '%c' unsupported.", direction);
    330         abort();
    331     }
    332 
    333     m_dataLoaded = true;
    334 #else
    335     fprintf(stderr, "WARNING: Executable was compiled without wavelet compressor support...\n");
    336 #endif /* _USE_CUBISMZ_ */
    337 }