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