Skip to content

Commit

Permalink
Merge pull request #485 from JulienDoerner/GridProps
Browse files Browse the repository at this point in the history 
Loading Grids with properties
  • Loading branch information
@lukasmerten
lukasmerten authored Jun 27, 2024
2 parents dc3299b + a1adcde commit fe44552
Show file tree
Hide file tree
Showing 5 changed files with 368 additions and 74 deletions.
113 changes: 47 additions & 66 deletions doc/pages/example_notebooks/density/density_grid_sampling.ipynb
View file

Large diffs are not rendered by default.

31 changes: 31 additions & 0 deletions include/crpropa/Grid.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@

#include <vector>
#include <type_traits>
#include <string>
#include <sstream>
#if HAVE_SIMD
#include <immintrin.h>
#include <smmintrin.h>
Expand Down Expand Up @@ -128,6 +130,20 @@ class GridProperties: public Referenced {
void setClipVolume(bool b) {
clipVolume = b;
}

/** show all GridProperty parameters
* @param unit unit for the lengthscale (origin, spacing). Default is 1 = SI units
*/
std::string getDescription(double unit = 1) const {
std::stringstream ss;
ss << "GridProperties:\torigin: " << origin / unit
<< "\t" << "Nx: " << Nx << " Ny: " << Ny << " Nz: " << Nz
<< "\t" << "spacing: " << spacing / unit
<< "\t" << "refletive: " << reflective
<< "\t" << "interpolation: " << ipol
<< "\n";
return ss.str();
}
};

/**
Expand Down Expand Up @@ -247,6 +263,21 @@ class Grid: public Referenced {
}
}

interpolationType getInterpolationType() {
return ipolType;
}

std::string getInterpolationTypeName() {
if (ipolType == TRILINEAR)
return "TRILINEAR";
if (ipolType == TRICUBIC)
return "TRICUBIC";
if (ipolType == NEAREST_NEIGHBOUR)
return "NEAREST_NEIGHBOUR";

return "NOT_UNDERSTOOD";
}

/** returns the position of the lower left front corner of the volume */
Vector3d getOrigin() const {
return origin;
Expand Down
18 changes: 16 additions & 2 deletions include/crpropa/GridTools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,12 @@ void dumpGrid(ref_ptr<Grid3f> grid, std::string filename,
void dumpGrid(ref_ptr<Grid1f> grid, std::string filename,
double conversion = 1);

/** Load a Grid3f from a plain text file based on the gridproperties stored in the header
@param filename name of the input file
@param conversion multiply every point in a grid by a conversion factor
*/
ref_ptr<Grid3f> loadGrid3fFromTxt(std::string filename, double conversion = 1);

/** Load a Grid3f grid from a plain text file.
@param grid a vector grid (Grid3f) to which the points will be loaded
@param filename name of input file
Expand All @@ -125,6 +131,12 @@ void dumpGrid(ref_ptr<Grid1f> grid, std::string filename,
void loadGridFromTxt(ref_ptr<Grid3f> grid, std::string filename,
double conversion = 1);

/** Load a Grid1f from a plain text file based on the gridproperties stored in the header
@param filename name of the input file
@param conversion multiply every point in a grid by a conversion factor
*/
ref_ptr<Grid1f> loadGrid1fFromTxt(std::string filename, double conversion = 1);

/** Load a Grid1f from a plain text file
@param grid a scalar grid (Grid1f) to which the points will be loaded
@param filename name of input file
Expand All @@ -137,17 +149,19 @@ void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename,
@param grid a vector grid (Grid3f)
@param filename name of output file
@param conversion multiply every point in grid by a conversion factor
@param storeProperties if true the grid properties are stored as a comment
*/
void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename,
double conversion = 1);
double conversion = 1, bool storeProperties = false);

/** Dump a Grid1f to a plain text file.
@param grid a scalar grid (Grid1f)
@param filename name of output file
@param conversion multiply every point in grid by a conversion factor
@param storeProperties if true the grid properties are stored as a comment
*/
void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename,
double conversion = 1);
double conversion = 1, bool storeProperties = false);

#ifdef CRPROPA_HAVE_FFTW3F
/**
Expand Down
166 changes: 163 additions & 3 deletions src/GridTools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,15 +260,84 @@ void loadGridFromTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
fin.close();
}

ref_ptr<Grid3f> loadGrid3fFromTxt(std::string filename, double c) {
std::ifstream fin(filename.c_str());
if (!fin) {
std::stringstream ss;
ss << "load Grid3f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// search in header lines for GridProperties
while (fin.peek() == '#') {
std::string line;
std::getline(fin, line);

// find gridproperties in the header
if (line.rfind("GridProperties:") == 2) {
GridProperties gp(Vector3d(0.), 1, 1, 1, 1.); // simple grid properties for default
std::stringstream ss(line);

// skip first names and check type
std::string name, type;
ss >> name >> name >> name >> type;
if (type != "Grid3f")
throw std::runtime_error("Tried to load Grid3f, but Gridproperties assume grid type " + type);

// grid origin
double x, y, z;
ss >> name >> x >> y >> z ;
gp.origin = Vector3d(x, y, z);

// grid size
ss >> name >> gp.Nx >> gp.Ny >> gp.Nz;

// spacing
double dX, dY, dZ;
ss >> name >> dX >> dY >> dZ;
gp.spacing = Vector3d(dX, dY, dZ);

// reflective
ss >> name >> gp.reflective;

// clip volume
bool clip;
ss >> name >> clip;
gp.setClipVolume(clip);

// interpolation type
ss >> name >> type;
if (type == "TRICUBIC")
gp.setInterpolationType(TRICUBIC);
else if (type == "NEAREST_NEIGHBOUR")
gp.setInterpolationType(NEAREST_NEIGHBOUR);
else
gp.setInterpolationType(TRILINEAR);

// create new grid
ref_ptr<Grid3f> grid = new Grid3f(gp);
fin.close();

// load data for grid
loadGridFromTxt(grid, filename, c);

return grid;
}
}
throw std::runtime_error("could not find GridProperties in file " + filename);
}


void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
std::ifstream fin(filename.c_str());
if (!fin) {
std::stringstream ss;
ss << "load Grid1f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// skip header lines
while (fin.peek() == '#')
while (fin.peek() == '#')
fin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');

for (int ix = 0; ix < grid->getNx(); ix++) {
Expand All @@ -285,13 +354,92 @@ void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
fin.close();
}

void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
ref_ptr<Grid1f> loadGrid1fFromTxt(std::string filename, double c) {
std::ifstream fin(filename.c_str());
if (!fin) {
std::stringstream ss;
ss << "load Grid1f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// search in header lines for GridProperties
while (fin.peek() == '#') {
std::string line;
std::getline(fin, line);

// find gridproperties in the header
if (line.rfind("GridProperties:") == 2) {
GridProperties gp(Vector3d(0.), 1, 1, 1, 1.); // simple grid properties for default
std::stringstream ss(line);

// skip first names and check type
std::string name, type;
ss >> name >> name >> name >> type;
if (type != "Grid1f")
throw std::runtime_error("Tried to load Grid1f, but Gridproperties assume grid type " + type);

// grid origin
double x, y, z;
ss >> name >> x >> y >> z ;
gp.origin = Vector3d(x, y, z);

// grid size
ss >> name >> gp.Nx >> gp.Ny >> gp.Nz;

// spacing
double dX, dY, dZ;
ss >> name >> dX >> dY >> dZ;
gp.spacing = Vector3d(dX, dY, dZ);

// reflective
ss >> name >> gp.reflective;

// clip volume
bool clip;
ss >> name >> clip;
gp.setClipVolume(clip);

// interpolation type
ss >> name >> type;
if (type == "TRICUBIC")
gp.setInterpolationType(TRICUBIC);
else if (type == "NEAREST_NEIGHBOUR")
gp.setInterpolationType(NEAREST_NEIGHBOUR);
else
gp.setInterpolationType(TRILINEAR);

// create new grid
ref_ptr<Grid1f> grid = new Grid1f(gp);
fin.close();

// load data for grid
loadGridFromTxt(grid, filename, c);

return grid;
}
}
throw std::runtime_error("could not find GridProperties in file " + filename);
}

void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c, bool saveProp) {
std::ofstream fout(filename.c_str());
if (!fout) {
std::stringstream ss;
ss << "dump Grid3f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// store the properties in the file as header information
if (saveProp) {
fout << "# GridProperties: Type Grid3f"
<< "\t" << "origin: " << grid -> getOrigin()
<< "\t" << "gridsize: " << grid -> getNx() << " " << grid -> getNy() << " " << grid -> getNz()
<< "\t" << "spacing: " << grid -> getSpacing ()
<< "\t" << "reflective: " << grid -> isReflective()
<< "\t" << "clipVolume: " << grid -> getClipVolume()
<< "\t" << "interpolation: " << grid -> getInterpolationTypeName() << "\n";
}

for (int ix = 0; ix < grid->getNx(); ix++) {
for (int iy = 0; iy < grid->getNy(); iy++) {
for (int iz = 0; iz < grid->getNz(); iz++) {
Expand All @@ -303,13 +451,25 @@ void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
fout.close();
}

void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename, double c, bool saveProp) {
std::ofstream fout(filename.c_str());
if (!fout) {
std::stringstream ss;
ss << "dump Grid1f: " << filename << " not found";
throw std::runtime_error(ss.str());
}

// save properties as header information
if (saveProp) {
fout << "# GridProperties: Type Grid1f"
<< "\t" << "origin: " << grid -> getOrigin()
<< "\t" << "gridsize: " << grid -> getNx() << " " << grid -> getNy() << " " << grid -> getNz()
<< "\t" << "spacing: " << grid -> getSpacing ()
<< "\t" << "reflective: " << grid -> isReflective()
<< "\t" << "clipVolume: " << grid -> getClipVolume()
<< "\t" << "interpolation: " << grid -> getInterpolationTypeName() << "\n";
}

for (int ix = 0; ix < grid->getNx(); ix++) {
for (int iy = 0; iy < grid->getNy(); iy++) {
for (int iz = 0; iz < grid->getNz(); iz++) {
Expand Down
Loading

0 comments on commit fe44552

Please sign in to comment.