C++ Surface::GetNJ方法代码示例

void
GridMapping::calculateSurfaceFromVelocity(NRLib::Grid<float> * velocity,
                                          const Simbox       * timeSimbox)
{
  if(z0Grid_==NULL || z1Grid_==NULL)
  {
    Surface * isochore;
    if(z0Grid_==NULL)
      isochore = new Surface(*z1Grid_);
    else
      isochore = new Surface(*z0Grid_);

    //
    // If a constant time surface has been used, it may have only four grid
    // nodes. To handle this situation we use the grid resolution whenever
    // this is larger than the surface resolution.
    //
    int maxNx = std::max(timeSimbox->getnx(), static_cast<int>(isochore->GetNI()));
    int maxNy = std::max(timeSimbox->getny(), static_cast<int>(isochore->GetNJ()));
    isochore->Resize(maxNx, maxNy);

    double dx = 0.5*isochore->GetDX();
    double dy = 0.5*isochore->GetDY();

    for(int j=0 ; j<static_cast<int>(isochore->GetNJ()) ; j++)
    {
      for(int i=0 ; i<static_cast<int>(isochore->GetNI()) ; i++)
      {
        double x, y;
        isochore->GetXY(i,j,x,y);
        double tTop  = timeSimbox->getTop(x,y);
        double tBase = timeSimbox->getBot(x,y);
        int    nz    = timeSimbox->getnz();
        double dt    = (tBase - tTop)/(2000.0*static_cast<double>(nz));

        int ii, jj;
        timeSimbox->getIndexes(x,y,ii,jj);

        if(ii!=IMISSING && jj!=IMISSING)
        {
          double sum = 0.0;
          for(int k=0 ; k<nz ; k++)
            sum += velocity->GetValue(ii,jj,k);
          (*isochore)(i,j) = sum*dt;
        }
        else
        {
          int i1, i2, i3, i4, j1, j2, j3, j4;
          timeSimbox->getIndexes(x+dx, y+dy, i1, j1);
          timeSimbox->getIndexes(x-dx, y-dy, i2, j2);
          timeSimbox->getIndexes(x+dx, y-dy, i3, j3);
          timeSimbox->getIndexes(x-dx, y+dy, i4, j4);

          int n = 0;
          if(i1!=IMISSING && j1!=IMISSING)
            n++;
          if(i2!=IMISSING && j2!=IMISSING)
            n++;
          if(i3!=IMISSING && j3!=IMISSING)
            n++;
          if(i4!=IMISSING && j4!=IMISSING)
            n++;

          if (n==0)
            isochore->SetMissing(i,j);
          else
          {
            float sum = 0.0f;
            for(int k=0 ; k<nz ; k++) {
              if(i1!=IMISSING && j1!=IMISSING)
                sum += velocity->GetValue(i1,j1,k);
              if(i2!=IMISSING && j2!=IMISSING)
                sum += velocity->GetValue(i2,j2,k);
              if(i3!=IMISSING && j3!=IMISSING)
                sum += velocity->GetValue(i3,j3,k);
              if(i4!=IMISSING && j4!=IMISSING)
                sum += velocity->GetValue(i4,j4,k);
            }
            (*isochore)(i,j) = static_cast<double>(sum)*dt/static_cast<double>(n);
          }
        }
      }
    }

    if(z0Grid_==NULL)
    {
      z0Grid_ = new Surface(*isochore);
      z0Grid_->Multiply(-1.0);
      z0Grid_->AddNonConform(z1Grid_);
    }
    else
    {
      z1Grid_ = new Surface(*isochore);
      z1Grid_->AddNonConform(z0Grid_);
    }
    delete isochore;
  }
}

// Constructor with two correlation surfaces -------------------------
Simbox::Simbox(const Simbox         * simbox,
               const std::string    & interval_name,
               int                    n_layers,
               double                 ext_dz,
               double                 lz_limit,
               const Surface        & top_surface,
               const Surface        & base_surface,
               const Surface        * top_corr_surface,
               const Surface        * base_corr_surface,
               int                    other_output,
               int                    output_domain,
               int                    output_format,
               std::string          & err_text,
               bool                 & failed)
: Volume(*simbox),
top_eroded_surface_(NULL),
base_eroded_surface_(NULL)
{
  std::string output_name = "";
  if (interval_name != "")
    output_name = " for interval \'" + interval_name + "\'";
  LogKit::LogFormatted(LogKit::Low,"\nCreating a temporary simbox (used for inversion grid) with two correlation surfaces" + output_name + ".\n");

  interval_name_  = interval_name;
  status_         = BOXOK;
  cosrot_         = cos(simbox->GetAngle());
  sinrot_         = sin(simbox->GetAngle());
  dx_             = simbox->getdx();
  dy_             = simbox->getdy();
  dz_             = -1;
  nx_             = simbox->getnx();
  ny_             = simbox->getny();
  nz_             = n_layers;
  nx_pad_         = nx_;
  ny_pad_         = ny_;
  nz_pad_         = nz_;
  x_pad_fac_      = (nx_pad_ / nx_);
  y_pad_fac_      = (ny_pad_ / ny_);
  z_pad_fac_      = (nz_pad_ / nz_);
  lz_eroded_      = 0;
  inLine0_        = simbox->getIL0();
  crossLine0_     = simbox->getXL0();
  ilStepX_        = simbox->getILStepX();
  ilStepY_        = simbox->getILStepY();
  xlStepX_        = simbox->getXLStepX();
  xlStepY_        = simbox->getXLStepY();
  constThick_     = simbox->getIsConstantThick();

  setDepth(top_surface, base_surface, n_layers);
  SetErodedSurfaces(top_surface, base_surface);
  //this->calculateDz(lz_limit, err_text);
  double old_dz = GetLZ() / nz_;

  //
  // Check that the two surfaces do not intersect
  //
  double x,y;
  for (size_t i = 0; i < base_surface.GetNI(); i++){
    for (size_t j = 0; j < base_surface.GetNJ(); j++){

      base_surface.GetXY(i,j,x,y);
      double z_top_corr   = top_corr_surface->GetZ(x,y);
      double z_base_corr  = base_corr_surface->GetZ(x,y);

      if (!top_corr_surface->IsMissing(z_top_corr) && !base_corr_surface->IsMissing(z_base_corr) && z_top_corr > z_base_corr) {
        std::string interval_text = "";
        if (interval_name != "")
          interval_text = " for interval " + interval_name;
        err_text += "Error: The top correlation surface crosses the base correlation surface" + interval_text + ".\n";
        failed = true;
      }
    }
  }

  //
  // Check that the correlation surfaces cover the inversion surfaces
  //
  if (simbox->CheckSurface(*top_corr_surface) == false) {
    err_text += "Error: Top correlation surface "+ top_corr_surface->GetName() +"does not cover volume.\n";
    failed = true;
  }
  if (simbox->CheckSurface(*base_corr_surface) == false){
    err_text += "Error: Base correlation surface "+ base_corr_surface->GetName() +"does not cover volume.\n";
    failed = true;
  }

  //
  // Should the corr surfaces have the same resolution?
  //

  if(!failed){
    Surface * mean_corr_surface;
    //NRLib::Vector corr_plane_parameters_top = FindPlane(top_corr_surface);
    //NRLib::Vector corr_plane_parameters_base = FindPlane(base_corr_surface);

    Surface * mean_surface;
    // Use the finest resolution
    double resolution_top   = top_corr_surface->GetDX()*top_corr_surface->GetDY();
    double resolution_base  = base_corr_surface->GetDX()*base_corr_surface->GetDY();
//.........这里部分代码省略.........