/* Pixel.h */

/* 
   Copyright (C) 2006 Fabrice Ducos <fabrice.ducos@icare.univ-lille1.fr>
   This file is part of the Pixel Library.

   The Pixel 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 2.1 of the License, or (at your option) any later version.

   The Pixel 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 the Pixel Library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307 USA.

*/

#ifndef PIXEL_H
#define PIXEL_H

#include <cassert>
#include <set>
#include <algorithm>
#include <utility> // pair
#include <functional> // binary_function
#include <cmath>
#include <vector>
#include <iostream>
#include <sstream>

#define BENCH
#include "debug.h"

template<typename T>
inline T sqr(T x) { return x*x; }

template<typename T, typename V>
class Pixel_base;

template<typename Pixel_type>
class Nearer_from : public std::binary_function<Pixel_type, Pixel_type, bool> {
 public:
  Nearer_from(const Pixel_type &this_pixel) : this_pixel(this_pixel) { }

  bool operator()(const Pixel_type &pixel1, const Pixel_type &pixel2) const {
    return this_pixel.distance(pixel1) < this_pixel.distance(pixel2);
  }

 private:
  const Pixel_type this_pixel;

};

template<typename T, typename V>
class Pixel_base {
  
public:
  typedef Pixel_base<T, V> Pixel_type;
  typedef T coord_type;
  typedef T distance_type;
  typedef V value_type;
  
  static const coord_type DEFAULT_RESOLUTION; // km
  static const coord_type R_EARTH; // Earth's authalic radius ("equal aera"), from http://en.wikipedia.org/wiki/Earth_radius
  static const distance_type DEG2RAD;
  static const distance_type RAD2DEG;
  
  typedef enum {
    RADIANS,
    DEGREES
  } unit_type;
  
  static coord_type get_resolution() { return resolution; }
  static void set_resolution(distance_type new_resolution) { 
    assert(new_resolution > 0);
    resolution = new_resolution; 
  }

  Pixel_base() { /* do nothing */ }

  Pixel_base(const coord_type lat, const coord_type lon, const value_type &val, const unit_type unit = RADIANS) : 
    lat_(lat), 
    lon_(lon),
    val_(val)
  {
    if (unit == DEGREES) {
      lat_ *= DEG2RAD;
      lon_ *= DEG2RAD;
    }
    else {
      assert(unit == RADIANS);
    }

    if ((lat_ >= -M_PI && lat_ <= M_PI) == false || (lon_ >= -2*M_PI && lon_ <= 2*M_PI) == false) {
      std::ostringstream s;
      s << "class Pixel: could not instantiate Pixel(lat:" << lat_*RAD2DEG << ",lon:" << lon_*RAD2DEG << "): out of geographic bounds";
      throw s.str().c_str();
    }

	if (lat_ < -M_PI_2) { lat_ = -( M_PI + lat_); }
    else if (lat_ > M_PI_2) { lat_ = (M_PI - lat_); }
    
	if (lon_ < -M_PI) { lon_ = (2*M_PI + lon_); }
    else if (lon_ > M_PI) { lon_ = -(2*M_PI - lon_); }
  
    ilat = lround((lat_ + M_PI_2)*R_EARTH/(resolution));
    ilon = lround((lon_ + M_PI)*R_EARTH/(resolution)); 

  }

  bool operator< (const Pixel_type &other) const
  { 

    // this is a strict weak ordering

    if (*this == other) return false;
    if (ilat < other.ilat) return true;
    if (ilat > other.ilat) return false;
    return (ilon < other.ilon);
 
  }
  
  bool operator== (const Pixel_type &other) const 
  {
    return (ilat == other.ilat && ilon == other.ilon);
  }

  coord_type lat() const { return lat_; }
  coord_type lon() const { return lon_; }
  value_type val() const { return val_; }
  
  distance_type distance(const Pixel_type &other) const 
  {  
    // Haversine formulation of distances (better than law of cosines formulation which is ill-conditioned for small distances)
    // http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
    distance_type dlat;
    distance_type dlon;
    distance_type a;
    distance_type c;
    distance_type sin_dlat_2;
    distance_type sin_dlon_2;
    
    dlat = other.lat_ - lat_;
    dlon = other.lon_ - lon_;
    
    sin_dlat_2 = sin(dlat/2.);
    sin_dlon_2 = sin(dlon/2.);
    a = sin_dlat_2*sin_dlat_2 + cos(lat_)*cos(other.lat_)*sin_dlon_2*sin_dlon_2;
    c = 2 * asin(std::min(static_cast<T>(1.), static_cast<T>(sqrt(a)) ));

    /*
    distance_type lat = (lat_ + other.lat_)/2.;
    c = sqrt(sqr(dlat) + sqr(dlon * sin(lat)));
    */
    return R_EARTH*c;
   
  }

  void get_neighbours(std::vector<Pixel_type> &neighbours, std::multiset<Pixel_type> &pixels, distance_type resolution,
		      bool sorted = false) const  {

    using std::multiset;
    using std::pair;
    using std::vector;

    typedef multiset<Pixel_type> multiset_type;
    typedef typename multiset_type::iterator pixel_iterator;
    const coord_type dlat = resolution/R_EARTH;
    const coord_type dlon = resolution/R_EARTH;
    
    pair<pixel_iterator, pixel_iterator> range;

    for (coord_type xlat = lat_ - 2*dlat ; xlat <= lat_ + 2*dlat ; xlat += dlat ) {
      for (coord_type xlon = lon_ - 2*dlon ; xlon <= lon_ + 2*dlon ; xlon += dlon) {
	Pixel_type pixel(xlat, xlon, val_);
	range = pixels.equal_range(pixel);
	for (pixel_iterator ppixel = range.first ; ppixel != range.second ; ppixel++) {
	  if (distance(*ppixel) < resolution) { neighbours.push_back(*ppixel); }
	}
      } // xlon
    } // xlat
  
    if (sorted) {
      sort(neighbours.begin(), neighbours.end(), Nearer_from<Pixel_type>(*this));
    }

  }
  

  friend std::ostream &operator<< (std::ostream &os, const Pixel_type& pixel)
  {
    os << "{" << pixel.lat_*RAD2DEG << "°," << pixel.lon_*RAD2DEG << "°}";
    return os;
  }

private:
  static distance_type resolution; /* km */

  coord_type lat_; /* radians */
  coord_type lon_; /* radians */
  value_type val_;


  /* Earth globe is split into small square elements of resolution km along parallels and meridians.
   * Each element is indexed by (ilat, ilon).
   * All the (lat, lon) pairs that fall into the same (ilat, ilon) are considered equivalent.
   * (see operator==)
   */
  int ilat;
  int ilon;
};

template<typename T, typename V>
const typename Pixel_base<T, V>::coord_type Pixel_base<T,V>::DEFAULT_RESOLUTION = 1.; // km
template<typename T, typename V>
const typename Pixel_base<T, V>::coord_type Pixel_base<T,V>::R_EARTH = 6371.005076; // Earth's authalic radius ("equal aera"), from http://en.wikipedia.org/wiki/Earth_radius
template<typename T, typename V>
const typename Pixel_base<T, V>::distance_type Pixel_base<T,V>::DEG2RAD = M_PI/180.;
template<typename T, typename V>
const typename Pixel_base<T, V>::distance_type Pixel_base<T,V>::RAD2DEG = 180.*M_1_PI;
  
template<typename T, typename V> typename
Pixel_base<T, V>::distance_type Pixel_base<T, V>::resolution = DEFAULT_RESOLUTION;


#endif // PIXEL_H