Skip to content
Snippets Groups Projects
Forked from agnumpde / dune-fufem
201 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
portablegreymap.hh 20.03 KiB
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
#ifndef PORTABLEGREYMAP_HH
#define PORTABLEGREYMAP_HH

#include <fstream>
#include <iostream>
#include <string>
#include <vector>
#include <array>


#include <dune/grid/yaspgrid.hh>

#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functionspacebases/q1nodalbasis.hh>

/** \brief  Class wrapping a greymap (acquired from PGM ASCII (P2)-files) providing it as a piecewise bi-linear GridFunction 
 *
 *  A rectangular domain \f$ \Omega=[xDomainMin,xDomainMax]\times[yDomainMin,yDomainMax] \f$ is assumed.
 *
 *  \todo \li add support for binary read/write
 *  \li unify PPM,PGM -> PAM
 */
class PortableGreyMap:
    virtual public VirtualDifferentiableFunction<Dune::YaspGrid<2>::Codim<0>::Geometry::GlobalCoordinate, Dune::FieldVector<double, 1> >
{
    public:
        typedef Dune::YaspGrid<2> GridType;

        typedef VirtualDifferentiableFunction<GridType::Codim<0>::Geometry::GlobalCoordinate, Dune::FieldVector<double, 1> > BaseType;
        typedef GridType::Codim<0>::Geometry::LocalCoordinate LocalDomainType;
        typedef GridType::Codim<0>::Entity Element;

        typedef BaseType::DomainType DomainType;
        typedef BaseType::RangeType RangeType;
        typedef BaseType::DerivativeType DerivativeType;

        typedef Q1NodalBasis<GridType::LeafGridView> Basis;

        typedef Dune::BlockVector<RangeType> CoeffType;
        typedef BasisGridFunction<Basis, CoeffType> FunctionType;

        /** An enum identifying the employed color scheme.
         */
        enum ColorScheme {
                DEFAULT,    /**< The default scheme is BLACK <=> fRangeMin, WHITE <=> fRangeMax */
                REVERSE     /**< The reverse scheme is BLACK <=> fRangeMax, WHITE <=> fRangeMin */
        };

    private:
//        typedef GridType::Traits::Codim<0>::Entity Entity;

        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_;

    public:
        /** \brief Constructor
         *
         * \param xDomainMin  lower x-interval boundary of computational domain
         * \param xDomainMax  upper x-interval boundary of computational domain
         * \param yDomainMin  lower y-interval boundary of computational domain
         * \param yDomainMax  upper y-interval boundary of computational domain
         * \param fRangeMin   lowest function value that shall be mapped to a grey value; DEFAULT = 0
         * \param fRangeMax   largest function value that shall be mapped to a grey value; DEFAULT = 1
         * \param colorscheme the employed PortableGreyMap::ColorScheme; DEFAULT = DEFAULT
         */
        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.");
        }

        /** \brief Constructor
         *
         * \param fRangeMin   lowest function value that shall be mapped to a grey value; DEFAULT = 0
         * \param fRangeMax   largest function value that shall be mapped to a grey value; DEFAULT = 1
         * \param colorscheme the employed PortableGreyMap::ColorScheme; DEFAULT = DEFAULT
         */
        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.");
        }

        /** \brief reads the specified rectangular section of a greymap (PGM-ASCII)
         *
         * Method creates a YaspGrid<2> (structured rectangular grid) with width times height nodes. The greyvalues of the corresponding pixels are taken to be the
         * nodal values of a Q1Function on this grid. Thus we get a pixel centered linear interpolation of the greymap section.
         *
         * \param filename    the name of the PGM file to be read
         * \param width       the width (pixels) of the section to be imported (for negative values the width of the imported image is assumed); DEFAULT=-1
         * \param height      the height (pixels) of the section to be imported (for negative values the height of the imported image is assumed); DEFAULT=-1
         * \param xoffset     the x-position (in pixels) of the upper left corner of the importsection in the greymap counted from upper left; DEFAULT=0
         * \param yoffset     the y-position (in pixels) of the upper left corner of the importsection in the greymap counted from upper left; DEFAULT=0
         */
        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(*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.");

        }

        /** \brief creates a PGM (P2) file from a given scalar function
         *
         * Method creates a greymap (PGM P2) with width times height pixels whose grey values
         * correspond to the given function's domain, range, and values.
         *
         * \param filename      the name of the PGM file to be read
         * \param function      the function to be converted to a greymap
         * \param fRangeMin     lower range interval boundary (function value that will correspond to black (white)); lower function values will be "cut off"
         * \param fRangeMax     upper range interval boundary (function value that will correspond to white (black)); larger function values will be "cut off"
         * \param xDomainMin    lower x-interval boundary of function's domain (or rather the part of it to be exported as GreyMap)
         * \param xDomainMax    upper x-interval boundary of function's domain (or rather the part of it to be exported as GreyMap)
         * \param yDomainMin    lower y-interval boundary of function's domain (or rather the part of it to be exported as GreyMap)
         * \param yDomainMax    upper y-interval boundary of function's domain (or rather the part of it to be exported as GreyMap)
         * \param width         the horizontal resolution (absolute number of pixels) of the greymap to be created
         * \param height        the vertical resolution (absolute number of pixels) of the greymap to be created
         * \param info          some information concerning the function; DEFAULT=""
         * \param maxGreyVal    the value that shall correspond to white in the PGM (0 is black), has to satisfy 0 < maxGreyVal < 65536; determines the number of greyshades used for the greymap; DEFAULT = 255
         * \param colorscheme the employed PortableGreyMap::ColorScheme; DEFAULT = DEFAULT
         */
        template <typename FType>
        static void exportGreyMap(  const char* filename,
                                    const FType& function,
                                    typename FType::RangeType::field_type fRangeMin, typename FType::RangeType::field_type fRangeMax,
                                    const typename FType::DomainType::field_type xDomainMin, const typename FType::DomainType::field_type xDomainMax,
                                    const typename FType::DomainType::field_type yDomainMin, const typename FType::DomainType::field_type yDomainMax,
                                    const unsigned int width, const unsigned int height, 
                                    const std::string info="",
                                    const unsigned int maxGreyVal=255, const ColorScheme colorscheme=DEFAULT)
        {
            /* We could elegantly evade the following check by not making the function type the template parameter but solely the domain and range (field) types, hardwiring
             * DimRange and DimDomain and have the method take a Dune::VirtualFunction reference. However I find it to be more userfriendly as is for the following reason:
             * Having the function type as template parameter it can be inferred from the argument and thus does not have to be explicitely given whereas one would
             * have to specify both the DomainType::field_type and RangeType::field_type explicitely each time.
             * Another way would be to hardwire the FunctionType altogether, but then one would be restricted to e.g. double as field_type in (Range|Domain)Type which seems unnecessarily restrictive
             */
            if (FType::DomainType::dimension != 2 or FType::RangeType::dimension != 1)
                DUNE_THROW(Dune::Exception, "Given function is not a scalar function on a 2D domain (thrown by PortableGreyMap::exportGreyMap).");

            typedef typename FType::DomainType::field_type DFT;
            typedef typename FType::RangeType::field_type  RFT;

            DFT xLength = xDomainMax - xDomainMin,
                yLength = yDomainMax - yDomainMin;

            RFT 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;

            typename FType::DomainType r(0.0);
            typename FType::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)
                {
//                    std::cout << r << std::endl;
                    function.evaluate(r,fval);

                     // 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
         * \param y     will contain the function value
         */
        void evaluate(const DomainType& x, RangeType& y) 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 = 0.0;
            greyMapFunction_->evaluate(xx,r);

            y = fRangeMin_ + (colorscheme_==DEFAULT ? r/imgParams_[2] : (1-r/imgParams_[2]))*(fRangeMax_-fRangeMin_);
            return;
        }

        /* \brief evaluates interpolated greymap at local point
         *
         *  The greymap is evaluated at local coordinates x in the given element.
         *
         * \param ntt   element in which to evaluate
         * \param x     point in \f$ \Omega \f$ at which to evaluate
         * \param y     will contain the function value
         */
        void evaluateLocal(const Element& ntt, const LocalDomainType& x, RangeType& y) const
        {
            if (greyMapFunction_==0)
                DUNE_THROW(Dune::Exception,"No greymap has been loaded. Use readGreyMap() first.");

            RangeType r = 0.0;
            greyMapFunction_->evaluateLocal(ntt,x,r);

            y = fRangeMin_ + (colorscheme_==DEFAULT ? r/imgParams_[2] : (1-r/imgParams_[2]))*(fRangeMax_-fRangeMin_);
            return;
        }

        /** \brief evaluation of 1st derivative in global coordinates
         *
         * Default implementation using hierarchic search and evaluateDerivativeLocal
         *
         * \warning may be VERY SLOW due to hierarchical search
         *
         *  \param x point (global coordinates) at which to evaluate the derivative
         *  \param d will contain the derivative at x after return
         */
        virtual void evaluateDerivative(const DomainType& x DUNE_UNUSED, DerivativeType& d DUNE_UNUSED) const
        {
          DUNE_THROW(Dune::NotImplemented, "Derivative not implemented");
        }

        /** \brief evaluation of derivative in local coordinates
         *
         *  \param e Evaluate in local coordinates with respect to this element.
         *  \param x point in local coordinates at which to evaluate the derivative
         *  \param d will contain the derivative at x after return
         */
        virtual void evaluateDerivativeLocal(const Element& e DUNE_UNUSED, const LocalDomainType& x DUNE_UNUSED, DerivativeType& d DUNE_UNUSED) const
        {
          DUNE_THROW(Dune::NotImplemented, "Derivative not implemented");
        }

        /* \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