portablegreymap.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 PORTABLEGREYMAP_HH
#define PORTABLEGREYMAP_HH
 
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
#include <array>
 
#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 PortableGreyMap
{
    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 = typename Dune::Functions::LagrangeBasis<GridType::LeafGridView, 1>;
 
        using CoeffType = typename Dune::BlockVector<RangeType>;
        using FunctionType = Dune::Functions::DiscreteGlobalBasisFunction<Basis,CoeffType,Dune::Functions::HierarchicNodeToRangeMap,RangeType>;
 
        using EntitySet = typename FunctionType::EntitySet;
 
        enum ColorScheme {
                DEFAULT,    
                REVERSE     
        };
 
    private:
 
        unsigned int
            width_,
            height_;
 
        DomainType::field_type
            xDomainMin_,
            xDomainMax_,
            yDomainMin_,
            yDomainMax_;
 
        const RangeType
            fRangeMin_,
            fRangeMax_;
 
        std::vector<unsigned int> imgParams_;
 
        const ColorScheme colorscheme_;
 
        Basis* basis_;
        CoeffType coeffs_;
        FunctionType* greyMapFunction_;
        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:
        PortableGreyMap(const DomainType::field_type xDomainMin, const DomainType::field_type xDomainMax,
                        const DomainType::field_type yDomainMin, const DomainType::field_type yDomainMax,
                        const RangeType fRangeMin=0, const RangeType fRangeMax=1,
                        const ColorScheme colorscheme=DEFAULT):
            xDomainMin_(xDomainMin),
            xDomainMax_(xDomainMax),
            yDomainMin_(yDomainMin),
            yDomainMax_(yDomainMax),
            fRangeMin_(fRangeMin),
            fRangeMax_(fRangeMax),
            colorscheme_(colorscheme)
        {
            if (fRangeMax<=fRangeMin)
                DUNE_THROW(Dune::RangeError, "Range of function values has inverse bounds or zero length.");
        }
 
        PortableGreyMap(const RangeType fRangeMin=0, const RangeType fRangeMax=1, const ColorScheme colorscheme=DEFAULT):
            xDomainMin_(0),
            xDomainMax_(0),
            yDomainMin_(0),
            yDomainMax_(0),
            fRangeMin_(fRangeMin),
            fRangeMax_(fRangeMax),
            colorscheme_(colorscheme)
        {
            if (fRangeMax<=fRangeMin)
                DUNE_THROW(Dune::RangeError, "Range of function values has inverse bounds or zero length.");
        }
 
        void readGreyMap(const char* filename, const int width=-1, const int height=-1, const int xoffset=0, const int yoffset=0)
        {
            std::string str;
            std::string pgmID;
            std::ifstream pgmfile;
 
            pgmfile.open(filename);
 
            // get rid of leading comment lines
            do {
                std::getline(pgmfile, pgmID);
            } while (str[0] == 35); // ASCII(35) = #
 
            if (pgmID[0]!=80) // ASCII(80) = P
                DUNE_THROW(Dune::Exception,"File is not a PGM-File.");
 
            int count = 0;
            while (count < 3)
            {
                pgmfile >> std::skipws >> str;
                if (str[0] == 35)
                {
                    pgmfile.ignore(1024,10);
                    continue;
                }
                else
                {
                    imgParams_.push_back(std::atoi(str.c_str()));
                    ++count;
                }
            }
 
            int imgWidth = imgParams_[0],
                imgHeight = imgParams_[1];
                //maxGreyVal = imgParams_[2];
 
            count = 0;
 
 
            width_  = (width  < 0) ? imgWidth  : width;
            height_ = (height < 0) ? imgHeight : height;
 
            if(xDomainMax_-xDomainMin_ <= 0.0)
            {
                xDomainMin_ = 0;
                xDomainMax_ = width_-1;
 
                yDomainMin_ = 0;
                yDomainMax_ = height_-1;
            }
 
            int value=0;
 
            int p1 = ceil(log(width_ -1)/log(2)),
                p2 = ceil(log(height_-1)/log(2));
            unsigned int i = (p1-p2 >= 0)? p1-p2 : 0;
            unsigned int 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;
            greyMapFunction_ = new FunctionType(Dune::Functions::makeDiscreteGlobalBasisFunction<RangeType>(*basis_, coeffs_));
 
            if (pgmID == "P2")
            {
                for (int k=0; k < yoffset*imgWidth+xoffset; ++k)
                    pgmfile >> std::skipws >> value;
 
                for (i=0; i<height_; ++i)
                {
                    for(j=0; j<width_; ++j)
                    {
                        pgmfile >> std::skipws >> value;
 
                        if (pgmfile.eof())
                            break;
 
                        coeffs_[(height_-1-i)*(L[0]+1)+j] = static_cast<double>(value);
                    }
 
                    for (unsigned int k=0; k < imgWidth-width_; ++k)
                        pgmfile >> std::skipws >> value;
                }
            }
            else
                DUNE_THROW(Dune::Exception,"File is not an ASCII PGM-File.");
 
        }
 
        template <typename FType>
        static void exportGreyMap(  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 unsigned int maxGreyVal=255, const ColorScheme colorscheme=DEFAULT)
        {
            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 PortableGreyMap::exportGreyMap.");
 
            double tFactor = maxGreyVal/rangeInterval;
 
            DomainType r(0.0);
            RangeType  fval(0.0);
 
            std::ofstream pgmfile(filename, std::ios::trunc | std::ios::out);
 
            // write preliminaries ("P2", some comment, height, width, maxGreyVal) to stream
            pgmfile << "P2" << std::endl << "# file created by PortableGreyMap (dune-fufem)" << std::endl << "# " << info << std::endl << height << " " << width << std::endl << maxGreyVal << std::endl;
 
            r[0] = xDomainMin;
            r[1] = yDomainMax;
 
            for (unsigned int j= 0; j<height; ++j)
            {
                r[0] = xDomainMin;
                for (unsigned int i= 0; i<width; ++i)
                {
                    fval = function(r);
 
                     // we transform the function value to the corresponding point in [0,maxGreyVal] and round to the next integer. According to the colorscheme (whether fRangeMin or fRangemax is white) we have to choose the greyvalue
                    int greyval = (colorscheme==DEFAULT ? std::floor((fval-fRangeMin)*tFactor+0.5) : maxGreyVal-std::floor((fval-fRangeMin)*tFactor+0.5) );
 
                    // greyval >> endl >> PGM
                    pgmfile << greyval << std::endl;
                    r[0] = xDomainMin + (i+1)*xLength/(width-1);
                }
                r[1] = yDomainMax - (j+1)*yLength/(height-1);
            }
 
            // close PGM (?)
            pgmfile.close();
        }
 
        /* \brief evaluates interpolated greymap at global point
         *
         *  The greymap 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 whicnh to evaluate
         * \returns     the function value
         */
        RangeType operator()(const DomainType& x) const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            DomainType xx = x;
 
            xx[0] = (x[0]-xDomainMin_)/(xDomainMax_-xDomainMin_)*(width_-1);
            xx[1] = (x[1]-yDomainMin_)/(yDomainMax_-yDomainMin_)*(height_-1);
            RangeType r = (*greyMapFunction_)(xx);
 
            return fRangeMin_ + (colorscheme_==DEFAULT ? r/imgParams_[2] : (1-r/imgParams_[2]))*(fRangeMax_-fRangeMin_);
        }
 
        friend auto localFunction(const PortableGreyMap& pgm)
        {
            if (pgm.greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
            auto lf = localFunction(*(pgm.greyMapFunction_));
            return [lf, &pgm](const auto& x) {
                RangeType r = localGreyMapFunction(x);
                return pgm.fRangeMin_ + (pgm.colorscheme_==DEFAULT ? r/pgm.imgParams_[2] : (1-r/pgm.imgParams_[2]))*(pgm.fRangeMax_-pgm.fRangeMin_);
            };
        }
 
        const EntitySet& entitySet() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
            return greyMapFunction_->entitySet();
        }
 
        /* \brief provides read access to the coefficients
         */
        const CoeffType& getCoeffs() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            return coeffs_;
        }
 
        /* \brief provides access to the grid on which the interpolated greymap lives
         */
        const GridType& getGrid() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            return *grid_;
        }
 
        /* \brief provides access to the width of the imported greymap section
         */
        unsigned int getWidth() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            return width_;
        }
 
        /* \brief provides access to the height of the imported greymap section
         */
        unsigned int getHeight() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            return height_;
        }
 
        /* \brief provides access to the width of the full (original) greymap image
         */
        unsigned int getImgWidth() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            return imgParams_[0];
        }
 
        /* \brief provides access to the height of the full (original) greymap image
         */
        unsigned int getImgHeight() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            return imgParams_[1];
        }
 
        /* \brief provides access to the used maximal value in the (original and imported) greymap corresponding to white
         */
        unsigned int getWhiteValue() const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");
 
            return imgParams_[2];
        }
 
        ~PortableGreyMap()
        {
            if (grid_!=0)
                delete grid_;
            if (greyMapFunction_!=0)
                delete greyMapFunction_;
            if (basis_!=0)
                delete basis_;
        }
};
 
#endif