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