nonuniform_spectrum_sampling.cc

Go to the documentation of this file.

#include "nonuniform_spectrum_sampling.h"
#include "ostream_pad.h"
#include <algorithm>

using namespace FullPhysics;
using namespace std;
using namespace blitz;

#ifdef HAVE_LUA
#include "register_lua.h"
REGISTER_LUA_DERIVED_CLASS(NonuniformSpectrumSampling, SpectrumSampling)
.def(luabind::constructor<const SpectralDomain&,
                          const SpectralDomain&,
                          const SpectralDomain&,
                          const boost::shared_ptr<SpectrumSampling>&>())
REGISTER_LUA_END()
#endif



//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

NonuniformSpectrumSampling::NonuniformSpectrumSampling(
const SpectralDomain& Grid,
const boost::shared_ptr<SpectrumSampling>& Interpolated_sampling
)
: SpectrumSampling(1), interpolated_sampling(Interpolated_sampling)
{
  spec_domain.push_back(sort_sd(Grid));
}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

NonuniformSpectrumSampling::NonuniformSpectrumSampling(
const SpectralDomain& Grid1,
const SpectralDomain& Grid2,
const SpectralDomain& Grid3,
const boost::shared_ptr<SpectrumSampling>& Interpolated_sampling
)
: SpectrumSampling(3), interpolated_sampling(Interpolated_sampling)
{
  spec_domain.push_back(sort_sd(Grid1));
  spec_domain.push_back(sort_sd(Grid2));
  spec_domain.push_back(sort_sd(Grid3));
}

SpectralDomain NonuniformSpectrumSampling::sort_sd
(const SpectralDomain& In) const
{
  Array<double, 1> res(In.data().copy());
  sort(res.data(), res.data()+res.size());
  return SpectralDomain(res, In.units());
}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

void NonuniformSpectrumSampling::print(std::ostream& Os) const 
{ 
  OstreamPad opad(Os, "    ");
  Os << "NonuniformSpectrumSampling\n";
  Os << "  Uniform spectrum sampling:\n";
  opad << *interpolated_sampling << "\n";
  opad.strict_sync();
  for(int i = 0; i < number_spectrometer(); ++i)
    if(spec_domain[i].data().rows() <= 0)
      Os << "  Band " << i + 1 << ":\n"
         << "     Nonuniform RT grid not provided, using uniform grid for RT\n";
    else
      Os << "  Band " << i + 1 << ":\n"
         << "     grid_start: " << spec_domain[i].data()(0) << "\n"
         << "     grid_end:   " << spec_domain[i].data()(spec_domain[i].data().rows() - 1) << "\n"
         << "     grid_points: " << spec_domain[i].data().rows() << "\n";
}

// See base class for description
SpectralDomain NonuniformSpectrumSampling::spectral_domain
(int spec_index,
 const SpectralDomain& Lowres_grid, 
 const DoubleWithUnit& Ils_half_width) const
{
  range_check(spec_index, 0, number_spectrometer());
  SpectralDomain ispec = interpolated_sampling->
    spectral_domain(spec_index, Lowres_grid, Ils_half_width);
  // If we aren't actually doing nonuniform sampling, just return the
  // spectral domain we want to sample to.
  if(spec_domain[spec_index].data().rows() == 0)
    return ispec;
  // Find all points in ispec that are closest to one of the grid points.
  std::map<double, bool> res;
  typedef std::map<double, bool>::iterator it_type;
  typedef std::map<double, bool>::value_type val_type;
  BOOST_FOREACH(double x, ispec.data())
    res[x] = false;
  BOOST_FOREACH(double x, spec_domain[spec_index].convert_wave(ispec.units())) {
    it_type i = res.lower_bound(x);
    if(i != res.begin()) {
      it_type i2 = i;
      --i2;
      if(fabs(x - (*i2).first) < fabs(x - (*i).first))
        i = i2;
    }
    (*i).second = true;
  }
  // And use those points to make up the spectral domain.
  std::vector<double> resx;
  BOOST_FOREACH(const val_type& i, res)
    if(i.second)
      resx.push_back(i.first);
  Array<double, 1> sd(&resx[0], shape((int) resx.size()), duplicateData);
  return SpectralDomain(sd, ispec.units());
}