polaroid-pp

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

SchlierenCartridge.h (3315B)

      1 // File       : SchlierenCartridge.h
      2 // Date       : Thu 28 Apr 2016 04:45:20 PM CEST
      3 // Author     : Fabian Wermelinger
      4 // Description: Schlieren Image Cartridge
      5 // Copyright 2016 ETH Zurich. All Rights Reserved.
      6 #ifndef SCHLIERENCARTRIDGE_H_AUHX7ILM
      7 #define SCHLIERENCARTRIDGE_H_AUHX7ILM
      8 
      9 #include <cstdio>
     10 #include <algorithm>
     11 #include <cmath>
     12 #include <functional>
     13 #include "Cartridge.h"
     14 #include "GammaCorrection.h"
     15 
     16 class SchlierenCartridge : public Cartridge
     17 {
     18     Real m_mean;
     19     Real m_std;
     20     size_t m_k;
     21     Slice m_gradX;
     22     Slice m_gradY;
     23 
     24     // helper
     25     void _compute(Slice& data);
     26     void _gradX(const Slice& data);
     27     void _gradY(const Slice& data);
     28 
     29 public:
     30     SchlierenCartridge(ArgumentParser& parser) : Cartridge(parser), m_mean(0.0), m_std(0.0), m_k(0) {}
     31 
     32     virtual void capture(PhotoPaper& photo, Slice& data)
     33     {
     34         const Real ka0 = m_parser("-ka0").asDouble(1.0);
     35         const Real ka1 = m_parser("-ka1").asDouble(1.0);
     36         // const Real k0 = m_parser("-k0").asDouble(0.0);
     37         // const Real k1 = m_parser("-k1").asDouble(1.0);
     38         GammaCorrection gammaCorrect(m_parser);
     39 
     40         photo.make_new(photo.get_name()+"-schlieren", data.width(), data.height());
     41 
     42         // set description
     43         string desc("2D_Schlieren");
     44         photo.set_description(desc.c_str());
     45 
     46         // compute schlieren
     47         if (!m_bComputed)
     48         {
     49             _compute(data);
     50             m_dataMin = data.min();
     51             m_dataMax = data.max();
     52         }
     53         const Real dataMaxInv = 1.0/m_dataMax;
     54         const Real dataMaxInv_log = 1.0/std::log(m_dataMax+1.0);
     55 
     56         // const Real fac = -ka0/(k1 - k0);
     57         // auto exp_shader     = [&](const Real x) { return std::exp(fac*(x*dataMaxInv - k0)); };
     58 
     59         auto exp_shader     = [&](const Real x) { return std::exp(-ka0*x*dataMaxInv); };
     60 
     61         auto exp_shader_log = [&](const Real x)
     62         {
     63             const Real psi = std::log(x+1.0) * dataMaxInv_log;
     64             return std::exp(-ka1*psi);
     65         };
     66 
     67         auto blend_shaders = [&](const Real x)
     68         {
     69             const Real beta = (x-m_dataMin)/(m_dataMax-m_dataMin);
     70             return beta*exp_shader(x) + (1.0-beta)*exp_shader_log(x);
     71         };
     72 
     73         std::function<Real(const Real)> _shader = exp_shader;
     74         if (m_parser("-blend").asBool(false))
     75             _shader = blend_shaders;
     76 
     77         // apply shader
     78         for (int h=0; h < data.height(); ++h)
     79             for (int w=0; w < data.width(); ++w)
     80             {
     81                 const Real phi = _shader(data(w,h));
     82                 photo.set_pixel(gammaCorrect(phi), w, h);
     83             }
     84 
     85         photo.write();
     86     }
     87 
     88     virtual void compute(Slice& data)
     89     {
     90         _compute(data);
     91         m_dataMin = std::min(m_dataMin, data.min());
     92         m_dataMax = std::max(m_dataMax, data.max());
     93         m_bComputed = true;
     94         for (int h=0; h < data.height(); ++h)
     95             for (int w=0; w < data.width(); ++w)
     96             {
     97                 ++m_k;
     98                 const Real delta = data(w,h) - m_mean;
     99                 m_mean += delta/m_k;
    100                 m_std  += delta*(data(w,h) - m_mean);
    101             }
    102     }
    103 
    104     virtual void reset()
    105     {
    106         Cartridge::reset();
    107         m_mean = 0.0;
    108         m_std  = 0.0;
    109     }
    110 };
    111 
    112 #endif /* SCHLIERENCARTRIDGE_H_AUHX7ILM */