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