Line data Source code
1 : //
2 : // Initialize a field using a bitmap image. (2D only)
3 : //
4 : // Note that in GIMP, you must select "do not write color space information"
5 : // and "24 bit R8 G8 B8" when exporting the BMP file.
6 : //
7 :
8 : #ifndef IC_BMP_H
9 : #define IC_BMP_H
10 : #include <cmath>
11 :
12 : #include "IC/IC.H"
13 : #include "Util/Util.H"
14 : #include "Util/BMP.H"
15 : #include "Unit/Unit.H"
16 : #include "Set/Set.H"
17 : #include "IO/ParmParse.H"
18 :
19 : namespace IC
20 : {
21 : class BMP : public IC<Set::Scalar>
22 : {
23 : public:
24 : static constexpr const char* name = "bmp";
25 :
26 : //enum Type {XYZ, XY, YZ, XZ};
27 : enum Fit {Stretch,FitWidth,FitHeight,Coord};
28 : enum Channel {R=0, G=1, B=2};
29 :
30 2 : BMP (amrex::Vector<amrex::Geometry> &_geom) : IC(_geom) {}
31 :
32 2 : BMP (amrex::Vector<amrex::Geometry> &_geom, IO::ParmParse &pp, std::string name) : BMP(_geom)
33 : {
34 2 : pp_queryclass(name,*this);
35 2 : }
36 0 : BMP (amrex::Vector<amrex::Geometry> &_geom, Unit a_unit, IO::ParmParse &pp, std::string name) : BMP(_geom)
37 : {
38 0 : unit = a_unit;
39 0 : pp_queryclass(name,*this);
40 0 : }
41 :
42 : void Define(std::string bmpfilename)
43 : {
44 : bmp.Define(bmpfilename);//"Interface_Blur2.bmp");
45 : }
46 :
47 12 : void Add(const int &lev, Set::Field<Set::Scalar> &a_field, Set::Scalar)
48 : {
49 :
50 12 : Set::Vector DX(geom[lev].CellSize());
51 12 : amrex::Box domain = geom[lev].Domain();
52 :
53 : //Set::Scalar width = geom[lev].ProbHi()[0] - geom[lev].ProbLo()[0];
54 : //Set::Scalar height = geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1];
55 :
56 12 : Set::Scalar img_width = (Set::Scalar)(bmp.nx-1);
57 12 : Set::Scalar img_height = (Set::Scalar)(bmp.ny-1);
58 12 : Set::Scalar img_dx = 1.0;
59 12 : Set::Scalar img_dy = 1.0;
60 :
61 12 : amrex::IndexType type = a_field[lev]->ixType();
62 12 : domain.convert(type);
63 :
64 12 : Set::Vector domlo(AMREX_D_DECL(geom[lev].ProbLo()[0],geom[lev].ProbLo()[1],0.0));
65 12 : Set::Vector domhi(AMREX_D_DECL(geom[lev].ProbHi()[0],geom[lev].ProbHi()[1],0.0));
66 :
67 24 : for (amrex::MFIter mfi(*a_field[lev],amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
68 : {
69 12 : amrex::Box bx;
70 24 : if (type == amrex::IndexType::TheNodeType()) bx = mfi.grownnodaltilebox();
71 24 : if (type == amrex::IndexType::TheCellType()) bx = mfi.growntilebox();
72 12 : bx = bx & domain;
73 :
74 :
75 12 : amrex::Array4<Set::Scalar> const& field = a_field[lev]->array(mfi);
76 12 : amrex::ParallelFor (bx,[=] AMREX_GPU_DEVICE(int i, int j, int k)
77 : {
78 241664 : Set::Vector x = Set::Vector::Zero();
79 : // NODE
80 483328 : if (type == amrex::IndexType::TheNodeType())
81 : {
82 0 : x(0) = domlo(0) + ((amrex::Real)(i)) * geom[lev].CellSize()[0];
83 0 : x(1) = domlo(1) + ((amrex::Real)(j)) * geom[lev].CellSize()[1];
84 : }
85 483328 : else if (type == amrex::IndexType::TheCellType())
86 : {
87 241664 : x(0) = domlo(0) + ((amrex::Real)(i) + 0.5) * geom[lev].CellSize()[0];
88 241664 : x(1) = domlo(1) + ((amrex::Real)(j) + 0.5) * geom[lev].CellSize()[1];
89 : }
90 :
91 : //Set::Scalar x, y;
92 241664 : Set::Vector ximg;
93 :
94 241664 : if (fit == Fit::Stretch)
95 : {
96 0 : ximg(0) = (x(0) - domlo(0)) / (domhi(0) - domlo(0));
97 0 : ximg(1) = (x(1) - domlo(1)) / (domhi(1) - domlo(1));
98 : }
99 241664 : else if (fit == Fit::FitWidth)
100 : {
101 241664 : Set::Scalar aspect_ratio = img_width / img_height;
102 241664 : ximg(0) = (x(0) - domlo(0)) / (domhi(0) - domlo(0));
103 241664 : ximg(1) = (x(1) - domlo(1)) / (domhi(1) - domlo(1));
104 241664 : ximg(1) -= 0.5 - 0.5 / aspect_ratio;
105 241664 : ximg(1) *= aspect_ratio;
106 : }
107 0 : else if (fit == Fit::FitHeight)
108 : {
109 0 : Set::Scalar aspect_ratio = img_height / img_width;
110 0 : ximg(0) = (x(0) - domlo(0)) / (domhi(0) - domlo(0));
111 0 : ximg(1) = (x(1) - domlo(1)) / (domhi(1) - domlo(1));
112 0 : ximg(0) -= 0.5 - 0.5 / aspect_ratio;
113 0 : ximg(0) *= aspect_ratio;
114 : }
115 0 : else if (fit == Fit::Coord)
116 : {
117 0 : ximg(0) = (x(0) - coord_lo(0)) / (coord_hi(0)-coord_lo(0));
118 0 : ximg(1) = (x(1) - coord_lo(1)) / (coord_hi(1)-coord_lo(1));
119 : }
120 :
121 241664 : ximg(0) = std::min(ximg(0),1.0); ximg(1) = std::min(ximg(1),1.0);
122 241664 : ximg(0) = std::max(ximg(0),0.0); ximg(1) = std::max(ximg(1),0.0);
123 :
124 241664 : ximg(0) *= img_width;
125 241664 : ximg(1) *= img_height;
126 :
127 241664 : int I = (int)(ximg(0));
128 241664 : int J = (int)(ximg(1));
129 :
130 241664 : Set::Scalar x1 = I*img_dx, x2 = (I+1)*img_dx;
131 241664 : Set::Scalar y1 = J*img_dy, y2 = (J+1)*img_dy;
132 :
133 241664 : if (I < bmp.nx-1 && J < bmp.ny - 1)
134 : {
135 483328 : Set::Scalar fQ11 = ((Set::Scalar)(bmp(I, J )[channel]) - min) / (max - min);
136 483328 : Set::Scalar fQ12 = ((Set::Scalar)(bmp(I, J+1)[channel]) - min) / (max - min);
137 483328 : Set::Scalar fQ21 = ((Set::Scalar)(bmp(I+1,J )[channel]) - min) / (max - min);
138 483328 : Set::Scalar fQ22 = ((Set::Scalar)(bmp(I+1,J+1)[channel]) - min) / (max - min);
139 :
140 483328 : field(i,j,k) =
141 241664 : (fQ11*(x2-ximg(0))*(y2-ximg(1)) + fQ21*(ximg(0)-x1)*(y2-ximg(1)) + fQ12*(x2-ximg(0))*(ximg(1)-y1) + fQ22*(ximg(0)-x1)*(ximg(1)-y1)) / (img_dx * img_dy);
142 241664 : }
143 0 : else if (I == bmp.nx-1 && J < bmp.ny - 1)
144 : {
145 0 : Set::Scalar fQ11 = ((Set::Scalar)(bmp(I ,J )[channel]) - min) / (max - min);
146 0 : Set::Scalar fQ12 = ((Set::Scalar)(bmp(I ,J+1)[channel]) - min) / (max - min);
147 0 : field(i,j,k) = fQ11 + (fQ12-fQ11) * (ximg(1) - y1);
148 0 : }
149 0 : else if (I < bmp.nx-1 && J == bmp.ny - 1)
150 : {
151 0 : Set::Scalar fQ11 = ((Set::Scalar)(bmp(I ,J )[channel]) - min) / (max - min);
152 0 : Set::Scalar fQ21 = ((Set::Scalar)(bmp(I+1,J )[channel]) - min) / (max - min);
153 0 : field(i,j,k) = fQ11 + (fQ21-fQ11) * (ximg(0) - x1);
154 0 : }
155 0 : else if (I == bmp.nx-1 && J == bmp.ny - 1)
156 : {
157 0 : Set::Scalar fQ11 = ((Set::Scalar)(bmp(I ,J )[channel]) - min) / (max - min);
158 0 : field(i,j,k) = fQ11;
159 0 : }
160 : else
161 : {
162 0 : field(i,j,k) = 0.0;
163 : }
164 :
165 483328 : if (field.nComp() > 1) field(i,j,k,1) = 1.0 - field(i,j,k,0);
166 :
167 241664 : field(i,j,k,0) *= scaling_factor;
168 483328 : if (field.nComp() > 1)
169 0 : field(i,j,k,1) *= scaling_factor;
170 :
171 241664 : });
172 12 : }
173 12 : a_field[lev]->FillBoundary();
174 12 : };
175 :
176 : private:
177 : Util::BMP bmp;
178 : Fit fit = Fit::Stretch;
179 : Channel channel = Channel::G;
180 : Set::Scalar min=NAN, max=NAN;
181 : Set::Vector coord_lo = Set::Vector::Zero();
182 : Set::Vector coord_hi = Set::Vector::Zero();
183 : Unit unit = Unit::Less();
184 : Set::Scalar scaling_factor = NAN;
185 :
186 : public:
187 2 : static void Parse(BMP & value, IO::ParmParse & pp)
188 : {
189 2 : std::string filename;
190 2 : pp_query_file("filename",filename); // BMP filename.
191 2 : value.bmp.Define(filename);
192 :
193 2 : std::string fit;
194 : // How to position image in space
195 8 : pp_query_validate("fit",fit,{"stretch","fitheight","fitwidth","coord"});
196 2 : if (fit=="stretch") value.fit = Fit::Stretch;
197 2 : else if (fit=="fitheight") value.fit = Fit::FitHeight;
198 2 : else if (fit=="fitwidth") value.fit = Fit::FitWidth;
199 0 : else if (fit=="coord")
200 : {
201 0 : value.fit = Fit::Coord;
202 0 : pp_queryarr("coord.lo",value.coord_lo); // Location of lower-left corner in the domain
203 0 : pp_queryarr("coord.hi",value.coord_hi); // Location of upper-right corner in the domain
204 : }
205 0 : else Util::Abort(INFO,"Invalid value for bmp fit - should be stretch/fitheight/fitwidth but received '",fit,"'");
206 :
207 2 : std::string channel;
208 : // Color channel to use
209 8 : pp_query_validate("channel",channel,{"r","g","b","R","G","B"});
210 2 : if (channel=="r" || channel=="R") value.channel = Channel::R;
211 2 : else if (channel=="g" || channel=="G") value.channel = Channel::G;
212 0 : else if (channel=="b" || channel=="B") value.channel = Channel::B;
213 0 : else Util::Abort(INFO,"Invalid value for bmp channel - should be r/g/b but received '",channel,"'");
214 :
215 2 : value.min = (Set::Scalar) value.bmp.min()[value.channel];
216 2 : value.max = (Set::Scalar) value.bmp.max()[value.channel];
217 4 : pp_query_default("min",value.min,0.0); // Scaling value - minimum
218 4 : pp_query_default("max",value.max,255.0); // Scaling value - maximum
219 :
220 2 : value.scaling_factor = value.unit.normalized_value();
221 2 : }
222 : };
223 : }
224 : #endif
|
