#ifndef RSTBX_SIMAGE_IMAGE_SIMPLE_HPP
#define RSTBX_SIMAGE_IMAGE_SIMPLE_HPP

#include <rstbx/import_scitbx_af.hpp>
#include <cctbx/miller.h>
#include <cctbx/uctbx.h>
#include <scitbx/mat3.h>
#include <scitbx/vec2.h>
#include <scitbx/array_family/versa.h>
#include <tbxx/error_utils.hpp>

namespace rstbx { namespace simage {

  struct image_simple
  {
    bool apply_proximity_filter;
    bool apply_detector_clipping;
    bool apply_proximity_factor;
    bool store_miller_index_i_seqs;
    bool store_spots;
    bool store_signals;
    bool set_pixels;
    af::shared<std::size_t> miller_index_i_seqs;
    af::shared<scitbx::vec3<double> > spots;
    af::shared<double> signals;
    af::versa<int, af::flex_grid<> > pixels;

    image_simple(
      bool apply_proximity_filter_,
      bool apply_detector_clipping_,
      bool apply_proximity_factor_,
      bool store_miller_index_i_seqs_,
      bool store_spots_,
      bool store_signals_,
      bool set_pixels_)
    :
      apply_proximity_filter(apply_proximity_filter_),
      apply_detector_clipping(apply_detector_clipping_),
      apply_proximity_factor(apply_proximity_factor_),
      store_miller_index_i_seqs(store_miller_index_i_seqs_),
      store_spots(store_spots_),
      store_signals(store_signals_),
      set_pixels(set_pixels_)
    {}

    image_simple&
    compute(
      cctbx::uctbx::unit_cell const& unit_cell,
      af::const_ref<cctbx::miller::index<> > const& miller_indices,
      af::const_ref<double> const& spot_intensity_factors,
      scitbx::mat3<double> const& crystal_rotation_matrix,
      double ewald_radius,
      double ewald_proximity,
      int signal_max,
      double detector_distance,
      scitbx::vec2<double> detector_size,
      scitbx::vec2<int> detector_pixels,
      unsigned point_spread,
      double gaussian_falloff_scale)
    {
      if (spot_intensity_factors.size() != 0) {
        TBXX_ASSERT(spot_intensity_factors.size() == miller_indices.size());
        TBXX_ASSERT(spot_intensity_factors.all_ge(0));
        TBXX_ASSERT(spot_intensity_factors.all_le(1));
      }
      TBXX_ASSERT(ewald_proximity > 0);
      TBXX_ASSERT(ewald_radius > 0);
      TBXX_ASSERT(detector_size.const_ref().all_gt(0));
      TBXX_ASSERT(detector_pixels.const_ref().all_gt(0));
      TBXX_ASSERT(point_spread > 0);
      TBXX_ASSERT(gaussian_falloff_scale >= 0);
      int dpx = detector_pixels[0];
      int dpy = detector_pixels[1];
      if (set_pixels) {
        pixels.resize(af::flex_grid<>(dpx, dpy), 0);
      }
      int* pixels_beg = pixels.begin();
      double dsx = detector_size[0];
      double dsy = detector_size[1];
      unsigned point_spread_half = point_spread / 2;
      bool point_spread_is_even_value = (point_spread % 2 == 0);
      double circle_radius_sq = point_spread * std::max(dsx/dpx, dsy/dpy) / 2;
      circle_radius_sq *= circle_radius_sq;
      TBXX_ASSERT(circle_radius_sq != 0);
      typedef scitbx::vec3<double> v3d;
      for(std::size_t ih=0;ih<miller_indices.size();ih++) {
        // compare with ewald_proximity() code in explore_completeness.py
        v3d rv = unit_cell.reciprocal_space_vector(miller_indices[ih]);
        v3d rvre = crystal_rotation_matrix * rv;
        rvre[2] += ewald_radius; // direct beam anti-parallel (0,0,1)
        double rvre_len = rvre.length();
        double rvre_proximity = std::abs(1 - rvre_len / ewald_radius);
        if (!apply_proximity_filter || rvre_proximity < ewald_proximity) {
          if (rvre[2] <= 0) {
            TBXX_ASSERT(apply_proximity_filter);
          }
          else {
            // http://en.wikipedia.org/wiki/Line-plane_intersection
            double d = -detector_distance / rvre[2];
            double dx = rvre[0] * d;
            double dy = rvre[1] * d;
            bool off_detector = (
                 std::abs(dx) > dsx/2
              || std::abs(dy) > dsy/2);
            if (off_detector && apply_detector_clipping) {
              continue;
            }
            double pxf = (dx/dsx + 0.5) * dpx;
            double pyf = (dy/dsy + 0.5) * dpy;
            if (store_miller_index_i_seqs) {
              miller_index_i_seqs.push_back(ih);
            }
            if (store_spots) {
              spots.push_back(scitbx::vec3<double>(pxf, pyf, 0));
              if (!store_signals && !set_pixels) continue;
            }
            double proximity_factor = 1;
            if (apply_proximity_factor) {
              proximity_factor -= scitbx::fn::pow2(
                rvre_proximity / ewald_proximity);
              if (proximity_factor < 0) proximity_factor = 0;
            }
            double signal_at_center =
                signal_max
              * (spot_intensity_factors.size() == 0 ? 1 :
                 spot_intensity_factors[ih])
              * proximity_factor;
            if (store_signals) {
              signals.push_back(signal_at_center);
            }
            if (off_detector || !set_pixels || proximity_factor == 0) continue;
            int signal = static_cast<int>(signal_at_center + 0.5);
            double gauss_arg_term = -gaussian_falloff_scale
                                  / circle_radius_sq;
            using scitbx::math::ifloor;
            int pxi = ifloor(pxf);
            int pyi = ifloor(pyf);
            int pxb = pxi - point_spread_half;
            int pyb = pyi - point_spread_half;
            if (point_spread_is_even_value) {
              if (pxf - pxi > 0.5) pxb++;
              if (pyf - pyi > 0.5) pyb++;
            }
            for(int i=0;i<=point_spread;i++) {
              int pi = pxb + i;
              if (pi < 0 || pi >= dpx) continue;
              int pi0 = pi * dpy;
              for(int j=0;j<=point_spread;j++) {
                int pj = pyb + j;
                if (pj < 0 || pj >= dpy) continue;
                if (point_spread > 2) {
                  double pcx = ((pi + 0.5) / dpx - 0.5) * dsx - dx;
                  double pcy = ((pj + 0.5) / dpy - 0.5) * dsy - dy;
                  double pc_sq = pcx*pcx + pcy*pcy;
                  if (pc_sq > circle_radius_sq) continue;
                  if (gaussian_falloff_scale != 0) {
                    double falloff_factor = std::exp(pc_sq * gauss_arg_term);
                    signal = static_cast<int>(
                      signal_at_center * falloff_factor + 0.5);
                  }
                }
                pixels_beg[pi0+pj] = signal;
              }
            }
          }
        }
      }
      return *this;
    }
  };

}} // namespace rstbx::simage

#endif // GUARD