purify::operators Namespace Reference

Enumerations
enum class	fftw_plan { estimate , measure }
	enum for fftw plans More...

Functions
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	init_on_the_fly_gridding_matrix_2d (const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_complex > &weights, const t_uint &imsizey_, const t_uint &imsizex_, const t_real &oversample_ratio, const std::function< t_real(t_real)> &kernelu, const t_uint Ju, const t_int total_samples)
	on the fly application of the degridding operator using presampling More...
template<class T , class... ARGS>
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	init_gridding_matrix_2d (ARGS &&...args)
	constructs lambdas that apply degridding matrix with adjoint More...
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	init_zero_padding_2d (const Image< typename T::Scalar > &S, const t_real &oversample_ratio)
	Construsts zero padding operator. More...
template<class T >
sopt::OperatorFunction< T >	init_normalise (const t_real &op_norm)
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	init_weights_ (const Vector< t_complex > &weights)
	Construsts zero padding operator. More...
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	init_FFT_2d (const t_uint &imsizey_, const t_uint &imsizex_, const t_real &oversample_factor_, const fftw_plan fftw_plan_flag_=fftw_plan::measure)
	Construsts FFT operator. More...
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	base_padding_and_FFT_2d (const std::function< t_real(t_real)> &ftkernelu, const std::function< t_real(t_real)> &ftkernelv, const t_uint &imsizey, const t_uint &imsizex, const t_real &oversample_ratio=2, const fftw_plan &ft_plan=fftw_plan::measure, const t_real &w_mean=0, const t_real &cellx=1, const t_real &celly=1)
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	base_degrid_operator_2d (const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_real > &w, const Vector< t_complex > &weights, const t_uint &imsizey, const t_uint &imsizex, const t_real &oversample_ratio=2, const kernels::kernel kernel=kernels::kernel::kb, const t_uint Ju=4, const t_uint Jv=4, const fftw_plan &ft_plan=fftw_plan::measure, const bool w_stacking=false, const t_real &cellx=1, const t_real &celly=1, const bool on_the_fly=true)
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	base_padding_and_FFT_2d (const std::function< t_real(t_real)> &ftkerneluv, const t_uint imsizey, const t_uint imsizex, const t_real oversample_ratio, const fftw_plan ft_plan, const t_real w_mean, const t_real cellx, const t_real celly)
template<class T >
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > >	base_degrid_operator_2d (const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_real > &w, const Vector< t_complex > &weights, const t_uint &imsizey, const t_uint &imsizex, const t_real oversample_ratio, const kernels::kernel kernel, const t_uint Ju, const t_uint Jw, const fftw_plan ft_plan, const bool w_stacking, const t_real cellx, const t_real celly, const t_real absolute_error, const t_real relative_error, const dde_type dde)

Enumeration Type Documentation

fftw_plan

enum purify::operators::fftw_plan

strong

enum for fftw plans

Enumerator
estimate
measure

Definition at line 389 of file operators.h.

{ estimate, measure };

Function Documentation

base_degrid_operator_2d() [1/2]

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::base_degrid_operator_2d	(	const Vector< t_real > &	u,
		const Vector< t_real > &	v,
		const Vector< t_real > &	w,
		const Vector< t_complex > &	weights,
		const t_uint &	imsizey,
		const t_uint &	imsizex,
		const t_real &	oversample_ratio = 2,
		const kernels::kernel	kernel = kernels::kernel::kb,
		const t_uint	Ju = 4,
		const t_uint	Jv = 4,
		const fftw_plan &	ft_plan = fftw_plan::measure,
		const bool	w_stacking = false,
		const t_real &	cellx = 1,
		const t_real &	celly = 1,
		const bool	on_the_fly = true
	)

Definition at line 490 of file operators.h.

                                  {
  std::function<t_real(t_real)> kernelu, kernelv, ftkernelu, ftkernelv;
  std::tie(kernelu, kernelv, ftkernelu, ftkernelv) =
      purify::create_kernels(kernel, Ju, Jv, imsizey, imsizex, oversample_ratio);
  sopt::OperatorFunction<T> directFZ, indirectFZ;
  t_real const w_mean = w_stacking ? w.array().mean() : 0.;
  std::tie(directFZ, indirectFZ) = base_padding_and_FFT_2d<T>(
      ftkernelu, ftkernelv, imsizey, imsizex, oversample_ratio, ft_plan, w_mean, cellx, celly);
  sopt::OperatorFunction<T> directG, indirectG;
  PURIFY_MEDIUM_LOG("FoV (width, height): {} deg x {} deg", imsizex * cellx / (60. * 60.),
                    imsizey * celly / (60. * 60.));
  PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG");
  PURIFY_MEDIUM_LOG("Number of visibilities: {}", u.size());
  PURIFY_MEDIUM_LOG("Mean, w: {}, +/- {}", w_mean, (w.maxCoeff() - w.minCoeff()) * 0.5);
  std::tie(directG, indirectG) =
      (on_the_fly)
          ? purify::operators::init_on_the_fly_gridding_matrix_2d<T>(
                u, v, weights, imsizey, imsizex, oversample_ratio, kernelu, Ju, 4e5)
          : purify::operators::init_gridding_matrix_2d<T>(
                u, v, weights, imsizey, imsizex, oversample_ratio, kernelv, kernelu, Ju, Jv);
  auto direct = sopt::chained_operators<T>(directG, directFZ);
  auto indirect = sopt::chained_operators<T>(indirectFZ, indirectG);
  PURIFY_LOW_LOG("Finished consturction of Φ.");
  return std::make_tuple(direct, indirect);
}

References purify::create_kernels(), PURIFY_LOW_LOG, PURIFY_MEDIUM_LOG, operators_test::u, operators_test::v, and purify::utilities::w_stacking().

base_degrid_operator_2d() [2/2]

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::base_degrid_operator_2d	(	const Vector< t_real > &	u,
		const Vector< t_real > &	v,
		const Vector< t_real > &	w,
		const Vector< t_complex > &	weights,
		const t_uint &	imsizey,
		const t_uint &	imsizex,
		const t_real	oversample_ratio,
		const kernels::kernel	kernel,
		const t_uint	Ju,
		const t_uint	Jw,
		const fftw_plan	ft_plan,
		const bool	w_stacking,
		const t_real	cellx,
		const t_real	celly,
		const t_real	absolute_error,
		const t_real	relative_error,
		const dde_type	dde
	)

Definition at line 60 of file wproj_operators.h.

                                                                                  {
  sopt::OperatorFunction<T> directFZ, indirectFZ;
  sopt::OperatorFunction<T> directG, indirectG;
  t_real const w_mean = w_stacking ? w.array().mean() : 0;
  switch (dde) {
  case (dde_type::wkernel_radial): {
    auto const kerneluvs = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio);
    std::tie(directFZ, indirectFZ) = base_padding_and_FFT_2d<T>(
        std::get<0>(kerneluvs), imsizey, imsizex, oversample_ratio, ft_plan, w_mean, cellx, celly);
    PURIFY_MEDIUM_LOG("FoV (width, height): {} deg x {} deg", imsizex * cellx / (60. * 60.),
                      imsizey * celly / (60. * 60.));
    PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG");
    PURIFY_MEDIUM_LOG("Number of visibilities: {}", u.size());
    PURIFY_MEDIUM_LOG("Mean, w: {}, +/- {}", w.array().mean(), (w.maxCoeff() - w.minCoeff()) * 0.5);
    std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_2d<T>(
        u, v, w.array() - w_mean, weights, imsizey, imsizex, oversample_ratio,
        std::get<0>(kerneluvs), std::get<1>(kerneluvs), Ju, Jw, cellx, celly, absolute_error,
        relative_error, dde);
    break;
  }
  case (dde_type::wkernel_2d): {
    std::function<t_real(t_real)> kernelu, kernelv, ftkernelu, ftkernelv;
    std::tie(kernelu, kernelv, ftkernelu, ftkernelv) =
        purify::create_kernels(kernel, Ju, Ju, imsizey, imsizex, oversample_ratio);
    std::tie(directFZ, indirectFZ) = base_padding_and_FFT_2d<T>(
        ftkernelu, ftkernelv, imsizey, imsizex, oversample_ratio, ft_plan, w_mean, cellx, celly);
    PURIFY_MEDIUM_LOG("FoV (width, height): {} deg x {} deg", imsizex * cellx / (60. * 60.),
                      imsizey * celly / (60. * 60.));
    PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG");
    PURIFY_MEDIUM_LOG("Number of visibilities: {}", u.size());
    PURIFY_MEDIUM_LOG("Mean, w: {}, +/- {}", w.array().mean(), (w.maxCoeff() - w.minCoeff()) * 0.5);
    auto const kerneluvs = purify::create_radial_ftkernel(kernel, Ju, oversample_ratio);
    std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_2d<T>(
        u, v, w.array() - w_mean, weights, imsizey, imsizex, oversample_ratio,
        std::get<0>(kerneluvs), std::get<1>(kerneluvs), Ju, Jw, cellx, celly, absolute_error,
        relative_error, dde);
    break;
  }
  default:
    throw std::runtime_error("DDE is not recognised.");
  }
  auto direct = sopt::chained_operators<T>(directG, directFZ);
  auto indirect = sopt::chained_operators<T>(indirectFZ, indirectG);
  PURIFY_LOW_LOG("Finished consturction of Φ.");
  return std::make_tuple(direct, indirect);
}

References purify::create_kernels(), purify::create_radial_ftkernel(), PURIFY_LOW_LOG, PURIFY_MEDIUM_LOG, operators_test::u, operators_test::v, purify::utilities::w_stacking(), purify::wkernel_2d, and purify::wkernel_radial.

base_padding_and_FFT_2d() [1/2]

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::base_padding_and_FFT_2d	(	const std::function< t_real(t_real)> &	ftkernelu,
		const std::function< t_real(t_real)> &	ftkernelv,
		const t_uint &	imsizey,
		const t_uint &	imsizex,
		const t_real &	oversample_ratio = 2,
		const fftw_plan &	ft_plan = fftw_plan::measure,
		const t_real &	w_mean = 0,
		const t_real &	cellx = 1,
		const t_real &	celly = 1
	)

Definition at line 453 of file operators.h.

                                                      {
  sopt::OperatorFunction<T> directZ, indirectZ;
  sopt::OperatorFunction<T> directFFT, indirectFFT;
 
  const Image<t_complex> S =
      purify::details::init_correction2d(oversample_ratio, imsizey, imsizex, ftkernelu, ftkernelv,
                                         w_mean, cellx, celly) *
      std::sqrt(imsizex * imsizey) * oversample_ratio;
  PURIFY_LOW_LOG("Building Measurement Operator: WGFZDB");
  PURIFY_LOW_LOG(
      "Constructing Zero Padding "
      "and Correction Operator: "
      "ZDB");
  PURIFY_MEDIUM_LOG("Image size (width, height): {} x {}", imsizex, imsizey);
  PURIFY_MEDIUM_LOG("Oversampling Factor: {}", oversample_ratio);
  std::tie(directZ, indirectZ) = purify::operators::init_zero_padding_2d<T>(S, oversample_ratio);
  PURIFY_LOW_LOG("Constructing FFT operator: F");
  switch (ft_plan) {
  case fftw_plan::measure:
    PURIFY_MEDIUM_LOG("Measuring Plans...");
    break;
  case fftw_plan::estimate:
    PURIFY_MEDIUM_LOG("Estimating Plans...");
    break;
  }
  std::tie(directFFT, indirectFFT) =
      purify::operators::init_FFT_2d<T>(imsizey, imsizex, oversample_ratio, ft_plan);
  auto direct = sopt::chained_operators<T>(directFFT, directZ);
  auto indirect = sopt::chained_operators<T>(indirectZ, indirectFFT);
  return std::make_tuple(direct, indirect);
}

References estimate, purify::details::init_correction2d(), measure, PURIFY_LOW_LOG, and PURIFY_MEDIUM_LOG.

base_padding_and_FFT_2d() [2/2]

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::base_padding_and_FFT_2d	(	const std::function< t_real(t_real)> &	ftkerneluv,
		const t_uint	imsizey,
		const t_uint	imsizex,
		const t_real	oversample_ratio,
		const fftw_plan	ft_plan,
		const t_real	w_mean,
		const t_real	cellx,
		const t_real	celly
	)

Definition at line 25 of file wproj_operators.h.

                        {
  sopt::OperatorFunction<T> directZ, indirectZ;
  sopt::OperatorFunction<T> directFFT, indirectFFT;
  const Image<t_complex> S =
      purify::details::init_correction_radial_2d(oversample_ratio, imsizey, imsizex, ftkerneluv,
                                                 w_mean, cellx, celly) *
      std::sqrt(imsizex * imsizey) * oversample_ratio;
  PURIFY_LOW_LOG("Building Measurement Operator: WGFZDB");
  PURIFY_LOW_LOG(
      "Constructing Zero Padding "
      "and Correction Operator: "
      "ZDB");
  PURIFY_MEDIUM_LOG("Image size (width, height): {} x {}", imsizex, imsizey);
  PURIFY_MEDIUM_LOG("Oversampling Factor: {}", oversample_ratio);
  std::tie(directZ, indirectZ) = purify::operators::init_zero_padding_2d<T>(S, oversample_ratio);
  PURIFY_LOW_LOG("Constructing FFT operator: F");
  switch (ft_plan) {
  case fftw_plan::measure:
    PURIFY_MEDIUM_LOG("Measuring Plans...");
    break;
  case fftw_plan::estimate:
    PURIFY_MEDIUM_LOG("Estimating Plans...");
    break;
  }
  std::tie(directFFT, indirectFFT) =
      purify::operators::init_FFT_2d<T>(imsizey, imsizex, oversample_ratio, ft_plan);
  auto direct = sopt::chained_operators<T>(directFFT, directZ);
  auto indirect = sopt::chained_operators<T>(indirectZ, indirectFFT);
  return std::make_tuple(direct, indirect);
}

References estimate, purify::details::init_correction_radial_2d(), measure, PURIFY_LOW_LOG, and PURIFY_MEDIUM_LOG.

init_FFT_2d()

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::init_FFT_2d	(	const t_uint &	imsizey_,
		const t_uint &	imsizex_,
		const t_real &	oversample_factor_,
		const fftw_plan	fftw_plan_flag_ = fftw_plan::measure
	)

Construsts FFT operator.

Definition at line 392 of file operators.h.

                                                        {
  t_int const ftsizeu_ = std::floor(imsizex_ * oversample_factor_);
  t_int const ftsizev_ = std::floor(imsizey_ * oversample_factor_);
  t_int plan_flag = (FFTW_MEASURE | FFTW_PRESERVE_INPUT);
  switch (fftw_plan_flag_) {
  case (fftw_plan::measure):
    plan_flag = (FFTW_MEASURE | FFTW_PRESERVE_INPUT);
    break;
  case (fftw_plan::estimate):
    plan_flag = (FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
    break;
  }
 
#ifdef PURIFY_OPENMP_FFTW
  PURIFY_LOW_LOG("Using OpenMP threading with FFTW.");
  fftw_init_threads();
#endif
  Vector<typename T::Scalar> src = Vector<t_complex>::Zero(ftsizev_ * ftsizeu_);
  Vector<typename T::Scalar> dst = Vector<t_complex>::Zero(ftsizev_ * ftsizeu_);
  // creating plans
  const auto del = [](fftw_plan_s *plan) { fftw_destroy_plan(plan); };
  // fftw plan with threads needs to be used before each fftw_plan is created
#ifdef PURIFY_OPENMP_FFTW
  fftw_plan_with_nthreads(omp_get_max_threads());
#endif
  const std::shared_ptr<fftw_plan_s> m_plan_forward(
      fftw_plan_dft_2d(ftsizev_, ftsizeu_, reinterpret_cast<fftw_complex *>(src.data()),
                       reinterpret_cast<fftw_complex *>(dst.data()), FFTW_FORWARD, plan_flag),
      del);
  // fftw plan with threads needs to be used before each fftw_plan is created
#ifdef PURIFY_OPENMP_FFTW
  fftw_plan_with_nthreads(omp_get_max_threads());
#endif
  const std::shared_ptr<fftw_plan_s> m_plan_inverse(
      fftw_plan_dft_2d(ftsizev_, ftsizeu_, reinterpret_cast<fftw_complex *>(src.data()),
                       reinterpret_cast<fftw_complex *>(dst.data()), FFTW_BACKWARD, plan_flag),
      del);
  auto const direct = [m_plan_forward, ftsizeu_, ftsizev_](T &output, const T &input) {
    assert(input.size() == ftsizev_ * ftsizeu_);
    output = Matrix<typename T::Scalar>::Zero(input.rows(), input.cols());
    fftw_execute_dft(
        m_plan_forward.get(),
        const_cast<fftw_complex *>(reinterpret_cast<const fftw_complex *>(input.data())),
        reinterpret_cast<fftw_complex *>(output.data()));
    output /= std::sqrt(output.size());
  };
  auto const indirect = [m_plan_inverse, ftsizeu_, ftsizev_](T &output, const T &input) {
    assert(input.size() == ftsizev_ * ftsizeu_);
    output = Matrix<typename T::Scalar>::Zero(input.rows(), input.cols());
    fftw_execute_dft(
        m_plan_inverse.get(),
        const_cast<fftw_complex *>(reinterpret_cast<const fftw_complex *>(input.data())),
        reinterpret_cast<fftw_complex *>(output.data()));
    output /= std::sqrt(output.size());
  };
  return std::make_tuple(direct, indirect);
}

References estimate, measure, and PURIFY_LOW_LOG.

init_gridding_matrix_2d()

template<class T , class... ARGS>

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::init_gridding_matrix_2d

(

ARGS &&...

args

)

constructs lambdas that apply degridding matrix with adjoint

Definition at line 312 of file operators.h.

                    {
  const std::shared_ptr<const Sparse<t_complex>> interpolation_matrix =
      std::make_shared<const Sparse<t_complex>>(
          details::init_gridding_matrix_2d(std::forward<ARGS>(args)...));
  const std::shared_ptr<const Sparse<t_complex>> adjoint =
      std::make_shared<const Sparse<t_complex>>(interpolation_matrix->adjoint());
 
  return std::make_tuple(
      [=](T &output, const T &input) {
        output = utilities::sparse_multiply_matrix(*interpolation_matrix, input);
      },
      [=](T &output, const T &input) {
        output = utilities::sparse_multiply_matrix(*adjoint, input);
      });
}

References purify::details::init_gridding_matrix_2d(), and purify::utilities::sparse_multiply_matrix().

init_normalise()

template<class T >

sopt::OperatorFunction<T> purify::operators::init_normalise

(

const t_real &

op_norm

)

Definition at line 368 of file operators.h.

                                                              {
  if (not(op_norm > 0)) throw std::runtime_error("Operator norm is not greater than zero.");
  return [=](T &output, const T &x) { output = x / op_norm; };
}

init_on_the_fly_gridding_matrix_2d()

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::init_on_the_fly_gridding_matrix_2d	(	const Vector< t_real > &	u,
		const Vector< t_real > &	v,
		const Vector< t_complex > &	weights,
		const t_uint &	imsizey_,
		const t_uint &	imsizex_,
		const t_real &	oversample_ratio,
		const std::function< t_real(t_real)> &	kernelu,
		const t_uint	Ju,
		const t_int	total_samples
	)

on the fly application of the degridding operator using presampling

Definition at line 21 of file fly_operators.h.

                                                                                            {
  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
  const t_uint rows = u.size();
  const t_uint cols = ftsizeu_ * ftsizev_;
  if (u.size() != v.size())
    throw std::runtime_error(
        "Size of u and v vectors are not the same for creating gridding matrix.");
 
  const std::shared_ptr<Vector<t_real>> u_ptr = std::make_shared<Vector<t_real>>(u);
  const std::shared_ptr<Vector<t_real>> v_ptr = std::make_shared<Vector<t_real>>(v);
  const std::shared_ptr<Vector<t_complex>> weights_ptr =
      std::make_shared<Vector<t_complex>>(weights);
  const t_complex I(0, 1);
  const t_int ju_max = std::min(Ju, ftsizeu_);
  const t_int jv_max = std::min(Ju, ftsizev_);
  const auto samples = kernels::kernel_samples(
      total_samples, [&](const t_real x) { return kernelu(x * ju_max * 0.5); });
  std::set<t_int> nonZeros_set;
  for (t_int m = 0; m < rows; ++m) {
    t_complex result = 0;
    const t_real u_val = (*u_ptr)(m);
    const t_real v_val = (*v_ptr)(m);
    const t_real k_u = std::floor(u_val - ju_max * 0.5);
    const t_real k_v = std::floor(v_val - jv_max * 0.5);
    for (t_int jv = 1; jv < jv_max + 1; ++jv) {
      const t_uint p = utilities::mod(k_v + jv, ftsizev_);
      const t_real c_0 =
          static_cast<t_int>(std::floor(2 * std::abs(v_val - (k_v + jv)) * total_samples / jv_max));
      for (t_int ju = 1; ju < ju_max + 1; ++ju) {
        const t_uint q = utilities::mod(k_u + ju, ftsizeu_);
        const t_int i_0 = static_cast<t_int>(
            std::floor(2 * std::abs(u_val - (k_u + ju)) * total_samples / ju_max));
        const t_uint index = utilities::sub2ind(p, q, ftsizev_, ftsizeu_);
        nonZeros_set.insert(index);
      }
    }
  }
 
  std::vector<t_int> nonZeros_vec(nonZeros_set.begin(), nonZeros_set.end());
  std::sort(nonZeros_vec.data(), nonZeros_vec.data() + nonZeros_vec.size());
  SparseVector<t_int> mapping(ftsizev_ * ftsizeu_);
  mapping.reserve(nonZeros_vec.size());
  for (t_int index = 0; index < nonZeros_vec.size(); index++)
    mapping.coeffRef(nonZeros_vec[index]) = index;
  PURIFY_LOW_LOG("Non Zero grid locations: {} ", mapping.nonZeros());
 
  const auto degrid = [rows, ju_max, jv_max, I, u_ptr, v_ptr, weights_ptr, samples, total_samples,
                       ftsizeu_, ftsizev_](T &output, const T &input) {
    output = T::Zero(u_ptr->size());
    assert(input.size() == ftsizeu_ * ftsizev_);
#ifdef PURIFY_OPENMP
#pragma omp parallel for
#endif
    for (t_int m = 0; m < rows; ++m) {
      t_complex result = 0;
      const t_real u_val = (*u_ptr)(m);
      const t_real v_val = (*v_ptr)(m);
      const t_real k_u = std::floor(u_val - ju_max * 0.5);
      const t_real k_v = std::floor(v_val - jv_max * 0.5);
      for (t_int jv = 1; jv < jv_max + 1; ++jv) {
        const t_uint p = utilities::mod(k_v + jv, ftsizev_);
        const t_real c_0 = static_cast<t_int>(
            std::floor(2 * std::abs(v_val - (k_v + jv)) * total_samples / jv_max));
        assert(c_0 >= 0);
        assert(c_0 < total_samples);
        const t_real kernelv_val = samples[c_0] * (1. - (2 * (p % 2)));
        for (t_int ju = 1; ju < ju_max + 1; ++ju) {
          const t_uint q = utilities::mod(k_u + ju, ftsizeu_);
          const t_int i_0 = static_cast<t_int>(
              std::floor(2 * std::abs(u_val - (k_u + ju)) * total_samples / ju_max));
          assert(i_0 >= 0);
          assert(i_0 < total_samples);
          const t_real kernelu_val = samples[i_0] * (1. - (2 * (q % 2)));
          const t_uint index = utilities::sub2ind(p, q, ftsizev_, ftsizeu_);
          const t_real sign = kernelu_val * kernelv_val;
          result += input(index) * sign;
        }
      }
      output(m) = result;
    }
    output.array() *= (*weights_ptr).array();
  };
 
  const auto grid = [rows, ju_max, jv_max, I, u_ptr, v_ptr, weights_ptr, samples, total_samples,
                     ftsizeu_, ftsizev_, mapping, nonZeros_vec](T &output, const T &input) {
    const t_int N = ftsizeu_ * ftsizev_;
    output = T::Zero(N);
#ifdef PURIFY_OPENMP
    t_int const max_threads = omp_get_max_threads();
#else
    t_int const max_threads = 1;
#endif
    T output_compressed = T::Zero(nonZeros_vec.size() * max_threads);
    assert(output.size() == N);
#ifdef PURIFY_OPENMP
#pragma omp parallel for
#endif
    for (t_int m = 0; m < rows; ++m) {
      t_complex result = 0;
#ifdef PURIFY_OPENMP
      const t_int shift = omp_get_thread_num() * nonZeros_vec.size();
#else
      const t_int shift = 0;
#endif
      const t_real u_val = (*u_ptr)(m);
      const t_real v_val = (*v_ptr)(m);
      const t_real k_u = std::floor(u_val - ju_max * 0.5);
      const t_real k_v = std::floor(v_val - jv_max * 0.5);
      const t_complex vis = input(m) * std::conj((*weights_ptr)(m));
      for (t_int jv = 1; jv < jv_max + 1; ++jv) {
        const t_uint p = utilities::mod(k_v + jv, ftsizev_);
        const t_real c_0 = static_cast<t_int>(
            std::floor(2 * std::abs(v_val - (k_v + jv)) * total_samples / jv_max));
        assert(c_0 >= 0);
        assert(c_0 < total_samples);
        const t_real kernelv_val = samples[c_0] * (1. - (2 * (p % 2)));
        for (t_int ju = 1; ju < ju_max + 1; ++ju) {
          const t_uint q = utilities::mod(k_u + ju, ftsizeu_);
          const t_int i_0 = static_cast<t_int>(
              std::floor(2 * std::abs(u_val - (k_u + ju)) * total_samples / ju_max));
          assert(i_0 >= 0);
          assert(i_0 < total_samples);
          const t_real kernelu_val = samples[i_0] * (1. - (2 * (q % 2)));
          const t_int index = utilities::sub2ind(p, q, ftsizev_, ftsizeu_);
          const t_complex result = kernelu_val * kernelv_val * vis;
          output_compressed(mapping.coeff(index) + shift) += result;
        }
      }
    }
    for (t_int m = 1; m < max_threads; m++) {
      const t_int loop_shift = m * nonZeros_vec.size();
      output_compressed.segment(0, nonZeros_vec.size()) +=
          output_compressed.segment(loop_shift, nonZeros_vec.size());
    }
    for (t_int index = 0; index < nonZeros_vec.size(); index++)
      output(nonZeros_vec[index]) += output_compressed(index);
  };
  return std::make_tuple(degrid, grid);
}

References purify::I, purify::kernels::kernel_samples(), purify::utilities::mod(), PURIFY_LOW_LOG, purify::utilities::sub2ind(), operators_test::u, and operators_test::v.

init_weights_()

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::init_weights_

(

const Vector< t_complex > &

weights

)

Construsts zero padding operator.

Definition at line 375 of file operators.h.

                                      {
  PURIFY_DEBUG("Calculating weights: W");
  auto direct = [=](T &output, const T &x) {
    assert(weights.size() == x.size());
    output = weights.array() * x.array();
  };
  auto indirect = [=](T &output, const T &x) {
    assert(weights.size() == x.size());
    output = weights.conjugate().array() * x.array();
  };
  return std::make_tuple(direct, indirect);
}

References PURIFY_DEBUG.

init_zero_padding_2d()

template<class T >

std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T> > purify::operators::init_zero_padding_2d	(	const Image< typename T::Scalar > &	S,
		const t_real &	oversample_ratio
	)

Construsts zero padding operator.

Definition at line 331 of file operators.h.

                                                                      {
  const t_uint imsizex_ = S.cols();
  const t_uint imsizey_ = S.rows();
  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
  const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5);
  const t_uint y_start = std::floor(ftsizev_ * 0.5 - imsizey_ * 0.5);
  auto direct = [=](T &output, const T &x) {
    assert(x.size() == imsizex_ * imsizey_);
    output = Vector<t_complex>::Zero(ftsizeu_ * ftsizev_);
#pragma omp parallel for collapse(2)
    for (t_uint j = 0; j < imsizey_; j++) {
      for (t_uint i = 0; i < imsizex_; i++) {
        const t_uint input_index = utilities::sub2ind(j, i, imsizey_, imsizex_);
        const t_uint output_index =
            utilities::sub2ind(y_start + j, x_start + i, ftsizev_, ftsizeu_);
        output(output_index) = S(j, i) * x(input_index);
      }
    }
  };
  auto indirect = [=](T &output, const T &x) {
    assert(x.size() == ftsizeu_ * ftsizev_);
    output = T::Zero(imsizey_ * imsizex_);
#pragma omp parallel for collapse(2)
    for (t_uint j = 0; j < imsizey_; j++) {
      for (t_uint i = 0; i < imsizex_; i++) {
        const t_uint output_index = utilities::sub2ind(j, i, imsizey_, imsizex_);
        const t_uint input_index = utilities::sub2ind(y_start + j, x_start + i, ftsizev_, ftsizeu_);
        output(output_index) = std::conj(S(j, i)) * x(input_index);
      }
    }
  };
  return std::make_tuple(direct, indirect);
}

References purify::utilities::sub2ind().