/home/pascal/depot/filedata/src/pixel.h

/* 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 "tools.h"

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 coord_type R_EARTH_EQUATORIAL; // Earth's equatorial radius, from http://en.wikipedia.org/wiki/Earth_radius
  static const coord_type R_EARTH_POLAR; // Earth's polar radius, 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 = DEGREES) :
    lat_(lat),
    lon_(lon),
    val_(val)
  {
    if (unit == DEGREES) {
      lat_ *= DEG2RAD;
      lon_ *= DEG2RAD;
    }
    else {
      assert(unit == RADIANS);
    }

    assert(lat_ >= -M_PI_2 && lat_ <= M_PI_2);
    assert(lon_ >= -M_PI && lon_ <= M_PI);

   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 get_lat() const { return lat_; }
  coord_type get_lon() const { return lon_; }
  value_type get_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)) ));

    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 << " deg," << pixel.lon_*RAD2DEG << " deg}";
    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, class V>
class P_Pixel_base {
public:
  typedef P_Pixel_base<T, V> Pixel_type;
  typedef T coord_type;
  typedef T distance_type;
  typedef V value_type;
private:
  coord_type* lat; /* degrees */
  coord_type* lon; /* degrees */
  value_type val;
public:
  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 coord_type R_EARTH_EQUATORIAL; // Earth's equatorial radius, from http://en.wikipedia.org/wiki/Earth_radius
  static const coord_type R_EARTH_POLAR; // Earth's polar radius, from http://en.wikipedia.org/wiki/Earth_radius
  static const distance_type DEG2RAD;
  static const distance_type RAD2DEG;

  typedef enum {
    RADIANS,
    DEGREES
  } unit_type;

  P_Pixel_base() : lat((coord_type*)(0)),lon((coord_type*)(0)),val(value_type(0)) { /* do nothing */ }

  P_Pixel_base(const coord_type* lat, const coord_type* lon, const value_type &val=value_type(0)) :
    lat(const_cast<coord_type*>(lat)),lon(const_cast<coord_type*>(lon)),val(val){}

  P_Pixel_base(const P_Pixel_base &other):lat(const_cast<coord_type*>(other.lat)),lon(const_cast<coord_type*>(other.lon)),val(other.val){}

  P_Pixel_base & operator= (const P_Pixel_base &other) {
    if (&other == this) return *this; // protect against self assignation
    lat=const_cast<coord_type*>(other.lat);
    lon=const_cast<coord_type*>(other.lon);
    val=other.val;
    return *this;
  }

  bool operator< (const Pixel_type &other) const
  {
    // this is a strict weak ordering
    if (*this == other) return false;
    if ((*lat) < other.get_lat()) return true;
    if ((*lat) > other.get_lat()) return false;
    return ((*lon) < other.get_lon());

  }
  bool operator== (const Pixel_type &other) const
  {
    return (get_lat() == other.get_lat() && get_lon() == other.get_lon());
  }

  coord_type get_lat() const { return coord_type(*lat); }
  coord_type get_lon() const { return coord_type(*lon); }
  value_type get_val() const { return val; }

  distance_type distance(const Pixel_type &other) const
  {
    // use a plane approximation
    return static_cast<distance_type>(sqrt (pow(other.get_lat()-get_lat(),2) + pow(other.get_lon()-get_lon(),2)) );
  }

//   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));
//     }
//   }
//   template < class PPixelContainer >
  void get_neighbours( vector<Pixel_type> &neighbours, const vector<Pixel_type> &pixels,
                       distance_type resolution, bool sorted = false) {

    typedef typename vector<Pixel_type>::const_iterator pixel_iterator;

    /*** use a generic search on sorted datas ***/
    vector<int> default_value(1);
    default_value[0] = -1;

    // bounds of the colocation frame
    float lat_min=*lat-resolution; // lower latitude acceptable
    float lon_min=*lon-resolution; // lower longitude acceptable
    Pixel_type pix_min(&lat_min,&lon_min,default_value);
    float lat_max=*lat+resolution; // biggest latitude acceptable
    float lon_max=*lon+resolution; // biggest longitude acceptable
    Pixel_type pix_max(&lat_max,&lon_max,default_value);

    // set the range of colocated points candidates
    pixel_iterator i_pix_min = lower_bound(pixels.begin(), pixels.end(), pix_min);
    pixel_iterator i_pix_max = upper_bound(pixels.begin(), pixels.end(), pix_max);

//     cout<<"Distance "<<distance(i_pix_min,i_pix_max)<<endl;
    // cross the candidates to match the ones in the resolution frame
    for ( pixel_iterator i_pix = i_pix_min ; i_pix != i_pix_max ; ++i_pix ) {
        if ( i_pix != pixels.end() &&
             abs(i_pix->get_lon()-*lon) < resolution &&
             abs(i_pix->get_lat()-*lat) < resolution )
            // the pixel is in the colocalisation frame
            neighbours.push_back( (*i_pix) );
    }

    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.get_lat() << " deg," << pixel.get_lon() << " deg}";
    return os;
  }
};

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 P_Pixel_base<T, V>::coord_type P_Pixel_base<T,V>::R_EARTH=6371.005076;
template<typename T, typename V>
const typename Pixel_base<T, V>::coord_type Pixel_base<T,V>::R_EARTH_EQUATORIAL = 6378.137; // km
template<typename T, typename V>
const typename P_Pixel_base<T, V>::coord_type P_Pixel_base<T,V>::R_EARTH_EQUATORIAL = 6378.137; // km
template<typename T, typename V>
const typename Pixel_base<T, V>::coord_type Pixel_base<T,V>::R_EARTH_POLAR = 6356.7523; // km
template<typename T, typename V>
const typename P_Pixel_base<T, V>::coord_type P_Pixel_base<T,V>::R_EARTH_POLAR = 6356.7523; // km
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