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