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