Forked from
agnumpde / dune-fufem
201 commits behind the upstream repository.
-
Carsten Gräser authoredCarsten Gräser authored
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