LCOV - code coverage report
Current view: top level - src/IC - BMP.H (source / functions) Coverage Total Hit
Test: coverage_merged.info Lines: 65.8 % 117 77
Test Date: 2025-07-23 22:41:20 Functions: 83.3 % 6 5

            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
        

Generated by: LCOV version 2.0-1