Alamo
PNG.H
Go to the documentation of this file.
1 //
2 // Initialize a field using a PNG image. (2D only)
3 //
4 
5 #ifndef IC_PNG_H
6 #define IC_PNG_H
7 #include <cmath>
8 
9 #ifndef ALAMO_NOPNG
10 #include <stdarg.h>
11 #include <stddef.h>
12 #include <setjmp.h>
13 #include <png.h>
14 #endif
15 
16 #include "IC/IC.H"
17 #include "Util/Util.H"
18 #include "Util/BMP.H"
19 #include "Set/Set.H"
20 #include "IO/ParmParse.H"
21 
22 namespace IC
23 {
24 class PNG : public IC
25 {
26 public:
27  static constexpr const char* name = "png";
28 
29  //enum Type {XYZ, XY, YZ, XZ};
30  enum Fit { Stretch, FitWidth, FitHeight, Coord };
31  enum Channel { R = 0, G = 1, B = 2, A = 3 };
32 
33  PNG(amrex::Vector<amrex::Geometry>& _geom) : IC(_geom) {}
34 
35  PNG(amrex::Vector<amrex::Geometry>& _geom, IO::ParmParse& pp, std::string name) : PNG(_geom)
36  {
37  pp_queryclass(name, *this);
38  }
39 
40  void Define(std::string filename)
41  {
42 #ifndef ALAMO_NOPNG
43 
44  FILE* fp = std::fopen(filename.c_str(), "rb");
45 
46  if (fp == NULL) Util::Abort(INFO, "Cannot find file", filename);
47 
48  png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
49  if (!png) Util::Abort(INFO);
50 
51  png_infop info = png_create_info_struct(png);
52  if (!info) Util::Abort(INFO);
53 
54  if (setjmp(png_jmpbuf(png))) Util::Abort(INFO);
55 
56 
57  png_init_io(png, fp);
58 
59  png_read_info(png, info);
60 
61  png_width = png_get_image_width(png, info);
62  png_height = png_get_image_height(png, info);
63  color_type = png_get_color_type(png, info);
64  bit_depth = png_get_bit_depth(png, info);
65 
66  if (bit_depth == 16)
67  png_set_strip_16(png);
68 
69  if (color_type == PNG_COLOR_TYPE_PALETTE)
70  png_set_palette_to_rgb(png);
71 
72  if (color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8)
73  png_set_expand_gray_1_2_4_to_8(png);
74 
75  if (png_get_valid(png, info, PNG_INFO_tRNS))
76  png_set_tRNS_to_alpha(png);
77 
78  if (color_type == PNG_COLOR_TYPE_RGB ||
79  color_type == PNG_COLOR_TYPE_GRAY ||
80  color_type == PNG_COLOR_TYPE_PALETTE)
81  png_set_filler(png, 0xFF, PNG_FILLER_AFTER);
82 
83  if (color_type == PNG_COLOR_TYPE_GRAY ||
84  color_type == PNG_COLOR_TYPE_GRAY_ALPHA)
85  png_set_gray_to_rgb(png);
86 
87  png_read_update_info(png, info);
88 
89  if (row_pointers) Util::Abort(INFO);
90 
91  row_pointers = (png_bytep*)malloc(sizeof(png_bytep) * png_height);
92  for (int y = 0; y < png_height; y++) {
93  row_pointers[y] = (png_byte*)malloc(png_get_rowbytes(png, info));
94  }
95 
96  png_read_image(png, row_pointers);
97 
98  fclose(fp);
99 
100  png_destroy_read_struct(&png, &info, NULL);
101 #else
102  Util::Abort(INFO,"PNG is disabled");
103 #endif
104  }
105 
106  void Add(const int& lev, Set::Field<Set::Scalar>& a_field, Set::Scalar)
107  {
108 
109 #ifndef ALAMO_NOPNG
110  Set::Vector DX(geom[lev].CellSize());
111  amrex::Box domain = geom[lev].Domain();
112 
113  if (!row_pointers) Util::Abort(INFO, "Running IC without initialization...");
114 
115  Set::Scalar img_width = (Set::Scalar)(png_width - 1);
116  Set::Scalar img_height = (Set::Scalar)(png_height - 1);
117  Set::Scalar img_dx = 1.0;
118  Set::Scalar img_dy = 1.0;
119 
120  amrex::IndexType type = a_field[lev]->ixType();
121  domain.convert(type);
122 
123  Set::Vector domlo(AMREX_D_DECL(geom[lev].ProbLo()[0], geom[lev].ProbLo()[1], 0.0));
124  Set::Vector domhi(AMREX_D_DECL(geom[lev].ProbHi()[0], geom[lev].ProbHi()[1], 0.0));
125 
126  for (amrex::MFIter mfi(*a_field[lev], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
127  {
128  amrex::Box bx;
129  if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
130  if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
131  bx = bx & domain;
132 
133 
134  amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
135  amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
136  {
137  Set::Vector x = Set::Vector::Zero();
138  // NODE
139  if (type == amrex::IndexType::TheNodeType())
140  {
141  x(0) = domlo(0) + ((amrex::Real)(i)) * geom[lev].CellSize()[0];
142  x(1) = domlo(1) + ((amrex::Real)(j)) * geom[lev].CellSize()[1];
143  }
144  else if (type == amrex::IndexType::TheCellType())
145  {
146  x(0) = domlo(0) + ((amrex::Real)(i)+0.5) * geom[lev].CellSize()[0];
147  x(1) = domlo(1) + ((amrex::Real)(j)+0.5) * geom[lev].CellSize()[1];
148  }
149 
150  //Set::Scalar x, y;
151  Set::Vector ximg;
152 
153  if (fit == Fit::Stretch)
154  {
155  ximg(0) = (x(0) - domlo(0)) / (domhi(0) - domlo(0));
156  ximg(1) = (x(1) - domlo(1)) / (domhi(1) - domlo(1));
157  }
158  else if (fit == Fit::FitWidth)
159  {
160  Set::Scalar aspect_ratio = img_width / img_height;
161  ximg(0) = (x(0) - domlo(0)) / (domhi(0) - domlo(0));
162  ximg(1) = (x(1) - domlo(1)) / (domhi(1) - domlo(1));
163  ximg(1) -= 0.5 - 0.5 / aspect_ratio;
164  ximg(1) *= aspect_ratio;
165  }
166  else if (fit == Fit::FitHeight)
167  {
168  Set::Scalar aspect_ratio = img_height / img_width;
169  ximg(0) = (x(0) - domlo(0)) / (domhi(0) - domlo(0));
170  ximg(1) = (x(1) - domlo(1)) / (domhi(1) - domlo(1));
171  ximg(0) -= 0.5 - 0.5 / aspect_ratio;
172  ximg(0) *= aspect_ratio;
173  }
174  else if (fit == Fit::Coord)
175  {
176  ximg(0) = (x(0) - coord_lo(0)) / (coord_hi(0) - coord_lo(0));
177  ximg(1) = (x(1) - coord_lo(1)) / (coord_hi(1) - coord_lo(1));
178  }
179 
180  ximg(0) = std::min(ximg(0), 1.0); ximg(1) = std::min(ximg(1), 1.0);
181  ximg(0) = std::max(ximg(0), 0.0); ximg(1) = std::max(ximg(1), 0.0);
182 
183  ximg(0) *= img_width;
184  ximg(1) *= img_height;
185 
186  int I = (int)(ximg(0));
187  int J = (int)(ximg(1));
188 
189  Set::Scalar x1 = I * img_dx, x2 = (I + 1) * img_dx;
190  Set::Scalar y1 = J * img_dy, y2 = (J + 1) * img_dy;
191 
192 
193  if (I > 0 && I < png_width - 1 &&
194  J>0 && J < png_height - 1)
195  {
196  png_bytep px_sw = &(row_pointers[J][I * 4]);
197  png_bytep px_se = &(row_pointers[J][(I + 1) * 4]);
198  png_bytep px_nw = &(row_pointers[J + 1][I * 4]);
199  png_bytep px_ne = &(row_pointers[J + 1][(I + 1) * 4]);
200 
201 
202  Set::Scalar fQ11 = ((Set::Scalar)(px_sw[channel]) - min) / (max - min);
203  Set::Scalar fQ12 = ((Set::Scalar)(px_nw[channel]) - min) / (max - min);
204  Set::Scalar fQ21 = ((Set::Scalar)(px_se[channel]) - min) / (max - min);
205  Set::Scalar fQ22 = ((Set::Scalar)(px_ne[channel]) - min) / (max - min);
206 
207  field(i, j, k) =
208  (fQ11 * (x2 - ximg(0)) * (y2 - ximg(1)) +
209  fQ21 * (ximg(0) - x1) * (y2 - ximg(1)) +
210  fQ12 * (x2 - ximg(0)) * (ximg(1) - y1) +
211  fQ22 * (ximg(0) - x1) * (ximg(1) - y1)) / (img_dx * img_dy);
212 
213  }
214  else if ((I == 0 || I == png_width - 1) && J < png_height - 1)
215  {
216  png_bytep px_sw = &(row_pointers[J][I * 4]);
217  png_bytep px_nw = &(row_pointers[J + 1][I * 4]);
218 
219  Set::Scalar fQ11 = ((Set::Scalar)(px_sw[channel]) - min) / (max - min);
220  Set::Scalar fQ12 = ((Set::Scalar)(px_nw[channel]) - min) / (max - min);
221  field(i, j, k) = fQ11 + (fQ12 - fQ11) * (ximg(1) - y1);
222  }
223  else if (I < png_width - 1 && (J == 0 || J == png_height - 1))
224  {
225  png_bytep px_sw = &(row_pointers[J][I * 4]);
226  png_bytep px_se = &(row_pointers[J][(I + 1) * 4]);
227 
228  Set::Scalar fQ11 = ((Set::Scalar)(px_sw[channel]) - min) / (max - min);
229  Set::Scalar fQ21 = ((Set::Scalar)(px_se[channel]) - min) / (max - min);
230  field(i, j, k) = fQ11 + (fQ21 - fQ11) * (ximg(0) - x1);
231  }
232  else if (I == png_width - 1 && J == png_height - 1)
233  {
234  png_bytep px_sw = &(row_pointers[J][I * 4]);
235 
236  Set::Scalar fQ11 = ((Set::Scalar)(px_sw[channel]) - min) / (max - min);
237  field(i, j, k) = fQ11;
238  }
239  else
240  {
241  field(i, j, k) = 0.0;
242  }
243 
244  if (field.nComp() > 1) field(i, j, k, 1) = 1.0 - field(i, j, k, 0);
245 
246  });
247  }
248  a_field[lev]->FillBoundary();
249 #else
250  Util::Abort(INFO,"PNG is disabled");
251 #endif
252 
253  };
254 
255 private:
256 #ifndef ALAMO_NOPNG
257  int png_width, png_height;
258  png_byte color_type;
259  png_byte bit_depth;
260  png_bytep* row_pointers = NULL;
261  Set::Vector coord_lo = Set::Vector::Zero();
262  Set::Vector coord_hi = Set::Vector::Zero();
263 
264  //Util::BMP bmp;
265  Fit fit = Fit::Stretch;
266  Channel channel = Channel::G;
267  Set::Scalar min = NAN, max = NAN;
268 #endif
269 
270 public:
271  static void Parse(PNG& value, IO::ParmParse& pp)
272  {
273 #ifndef ALAMO_NOPNG
274  std::string filename;
275  pp_query_file("filename", filename); // BMP filename.
276  value.Define(filename);
277  //value.bmp.Define(filename);
278 
279  std::string fit = "stretch";
280  pp_query("fit", fit); // How to fit. {options: stretch, fitheight, fitwidth}
281  if (fit == "stretch") value.fit = Fit::Stretch;
282  else if (fit == "fitheight") value.fit = Fit::FitHeight;
283  else if (fit == "fitwidth") value.fit = Fit::FitWidth;
284  else if (fit == "coord")
285  {
286  value.fit = Fit::Coord;
287  pp_queryarr("coord.lo", value.coord_lo); // Lower-left coordinates of image in domain
288  pp_queryarr("coord.hi", value.coord_hi); // Upper-right coordinates of image in domain
289  }
290  else Util::Abort(INFO, "Invalid value for png fit - should be stretch/fitheight/fitwidth/coord but received '", fit, "'");
291 
292  std::string channel = "g";
293  //pp_query("channel", channel); // Color channel to use (options: r, R, g, G, b, B, a, A)
294 
295  pp_query_validate("channel", channel, {"r","g","b","a","R","G","B","A"}); // Color channel to use (options: r, R, g, G, b, B, a, A)
296  if (channel == "r" || channel == "R") value.channel = Channel::R;
297  else if (channel == "g" || channel == "G") value.channel = Channel::G;
298  else if (channel == "b" || channel == "B") value.channel = Channel::B;
299  else if (channel == "a" || channel == "A") value.channel = Channel::A;
300  else Util::Abort(INFO, "Invalid value for bmp channel - should be r/g/b/a but received '", channel, "'");
301 
302  value.min = 0.0; //(Set::Scalar) value.bmp.min()[value.channel];
303  value.max = 255.0; //(Set::Scalar) value.bmp.max()[value.channel];
304  pp_query_default("min", value.min, 0.0 ); // Scaling value - minimum
305  pp_query_default("max", value.max, 255.0); // Scaling value - maximum
306 #else
307  Util::Abort(INFO,"PNG is disabled");
308 #endif
309  }
310 };
311 }
312 #endif
Util::filename
std::string filename
Definition: Util.cpp:19
IC::PNG::row_pointers
png_bytep * row_pointers
Definition: PNG.H:260
IC::IC::geom
amrex::Vector< amrex::Geometry > & geom
Definition: IC.H:54
IC::PNG::Parse
static void Parse(PNG &value, IO::ParmParse &pp)
Definition: PNG.H:271
IC::PNG::Stretch
@ Stretch
Definition: PNG.H:30
IC::PNG::G
@ G
Definition: PNG.H:31
IC::PNG::fit
Fit fit
Definition: PNG.H:265
IC::PNG::png_height
int png_height
Definition: PNG.H:257
BC::AMREX_D_DECL
@ AMREX_D_DECL
Definition: BC.H:34
IC::PNG::min
Set::Scalar min
Definition: PNG.H:267
Util.H
IC::PNG::png_width
int png_width
Definition: PNG.H:253
IC::PNG::Channel
Channel
Definition: PNG.H:31
IC::PNG::Coord
@ Coord
Definition: PNG.H:30
Set::Field< Set::Scalar >
Definition: Set.H:236
IC::PNG::PNG
PNG(amrex::Vector< amrex::Geometry > &_geom, IO::ParmParse &pp, std::string name)
Definition: PNG.H:35
IC::PNG::Fit
Fit
Definition: PNG.H:30
IC::PNG::coord_lo
Set::Vector coord_lo
Definition: PNG.H:261
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition: Base.H:20
ParmParse.H
pp_query
#define pp_query(...)
Definition: ParmParse.H:105
IC::PNG::FitWidth
@ FitWidth
Definition: PNG.H:30
Set::Scalar
amrex::Real Scalar
Definition: Base.H:19
pp_queryclass
#define pp_queryclass(...)
Definition: ParmParse.H:106
IC::PNG
Definition: PNG.H:24
IC::PNG::Define
void Define(std::string filename)
Definition: PNG.H:40
IC::PNG::bit_depth
png_byte bit_depth
Definition: PNG.H:259
IC::PNG::B
@ B
Definition: PNG.H:31
IC::PNG::FitHeight
@ FitHeight
Definition: PNG.H:30
pp_query_file
#define pp_query_file(...)
Definition: ParmParse.H:102
IC::PNG::channel
Channel channel
Definition: PNG.H:266
pp_query_default
#define pp_query_default(...)
Definition: ParmParse.H:100
Util::Abort
void Abort(const char *msg)
Definition: Util.cpp:166
Set.H
IC
Definition: BMP.H:18
pp_query_validate
#define pp_query_validate(...)
Definition: ParmParse.H:101
IO::ParmParse
Definition: ParmParse.H:112
IC::PNG::PNG
PNG(amrex::Vector< amrex::Geometry > &_geom)
Definition: PNG.H:33
BMP.H
IC::PNG::name
static constexpr const char * name
Definition: PNG.H:27
IC::PNG::R
@ R
Definition: PNG.H:31
INFO
#define INFO
Definition: Util.H:20
IC.H
IC::PNG::Add
void Add(const int &lev, Set::Field< Set::Scalar > &a_field, Set::Scalar)
Definition: PNG.H:106
IC::PNG::A
@ A
Definition: PNG.H:31
IC::PNG::color_type
png_byte color_type
Definition: PNG.H:258
IC::PNG::coord_hi
Set::Vector coord_hi
Definition: PNG.H:262
IC::PNG::max
Set::Scalar max
Definition: PNG.H:267
pp_queryarr
#define pp_queryarr(...)
Definition: ParmParse.H:103
Generated on Sat Feb 8 2025 00:07:39 for Alamo by   doxygen 1.8.17