portablepixelmap.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef PORTABLEPIXELMAP_HH
#define PORTABLEPIXELMAP_HH
 
#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
#include <vector>
#include <array>
#include <map>
#include <utility>
 
 
#include <dune/common/fvector.hh>
 
#include <dune/istl/bvector.hh>
 
#include <dune/grid/yaspgrid.hh>
 
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
 
class ColorMap
{
    public:
        enum ColorMapFileFormat {
            CM
        };
 
        enum IODataFormat {
            ASCII,
            Bin
        };
 
//        typedef unsigned short int CType;
        typedef unsigned char CType;
        typedef unsigned short int KeyType;
 
        typedef Dune::FieldVector<CType, 3> Color;
 
        typedef std::map<KeyType, Color> ColorKey;
 
        ColorMap()
        {
            key_filled_ = false;
        }
 
        ColorMap(const std::string filename)
        {
            std::string str;
            std::ifstream file;
 
            file.open(filename.c_str());
 
            std::getline(file, str);
 
            if (str.find("dune-fufem colorkey")!=std::string::npos)
            {
                char char_arr[256];
 
                /* get rid of comment and empty lines */
                do {
                    file.getline(char_arr,256);
                    str = char_arr;                 // I do not use the method std::getline(std::istream,std::string) because it does not set the istream-gcount which will be needed later. Therefore I need to do this ugly conversion
                    str.erase(0,str.find_first_not_of(" \n\t"));
                } while (str[0] == 35 or str.empty()); // ASCII(35) = #
 
                file.seekg(-file.gcount(),std::ios_base::cur);
 
                unsigned short int value=0;
                Color keycolor(0);
                KeyType key=0;
 
                while (true)
                {
                    file >> std::skipws >> key;
                    for (int j=0; j<3; ++j)
                    {
                        file >> std::skipws >> value;
                        keycolor[j] = value;
                    }
                    insertInColorKey(key, keycolor);
                    if (file.eof() or file.fail())
                        break;
                }
            }
            else
                DUNE_THROW(Dune::Exception, "unknown file format");
        }
 
        bool insertInColorKey(const KeyType key, const Color color)
        {
            bool inserted = colorkey_.insert(std::make_pair(key,color)).second;
            key_filled_ = ( inserted ? false : key_filled_ );
            return inserted;
        }
 
        Color map(const double fMin, const double fMax, const double fVal) const
        {
            if (colorkey_.size()<2)
                DUNE_THROW(Dune::Exception, "No valid ColorKey has been specified");
 
            Dune::FieldVector<double,3> dcol(0.0);
            KeyType keymax = colorkey_.rbegin()->first;
            double keyval = (fVal - fMin)/(fMax - fMin) * keymax;
            if (keyval>keymax)
                keyval = keymax;
            else if (keyval < 0)
                keyval = 0.0;
 
            std::pair<CType,Color>   col2 = *colorkey_.lower_bound(floor(keyval + 0.5)),
                                    col1 = col2;
            if (col2.first!=floor(keyval+0.5)) col1 = *(--colorkey_.lower_bound(floor(keyval+ 0.5)));
 
            double intp = ( col2.first==col1.first ? 0 : (keyval - col1.first)/(col2.first - col1.first) );
 
            for (int i=0; i<3; ++i)
                dcol[i] = (1-intp)*col1.second[i] + intp*col2.second[i];
 
            Color ret_col(0.0);
 
            for (int i = 0; i<3; ++i)
                ret_col[i] = floor(dcol[i]+0.5);
 
 
            return ret_col;
        }
 
        double inv_map(const double fMin, const double fMax, const Color color)
        {
            if (colorkey_.size()<2)
                DUNE_THROW(Dune::Exception, "No valid ColorKey has been specified");
 
            if (not key_filled_)
                fill_in_colorkey();
 
            ColorKey::iterator it = colorkey_.begin();
            ColorKey::iterator endit = colorkey_.end();
 
            KeyType keymax = colorkey_.rbegin()->first;
 
            ColorKey::iterator keycolor;
            double mindiff = 1000.0;
 
            Dune::FieldVector<int,3> diff_col;
 
            for (;it != endit; ++it)
            {
                for (int j =0; j<3; ++j)
                    diff_col[j] = int(color[j])-int(it->second[j]);
 
                if (diff_col.one_norm() < mindiff)
                {
                    mindiff = diff_col.one_norm();
                    keycolor = it;
 
                    if (mindiff == 0)
                        return fMin + keycolor->first*(fMax-fMin)/keymax;
                }
            }
 
            Dune::FieldVector<int,3> difff_col;
            for (int j =0; j<3; ++j)
            {
                diff_col[j] = (keycolor->first != 0 ? int(color[j]) - int(colorkey_[(keycolor->first)-1][j]) : 255);
                difff_col[j] = (keycolor->first != 255 ? int(color[j]) - int(colorkey_[(keycolor->first)+1][j]) : 255);
            }
 
            double nextdiff=0.0;
            int modifier=0;
            if (diff_col.one_norm() <= difff_col.one_norm())
            {
                nextdiff = diff_col.one_norm();
                modifier = -1;
            }
            else
            {
                nextdiff = difff_col.one_norm();
                modifier = 1;
            }
 
            double intp = mindiff/(mindiff+nextdiff);
 
 
            return fMin + ((1 - intp)*keycolor->first + intp*(keycolor->first + modifier))*(fMax-fMin)/keymax;
        }
 
        void exportColorMapGraphic(const std::string filename="colormap.ppm", IODataFormat df=Bin) const
        {
            unsigned short int  width = 256,
                                height = 25;
 
            std::ofstream* ppmfile;
            if (df==Bin)
                ppmfile = new std::ofstream(filename.c_str(), std::ios::trunc | std::ios::out | std::ios::binary);
            else
                ppmfile = new std::ofstream(filename.c_str(), std::ios::trunc | std::ios::out);
 
            /* write preliminaries ("P3"/"P6", some comment, height, width, maxval) to stream */
            (*ppmfile) << (df==Bin ? "P6" : "P3") << std::endl << "# file created by PortablePixelMap (dune-fufem)" << std::endl << "# colormap" << std::endl << width << " " << height << std::endl << 255 << std::endl;
 
            for (int j= 0; j<height; ++j)
            {
                for (int i= 0; i<width; ++i)
                {
                    /* we transform the function value to the corresponding point in [0,255]^3 */
                    Color color = map(0, width-1, i);
 
                    /* write color to stream - take care: << operator for char writes the corresponding ascii character */
                    for (int c=0; c<3; ++c)
                        if (df==Bin)
                            ppmfile->write(reinterpret_cast<char*>(&(color[c])),sizeof(CType));
                        else
                             (*ppmfile) << int(color[c]) << " ";
 
                    if (df==ASCII)
                        (*ppmfile) << std::endl;
                }
            }
 
            /* close PPM */
            ppmfile->close();
            delete ppmfile;
            ppmfile=0;
        }
 
        void exportColorMap(std::string filename="colormap", ColorMapFileFormat format=CM)
        {
            if (format == CM)
            {
                std::string extension = (filename.rfind(".")!=std::string::npos ? filename.substr(filename.rfind(".")) : "");
                if (extension != ".cm")
                    filename.append(".cm");
 
                std::ofstream cmfile(filename.c_str(), std::ios::trunc | std::ios::out);
 
                /* write preliminaries to stream */
                cmfile << "dune-fufem colorkey format" << std::endl << "# file created by ColorMap (dune-fufem)" << std::endl;
 
                ColorKey::const_iterator it = colorkey_.begin();
                ColorKey::const_iterator endit = colorkey_.end();
 
                for (; it != endit; ++it)
                {
                    /* write key to stream */
                    cmfile << it->first << " ";
                    /* write color to stream - take care: << operator for char writes the corresponding ascii character */
                    for (int c=0; c<3; ++c)
                        cmfile << short(it->second[c]) << (c<2?" ":"");
                    cmfile << std::endl;
                }
 
                /* close CM */
                cmfile.close();
            }
        }
 
        int size()
        {
            return colorkey_.size();
        }
    private:
        static const CType MAXVAL_=255;
 
        ColorKey colorkey_;
 
        bool key_filled_;
 
        void fill_in_colorkey()
        {
            ColorKey::iterator it = colorkey_.begin();
            ColorKey::iterator endit = colorkey_.end();
            ColorKey::iterator nextit = it;
            ++nextit;
 
            ColorKey new_elements;
            ColorKey::iterator new_elements_pos = new_elements.end();
 
            while (nextit != endit)
            {
                for (KeyType i = it->first+1; i < nextit->first; ++i)
                {
                    KeyType spread = nextit->first - it->first;
                    double frac = double(i - it->first)/spread;
 
                    Dune::FieldVector<double,3> db_col(0.0);
                    for (int j=0; j<3; ++j)
                        db_col[j] = (1-frac)*it->second[j] + frac*nextit->second[j];
 
                    Color new_color;
                    new_color[0] = floor(db_col[0]+0.5);
                    new_color[1] = floor(db_col[1]+0.5);
                    new_color[2] = floor(db_col[2]+0.5);
                    if (new_elements_pos!=new_elements.end())
                        new_elements_pos = new_elements.insert(++new_elements_pos, std::make_pair(i,new_color));
                    else
                        new_elements_pos = (new_elements.insert(std::make_pair(i,new_color))).first;
                }
                it = nextit;
                ++nextit;
            }
 
            colorkey_.insert(new_elements.begin(),new_elements.end());
            key_filled_ = true;
        }
};
 
class PortablePixelMap
{
    public:
        typedef Dune::YaspGrid<2> GridType;
 
        typedef GridType::Codim<0>::Geometry::LocalCoordinate LocalDomainType;
        typedef GridType::Codim<0>::Entity Element;
 
        using DomainType = typename GridType::template Codim<0>::Geometry::GlobalCoordinate;
        using RangeType = Dune::FieldVector<double, 1>;
 
        using Basis = Dune::Functions::LagrangeBasis<GridType::LeafGridView, 1>;
 
        using CoeffType = Dune::BlockVector<RangeType>;
        using FunctionType = Dune::Functions::DiscreteGlobalBasisFunction<Basis,CoeffType,Dune::Functions::HierarchicNodeToRangeMap,RangeType>;
 
        using LocalFunction = typename FunctionType::LocalFunction;
        using EntitySet = typename FunctionType::EntitySet;
 
 
        enum IODataFormat {
            ASCII,
            Bin
        };
 
    private:
        typedef ColorMap::Color Color;
 
        unsigned int
            width_,
            height_;
 
        DomainType::field_type
            xDomainMin_,
            xDomainMax_,
            yDomainMin_,
            yDomainMax_;
 
        const RangeType
            fRangeMin_,
            fRangeMax_;
 
        std::vector<unsigned int> imgParams_;
 
        ColorMap colormap_;
 
        Basis* basis_;
        CoeffType coeffs_;
        FunctionType* pixMapFunction_;
        GridType* grid_;
 
        using DomainFieldType = typename Dune::template FieldTraits<DomainType>::field_type;
 
        template<class Function>
        using RangeFieldTypeFor = typename Dune::template FieldTraits<std::decay_t<std::invoke_result_t<Function, DomainType>>>::field_type;
 
    public:
        PortablePixelMap(const DomainType::field_type xDomainMin, const DomainType::field_type xDomainMax,
                         const DomainType::field_type yDomainMin, const DomainType::field_type yDomainMax,
                         const RangeType fRangeMin, const RangeType fRangeMax,
                         const ColorMap colormap):
            xDomainMin_(xDomainMin),
            xDomainMax_(xDomainMax),
            yDomainMin_(yDomainMin),
            yDomainMax_(yDomainMax),
            fRangeMin_(fRangeMin),
            fRangeMax_(fRangeMax),
            colormap_(colormap)
        {
            if (fRangeMax<=fRangeMin)
                DUNE_THROW(Dune::RangeError, "Range of function values has inverse bounds or zero length.");
        }
 
        PortablePixelMap(const RangeType fRangeMin, const RangeType fRangeMax, const ColorMap colormap):
            xDomainMin_(0),
            xDomainMax_(0),
            yDomainMin_(0),
            yDomainMax_(0),
            fRangeMin_(fRangeMin),
            fRangeMax_(fRangeMax),
            colormap_(colormap),
            grid_(NULL)
        {
            if (fRangeMax<=fRangeMin)
                DUNE_THROW(Dune::RangeError, "Range of function values has inverse bounds or zero length.");
        }
 
        void readPixelMap(const char* filename, const int width=-1, const int height=-1, const int xoffset=0, const int yoffset=0)
        {
            std::string str;
            std::string ppmID;
            std::ifstream ppmfile;
 
            ppmfile.open(filename);
 
            /* get rid of leading comment lines */
            do {
                std::getline(ppmfile, ppmID);
            } while (ppmID[0] == 35); // ASCII(35) = #
 
            if (ppmID!="P3" and ppmID!="P6")
                DUNE_THROW(Dune::Exception,"File is not a PPM-File.");
 
            int count = 0;
            while (count < 3)
            {
                ppmfile >> std::skipws >> str;
                if (str[0] == 35)
                {
                    ppmfile.ignore(1024,10);
                    continue;
                }
                else
                {
                    imgParams_.push_back(std::atoi(str.c_str()));
                    ++count;
                }
            }
 
            int imgWidth = imgParams_[0],
                imgHeight = imgParams_[1],
                maxval_ = imgParams_[2];
 
            unsigned short int datalength = maxval_ > 255 ? 2 : 1;
            if (datalength==2)
                DUNE_THROW(Dune::Exception, "At the moment only PixelMaps with 8Bit/color channel, i.e. a maxval of 255 may be processed due to a restriction of ColorMap::CType to unsigned char.");
 
            count = 0;
 
 
            width  < 0 ? width_  = imgWidth  : width_  = width;
            height < 0 ? height_ = imgHeight : height_ = height;
 
            if(xDomainMax_-xDomainMin_ <= 0.0)
            {
                xDomainMin_ = 0;
                xDomainMax_ = width_-1;
 
                yDomainMin_ = 0;
                yDomainMax_ = height_-1;
            }
 
            unsigned short int value=0;
 
            int p1 = ceil(log(width_ -1)/log(2)),
                p2 = ceil(log(height_-1)/log(2)),
                i = (p1-p2 >= 0)? p1-p2 : 0,
                j = (p1-p2 >= 0)? 0 : p2-p1;
 
            Dune::FieldVector<double,2> L;
            std::array<int,2> s;
 
            L[0] = pow(2, p1);
            L[1] = pow(2, p2);
            s[0] = pow(2, i);
            s[1] = pow(2, j);
 
            grid_ = new GridType(L,s);
            for (int k=0; k<std::min(p1,p2); ++k)
                grid_->globalRefine(1);
 
            basis_ = new Basis(grid_->leafGridView());
            coeffs_.resize(basis_->size());
            coeffs_ = 0.0;
            pixMapFunction_ = new FunctionType(Dune::Functions::makeDiscreteGlobalBasisFunction<RangeType>(*basis_, coeffs_));
 
            /* if it's the binary format we need to open in binary mode */
            if (ppmID == "P6")
            {
                ppmfile.close();
                ppmfile.open(filename,std::ios::binary);
            }
 
            /* Skip yoffset pixel lines and xoffset pixel columns */
            for (int k=0; k < 3*(yoffset*imgWidth+xoffset); ++k)
                if (ppmID=="P6")
                    ppmfile.read(reinterpret_cast<char*>(&value), datalength*sizeof(char));
                else
                    ppmfile >> std::skipws >> value;
 
            /* read width x height color values */
            for (std::size_t i=0; i<height_; ++i)
            {
                for(std::size_t j=0; j<width_; ++j)
                {
                    Color pix_color(0);
                    for (std::size_t c=0; c<3; ++c)
                    {
                        if (ppmID == "P6")
                            ppmfile.read(reinterpret_cast<char*>(&(pix_color[c])), datalength*sizeof(char));
                        else
                        {
                            ppmfile >> std::skipws >> value;
                            pix_color[c] = value;
                        }
 
                        if (ppmfile.eof())
                            break;
 
                        pix_color[c] = value;
                    }
                    coeffs_[(height_-1-i)*(L[0]+1)+j] = colormap_.inv_map(fRangeMin_, fRangeMax_, pix_color);
                }
 
                /* Skip rest of pixel line */
                for (std::size_t k=0; k < 3*(imgWidth-width_); ++k)
                    if (ppmID=="P6")
                        ppmfile.read(reinterpret_cast<char*>(&value), datalength*sizeof(char));
                    else
                        ppmfile >> std::skipws >> value;
            }
 
        }
 
        template <typename FType>
        static void exportPixelMap( const char* filename,
                                    const FType& function,
                                    RangeFieldTypeFor<FType> fRangeMin, RangeFieldTypeFor<FType> fRangeMax,
                                    const DomainFieldType xDomainMin, const DomainFieldType xDomainMax,
                                    const DomainFieldType yDomainMin, const DomainFieldType yDomainMax,
                                    const unsigned int width, const unsigned int height, 
                                    const std::string info,
                                    const ColorMap colormap, IODataFormat df=Bin)
        {
            using RangeType = std::decay_t<std::invoke_result_t<FType, DomainType>>;
            using RangeFieldType = RangeFieldTypeFor<FType>;
 
            DomainFieldType xLength = xDomainMax - xDomainMin;
            DomainFieldType yLength = yDomainMax - yDomainMin;
 
            RangeFieldType rangeInterval = fRangeMax-fRangeMin;
 
            if (xLength <= 0.0 or yLength <= 0.0 or rangeInterval <= 0.0)
                DUNE_THROW(Dune::RangeError,"Domain has zero or negative extension or range interval has zero length or negative extension(thrown by PortablePixelMap::exportPixelMap.");
 
            DomainType r(0.0);
            RangeType  fval(0.0);
 
            std::ofstream* ppmfile;
 
            if (df==Bin)
                ppmfile = new std::ofstream(filename, std::ios::trunc | std::ios::out | std::ios::binary);
            else
                ppmfile = new std::ofstream(filename, std::ios::trunc | std::ios::out);
 
            // write preliminaries ("P3", some comment, height, width, maxval) to stream
            (*ppmfile) << (df==Bin ? "P6" : "P3") << std::endl << "# file created by PortablePixelMap (dune-fufem)" << std::endl << "# " << info << std::endl << width << " " << height << std::endl << 255 << std::endl;
 
            r[0] = xDomainMin;
            r[1] = yDomainMax;
 
            for (std::size_t j= 0; j<height; ++j)
            {
                r[0] = xDomainMin;
                for (std::size_t i= 0; i<width; ++i)
                {
                    fval = function(r);
 
                     // we transform the function value to the corresponding point in [0,255]^3
                    Color color = colormap.map(fRangeMin, fRangeMax, fval);
 
                    // colorval >> endl >> PPM
                    for (int c=0; c<3; ++c)
                        if (df==Bin)
                            ppmfile->write(reinterpret_cast<char*>(&(color[c])),sizeof(ColorMap::CType));
                        else
                             (*ppmfile) << int(color[c]) << " ";
 
                    if (df==ASCII)
                        (*ppmfile) << std::endl;
 
                    r[0] = xDomainMin + (i+1)*xLength/(width-1);
                }
                r[1] = yDomainMax - (j+1)*yLength/(height-1);
            }
 
            // close PPM 
            ppmfile->close();
        }
 
        /* \brief evaluates interpolated pixelmap at global point
         *
         *  The pixelmap is evaluated at point x of the domain \f$ \Omega \f$
         *  assuming a linear transformation \f$ \Omega \longrightarrow [0,width_]\times[0,height_]\f$. The returnvalue is scaled
         *  according to \f$[fRangeMin_,fRangeMax_]\f$.
         *
         * \param x     point in \f$ \Omega \f$ at which to evaluate
         * \returns     the function value
         */
        RangeType operator()(const DomainType& x) const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            DomainType xx = x;
 
            xx[0] = (x[0]-xDomainMin_)/(xDomainMax_-xDomainMin_)*(width_-1);
            xx[1] = (x[1]-yDomainMin_)/(yDomainMax_-yDomainMin_)*(height_-1);
 
            return (*pixMapFunction_)(xx);
        }
 
        friend LocalFunction localFunction(const PortablePixelMap& ppm)
        {
            if (ppm.pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
            return localFunction(*(ppm.pixMapFunction_));
        }
 
        const EntitySet& entitySet() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
            return pixMapFunction_->entitySet();
        }
 
        /* \brief provides read access to the coefficient vector
         */
        const CoeffType& getCoeffs() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            return coeffs_;
        }
 
        /* \brief provides access to the grid on which the interpolated pixelmap lives
         */
        const GridType& getGrid() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            return *grid_;
        }
 
        /* \brief provides access to the width of the imported pixelmap section
         */
        unsigned int getWidth() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            return width_;
        }
 
        /* \brief provides access to the height of the imported pixelmap section
         */
        unsigned int getHeight() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            return height_;
        }
 
        /* \brief provides access to the width of the full (original) pixelmap image
         */
        unsigned int getImgWidth() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            return imgParams_[0];
        }
 
        /* \brief provides access to the height of the full (original) pixelmap image
         */
        unsigned int getImgHeight() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            return imgParams_[1];
        }
 
        /* \brief provides access to the used maximal value in the (original and imported) pixelmap corresponding to white
         */
        unsigned int getWhiteValue() const
        {
            if (pixMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No pixelmap has been loaded. Use readPixelMap() first.");
 
            return imgParams_[2];
        }
 
        void exportColormapGraphic(const std::string filename="colormap.ppm") const
        {
            colormap_.exportColorMapGraphic(filename);
        }
 
        ~PortablePixelMap()
        {
            if (grid_!=0)
                delete grid_;
            if (pixMapFunction_!=0)
                delete pixMapFunction_;
            if (basis_!=0)
                delete basis_;
        }
};
 
#endif