/*//////////////////////////////////////////////////////////////////////////////
// CMT Cosmic Muon Tomography project //////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
  All rights reserved

  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >

  ------------------------------------------------------------------
  This library is free software;  you  can  redistribute  it  and/or
  modify it  under the  terms  of  the  GNU  Lesser  General  Public
  License as published  by  the  Free  Software  Foundation;  either
  version 3.0 of the License, or (at your option) any later version.

  This library is  distributed in  the hope that it will  be useful,
  but  WITHOUT ANY WARRANTY;  without  even  the implied warranty of
  MERCHANTABILITY  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of  the GNU Lesser General  Public
  License along with this library.

//////////////////////////////////////////////////////////////////////////////*/

#ifndef VOXRAYTRACER_H
#define VOXRAYTRACER_H

#include <Core/DataAllocator.h>
#include <Core/Vector.h>
#include <math.h>
#include <vector>

#include "Math/StructuredGrid.h"
#include "Math/VoxImage.h"

#ifdef USE_CUDA
#include <cuda_runtime.h>
#endif

namespace uLib {

class VoxRaytracer {

public:
  class RayData {
  public:
    RayData() : m_TotalLength(0), m_Count(0) {}

    struct Element {
      Id_t vox_id;
      Scalarf L;
      ~Element() {}
    };

    inline void AddElement(Id_t id, float L);

    void AppendRay(const RayData &in);

    inline uLib::Vector<Element> &Data() { return this->m_Data; }

    inline const uLib::Vector<Element> &Data() const { return this->m_Data; }

    inline size_t Count() const { return this->m_Count; }

    inline size_t size() const { return this->m_Count; }

    inline const Scalarf &TotalLength() const { return this->m_TotalLength; }

    inline void SetCount(size_t c) {
      this->m_Count = c;
      if (this->m_Data.size() != c) {
        this->m_Data.resize(c);
      }
    }

    inline void SetTotalLength(Scalarf tl) { this->m_TotalLength = tl; }

    void PrintSelf(std::ostream &o);

  private:
    uLib::Vector<Element> m_Data;
    Scalarf m_TotalLength;
    size_t m_Count;
  };

public:
  VoxRaytracer(StructuredGrid &image) : m_Image(&image) {
    m_scale << (m_Image->GetWorldMatrix() * Vector4f(1, 0, 0, 0)).norm(),
        (m_Image->GetWorldMatrix() * Vector4f(0, 1, 0, 0)).norm(),
        (m_Image->GetWorldMatrix() * Vector4f(0, 0, 1, 0)).norm();
  }

  bool GetEntryPoint(const HLine3f &line, HPoint3f &pt);

  bool GetExitPoint(const HLine3f &line, HPoint3f &pt);

  RayData TraceBetweenPoints(const HPoint3f &in, const HPoint3f &out) const;

  RayData TraceLine(const HLine3f &line) const;

  inline StructuredGrid *GetImage() const { return this->m_Image; }

#ifdef USE_CUDA
  template <typename VoxelT>
  void AccumulateLinesCUDA(const HLine3f *lines, size_t num_lines,
                           VoxImage<VoxelT> &image);

  void TraceLineCUDA(const HLine3f *lines, size_t num_lines, RayData *out_rays,
                     int max_elements_per_ray = 128,
                     float *kernel_time_ms = nullptr);

  void TraceBetweenPointsCUDA(const HPoint3f *in_pts, const HPoint3f *out_pts,
                              size_t num_lines, RayData *out_rays,
                              int max_elements_per_ray = 128,
                              float *kernel_time_ms = nullptr);
#endif

private:
  StructuredGrid *m_Image;
  Vector3f m_scale;
};

} // namespace uLib

#ifdef USE_CUDA
#include "Math/VoxRaytracerCUDA.hpp"
#endif

#endif // VOXRAYTRACER_H