/**
 * @package de.atwillys.cc.swl.im
 * @license BSD (simplified)
 * @author Stefan Wilhelm (stfwi)
 *
 * @file cvimage.hh
 * @ccflags
 * @ldflags -lm -lopencv_core -lopencv_highgui -lopencv_imgproc
 * @platform linux, bsd, windows
 * @standard >= c++98
 *
 * -----------------------------------------------------------------------------
 *
 * Template based wrapper for OpenCV image processing functionality. Brief
 * overview (main image class is `cvimage`), example instance `img`. Most of
 * the methods have multiple overloads that are not shown here. The detailed
 * documentation about each method is in JavaDoc/Doxygen format in the source
 * code:
 *
 * - I/O, methods known from STL containers:
 *
 *  - img.load(filename, [as-is, gray, color]); // Load
 *  - img.save([filename]);                  // Save (file omitted===load name)
 *  - img.clear();                           // Reset image, clear matrix
 *  - if(img.empty()) ...                    // Empty image?
 *
 * - Basic information
 *
 *   - w = img.width();
 *   - h = img.height();
 *   - A = img.area();
 *   - n = img.channels();                    // Number of channels
 *   - if(img.is_rgb()) {}
 *   - if(img.is_gray()) {}
 *   - if(img.is_hsv()) {}
 *
 * - Information requiring calculations:
 *
 *   - n = img.num_nonzero_pixels();          // Count pixels that are not black
 *   - rct = img.bounding_rect();             // Gradient based border removal
 *   - my_color = img.min();                  // Minimum values of all channels
 *   - my_color = img.max();                  // Maximum values of all channels
 *   - my_color = img.mean();                 // Average values of all channels
 *   - my_color = img.color_at(x,y);          // Colour of pixel x/y
 *
 * - Color conversion, channel selection:
 *
 *   - img.to_grayscale();
 *   - img.to_rgb();
 *   - img.to_hsv();
 *
 *   - img2 = img.channel_red();              // Extract channel to a new gray
 *   - img2 = img.channel_blue();             // scale image ...
 *   - img2 = img.channel_green();
 *   - img2 = img.channel_hsv_hue()
 *   - img2 = img.channel_hsv_sat();
 *   - img2 = img.channel_hsv_val();
 *
 *   - imgs = img.split_channels();           // All channels as vector of images
 *   - imgs = img.split_channels_hsv();       // All channels as HSV
 *
 * - Methods returning new images:
 *
 *   - img2 = img.subimage(rect2d(x0,y0,x1,y1));
 *   - img2 = img.in_color_range(<color>,[max_difference]);
 *
 * - Methods modifying the image and returning a reference to itself:
 *
 *   - img.negate();                          // Invert colors
 *   - img.extend(<where?>, <how?>, <color>); // Extend the image area, make border
 *   - img.crop(<where (rectangle)>);         // Crop the image
 *
 *   - img.resample_up();                     // Make it double the image size
 *   - img.resample_down();                   // Make it half the image size
 *   - img.resize_stretch(<width>,<height>,[interpolation]); // Resize, stretch
 *   - img.resize_scale(<factor>);            // .. so that the edges are factor..
 *   - img.resize(<width>,<height>);   // To min width/height, keep aspect ratio
 *
 *   - img.dilate([radius]);                  // Diliate image
 *   - img.erode([radius]);                   // Erode image
 *   - img.blur([radius], [how?]);            // Blur with selected method
 *   - img.sobel(...);                        // Sobel filter
 *   - img.scharr(...);                       // Scharr filter
 *   - img.laplacian(...);                    // Laplace derivative
 *   - img.gradient(...);                     // Gradient (uses sobel() in x/y)
 *   - img.edges_canny(...);                  // Canny edge detection
 *   - img.corners_harris(...);               // Harris corner detector
 *
 *   - img.level([black-threshold],[white-threshold]); // Colour range --> 0..MAX
 *   - img.histogram_equalize();              // Similar effect as level()
 *   - img.normalize();                       // Normalise image
 *   - img.threshold([tsh]);                  // All pixels <tsh=0, all >=tsh =MAX
 *   - img.threshold_adaptive([radius]);      // Adaptive threshold
 *
 * - Contours, feature detection:
 *
 *   - cntrs = img.contours(...);       // Find contours, return them as plain list
 *   - cntrs = img.contours_ccomp(...); // Return them as externals / holes
 *   - cntrs = img.contours_tree(...);  // Return them as tree object
 *
 *   - lines = img.hough_lines(...);          // Hough line transform
 *   - lines = img.hough_circles(...);        // Hough circle transform
 *   - hist  = img.histogram(...);            // Calculate histogram
 *
 * - Drawing
 *
 *   - img.draw( <any geometry object> );     // Draw (contours) of ...
 *   - img.fill( <any geometry object> );     // Fill area of ... with color ...
 *
 * -----------------------------------------------------------------------------
 *
 * The class templates represent own data types that can be easily extended and
 * provide converters and accessors to OpenCV types, such as cv::Mat, cv::Scalar,
 * etc. The class templates are located in namespace `sw::im`, are build around
 * the "main class" ...
 *
 *  - Template  class        : template <...> class sw::im::basic_image;
 *
 *  - Default specialisation : class sw::im::cvimage;
 *
 * ... , can have the template typename arguments ...
 *
 *  - crd_t  : Typename for coordinates, default `int`
 *  - val_t  : Typename for the color representation, default `unsigned char`
 *  - float_t: Typename for used floating point values
 *  - char_t : Typename for characters, from which strings are represented
 *
 * ... and encompass the classes:
 *
 *  - template <typename crd_t, typename val_t, typename float_t> class basic_image;
 *  - template <typename char_type> class basic_image_exception;
 *  - template <typename val_t, typename float_t, val_t max_quantum> class basic_color;
 *  - template <typename crd_t, typename float_t> class basic_point2d;
 *  - template <typename crd_t, typename float_t> class basic_line2d;
 *  - template <typename crd_t, typename float_t> class basic_rect2d;
 *  - template <typename crd_t, typename float_t> class basic_circle2d;
 *  - template <typename crd_t, typename float_t> class basic_ellipse2d;
 *  - template <typename crd_t, typename float_t> class basic_rotated_rect2d;
 *  - template <typename crd_t, typename float_t> class basic_polynom2d;
 *  - template <typename crd_t, typename float_t> class basic_contour2d;
 *  - template <typename crd_t, typename val_t, typename float_t> class basic_contour_tree;
 *  - template <typename crd_t, typename val_t, typename float_t> class basic_histogram;
 *
 * The default specialisations are derived from the default specialisation of
 * `basic_image`, which holds the fitting (template-) types for its colors,
 * geometries etc:
 *
 *  - typedef detail::basic_image<int, unsigned char, double> cvimage;
 *  - typedef cvimage::exception image_exception;
 *  - typedef cvimage::color_t color;
 *  - typedef cvimage::point_t point2d;
 *  - typedef cvimage::line_t line2d;
 *  - typedef cvimage::line_list_t lines2d;
 *  - typedef cvimage::rect_t rect2d;
 *  - typedef cvimage::circle_t circle2d;
 *  - typedef cvimage::circle_list_t circles2d;
 *  - typedef cvimage::ellipse_t ellipse2d;
 *  - typedef cvimage::rotated_rect_t rotated_rect2d;
 *  - typedef cvimage::polynom_t polynom2d;
 *  - typedef cvimage::contour_t contour2d;
 *  - typedef cvimage::contour_list_t contours2d;
 *  - typedef cvimage::contour_nested_list_t nested_contours2d;
 *
 * -----------------------------------------------------------------------------
 *
 * Usage annotations:
 *
 *  - Link against: -lopencv_core -lopencv_highgui -lopencv_imgproc
 *
 *  - `cvimage`/`basic_image` work with copies when assigning (`image2=image1;`),
 *    not with references as cv::Mat.
 *
 *  - All methods that modify the image return a reference to itself, so that
 *    the methods can be chain-called, e.g.:
 *
 *      `image1.load("file.png", cvimage::load_gray).gradient().level(0.1,1)
 *             .threshold(0.2).dilate(3).erode(1).save();` , etc etc.
 *
 *  - OpenCV matrices can be accessed and used as usual using the method
 *    `cvimavge::mat()`.
 *
 *  - `cvimage`/`basic_image` provides drawing functions for each geometry class,
 *    as well as for the histogram (so that histograms can be drawn into the image):
 *
 *    - cvimage.draw(point2d/circle2d/ellipse2d ..);
 *    - cvimage.fill(point2d/circle2d/ellipse2d ..);
 *
 *  - The colour class provides reading a colour from hex string, and the compiler
 *    implicitly converts strings to colours for method arguments:
 *    - color light_gray("#aaaaaa");             // Define a light gray color
 *    - image.fill(circle2d(x,y,r), "#ff0000");  // Draw a red filled circle
 *
 *  - Colors images (if not converted to HSV) and class `color` represent the
 *    color as RGB. Because OpenCV normally uses BGR, the "internal" CV matrices
 *    are BGR. The implicit conversion operator from `color` to `cv::Scalar`
 *    returns a BGR(A) scalar as well, and the index accessors (operator[]) of
 *    color also arranges BGR, so mycolor[0] is the blue component. Normally the
 *    inline getters `r()`, `g()`, `b()` are used to access colour components.
 *
 *  - The color space used is 8 bit, t.m. CV_8UC1 for gryscale, CV_8UC3 for RGB
 *    and HSV. If you use CV functions that convert to another matrix element
 *    data type, keep in mind that `cvimage` methods might convert it back to
 *    CV_8U.
 *
 *  - Geometry objects provide short `inline` getters/setters, type converters
 *    from/to corresponding CV data types (e.g. cv::Vec2f etc), as well as the
 *    methods `area()` and `circumference()`.
 *
 *  - Image modification methods that deal with matrix element values (pixel
 *    colours) normally take `float_t` arguments and assume that the argument
 *    you pass is a float value between 0 and 1. Hence, you do not need to care
 *    about the colour quantum range (which can be defined as `color` template
 *    value argument and is normally 0 to 255). E.g.
 *
 *     - DON'T say: img.threshold(25);  // WRONG: Threshold at 100%, as 25 >= 1
 *     - but   say: img.threshold(0.1); // RIGHT; Threshold at 10%
 *
 *
 * -----------------------------------------------------------------------------
 * +++ BSD license header +++
 * Copyright (c) 2009-2014, Stefan Wilhelm (stfwi, <cerbero s@atwilly s.de>)
 * All rights reserved.
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met: (1) Redistributions
 * of source code must retain the above copyright notice, this list of conditions
 * and the following disclaimer. (2) Redistributions in binary form must reproduce
 * the above copyright notice, this list of conditions and the following disclaimer
 * in the documentation and/or other materials provided with the distribution.
 * (3) Neither the name of atwillys.de nor the names of its contributors may be
 * used to endorse or promote products derived from this software without specific
 * prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS
 * AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
 * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
 * OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY
 * WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 * DAMAGE.
 * -----------------------------------------------------------------------------
 */
#ifndef SW_CVIMG_HH
#define	SW_CVIMG_HH

// <editor-fold desc="preprocessor" defaultstate="collapsed">
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <string>
#include <exception>
#include <vector>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <limits>
#include <cstdio>
#include <cassert>
#define dbg(X) std::cerr << X << std::endl;
#define dbg2(X) std::cerr << X << std::endl;
#define iex(S) throw exception(S) /* macro undefined at EOF */
// </editor-fold>

namespace sw { namespace im { namespace detail {

// <editor-fold desc="forward decls" defaultstate="collapsed">

template <typename crd_t, typename float_t> class basic_point2d;
template <typename crd_t, typename float_t> class basic_line2d;
template <typename crd_t, typename float_t> class basic_rect2d;
template <typename crd_t, typename float_t> class basic_circle2d;
template <typename crd_t, typename float_t> class basic_ellipse2d;
template <typename crd_t, typename float_t> class basic_rotated_rect2d;
template <typename crd_t, typename float_t> class basic_polynom2d;
template <typename crd_t, typename float_t> class basic_contour2d;
template <typename crd_t, typename val_t, typename float_t> class basic_contour_tree;

template <typename val_t=unsigned char, typename float_t=double, val_t max_quantum=255>
class basic_color;

template <typename char_type> class basic_image_exception;
template <typename crd_t, typename val_t, typename float_t> class basic_image;
template <typename crd_t, typename val_t, typename float_t> class basic_histogram;

// </editor-fold>

// <editor-fold desc="basic_image_exception" defaultstate="collapsed">
/**
 * Class basic_image_exception
 */
template <typename char_type>
class basic_image_exception : public std::exception
{
public:

  basic_image_exception() throw()
  { what_ = "<Unspecified>"; }

  explicit basic_image_exception(const char_type *wht) throw()
  { what_ = wht ? wht : "<!NULL!>"; }

  virtual ~basic_image_exception() throw()
  { ; }

  inline const char_type * what() const throw()
  { return what_; }

protected:
  const char_type *what_;
};
// </editor-fold>

// <editor-fold desc="basic_color" defaultstate="collapsed">
/**
 * Color class template
 */
template <typename val_t, typename float_t, val_t max_quantum>
class basic_color
{
public:

  /**
   * Default color: Black: RGB/HSV=0, alpha=0
   */
  inline basic_color() throw()
  { c_[0]=c_[1]=c_[2]=c_[3]=0; }

  /**
   * Gray value color: RGB/HSV=value, alpha=0
   */
  inline basic_color(val_t v) throw()
  { c_[0]=c_[1]=c_[2]=v; c_[3]=0; }

  /**
   * RGB color: Default alpha=0
   * @param hex
   */
  inline basic_color(val_t r, val_t g, val_t b, val_t a=0) throw()
  { c_[2]=r; c_[1]=g; c_[0]=b; c_[3]=a; }

  /**
   * Hex defined RGB color, alpha=255
   * @param std::string hex
   */
  inline basic_color(std::string hex)
  { operator=(hex); }

  /**
   * Hex defined RGB color, alpha=255
   * @param const char *hex
   */
  inline basic_color(const char *hex)
  { operator=(hex); }

  /**
   * cv::Scalar defined color.
   * @param cv::Scalar s
   */
  inline basic_color(cv::Scalar s) throw()
  { c_[2]=s[0]; c_[1]=s[1]; c_[0]=s[2]; c_[3]=s[3]; }

  /**
   * d'tor
   */
  virtual ~basic_color()
  { ; }

  /**
   * Assign gray scale, alpha=255
   * @param val_t c
   * @return basic_color&
   */
  inline basic_color & operator=(val_t c) throw()
  { c_[0]=c_[1]=c_[2]=c; c_[3]=0; return *this; }

  /**
   * Assign copy color
   * @param const basic_color &c
   * @return basic_color&
   */
  inline basic_color & operator=(const basic_color &c) throw()
  { c_[0]=c.c_[0]; c_[1]=c.c_[1]; c_[2]=c.c_[2]; c_[3]=c.c_[3]; return *this; }

  /**
   * Assign hex RGB, alpha=255
   * @param const char *hex
   * @return basic_color&
   */
  inline basic_color & operator=(const char *hex)
  { return operator=(std::string(hex)); }

  /**
   * RGB color hex notation e.g. "#ff0000" for red.
   * @param std::string hex
   * @return basic_color&
   */
  inline basic_color & operator=(std::string hex)
  {
    const char ct[] = "0123456789abcdef";
    c_[0]=c_[1]=c_[2]=c_[3]=0;
    #define h2d(c) ((int)( \
        ((c>='a' && c<='f') ? (c-'a'+10) : \
        ((c>='A' && c<='F') ? (c-'A'+10) : \
        ((c>='0' && c<='9') ? (c-'0') : \
        ((0)))))))
    if(hex.find_first_of("#") == 0) hex.erase(0, 1);
    if(hex.find_first_not_of(ct) != std::string::npos) return *this;

    if(hex.length() >= 6) {
      c_[0] = (h2d(hex[4])<<4) + (h2d(hex[5])<<0);  // B
      c_[1] = (h2d(hex[2])<<4) + (h2d(hex[3])<<0);  // G
      c_[2] = (h2d(hex[0])<<4) + (h2d(hex[1])<<0);  // R
      if(hex.length() < 8) return *this;
      c_[3] = (h2d(hex[6])<<4) + (h2d(hex[7])<<0);  // A
    } else if(hex.length() == 3) {
      c_[0] = (h2d(hex[2])<<4) + (h2d(hex[2])<<0);  // B
      c_[1] = (h2d(hex[1])<<4) + (h2d(hex[1])<<0);  // G
      c_[2] = (h2d(hex[0])<<4) + (h2d(hex[0])<<0);  // R
    }
    return *this;
    #undef h2d
  }

  /**
   * Returns true if the color is black (all values 0, alpha channel ignored)
   * @return bool
   */
  inline bool operator!() const throw()
  { return c_[0]==0 && c_[1]==0 && c_[2]==0; }

  /**
   * Returns true if all components of the colour are equal to the given value
   * (gray value compare). Compare ignores alpha channel.
   * @return bool
   */
  inline bool operator==(val_t c) const throw()
  { return c_[0]==c && c_[1]==c && c_[2]==c; }

  /**
   * Returns true if the colours are in all channels identical. Compare ignores alpha channel.
   * @return bool
   */
  inline bool operator==(const basic_color &c) const throw()
  { return c_[0]==c.c_[0] && c_[1]==c.c_[1] && c_[2]==c.c_[2]; }

  /**
   * Returns true if the colour identical to the given hex value.
   * @return bool
   */
  inline bool operator==(const char *hex) const
  { return *this == basic_color(std::string(hex)); }

  /**
   * Returns true if the colour identical to the given hex value. Compare ignores alpha channel.
   * @return bool
   */
  inline bool operator==(std::string hex) const
  { basic_color c(hex); return c_[0]==c.c_[0] && c_[1]==c.c_[1] && c_[2]==c.c_[2]; }

  /**
   * Returns true if all components of the colour are equal to the given value
   * (gray value compare).
   * @return bool
   */
  inline bool operator!=(val_t c) const throw()
  { return !(*this == c); }

  /**
   * Returns true if the colours are in all channels identical.
   * @return bool
   */
  inline bool operator!=(const basic_color &c) const throw()
  { return !(*this == c); }

  /**
   * Returns true if the colour identical to the given hex value.
   * @return bool
   */
  inline bool operator!=(std::string hex) const
  { return !(*this == basic_color(hex)); }

  /**
   * Returns true if the colour identical to the given hex value.
   * @return bool
   */
  inline bool operator!=(const char *hex) const
  { return !(*this == basic_color(std::string(hex))); }

  /**
   * Blue channel
   * @return val_t
   */
  inline val_t b() const throw()
  { return c_[0]; }

  /**
   * Green channel
   * @return val_t
   */
  inline val_t g() const throw()
  { return c_[1]; }

  /**
   * Red channel
   * @return val_t
   */
  inline val_t r() const throw()
  { return c_[2]; }

  /**
   * ALpha channel
   * @return val_t
   */
  inline val_t a() const throw()
  { return c_[3]; }

  /**
   * Blue channel
   * @param val_t
   */
  inline void b(val_t v) throw()
  { c_[0] = v; }

  /**
   * Green channel
   * @param val_t
   */
  inline void g(val_t v) throw()
  { c_[1] = v; }

  /**
   * Red channel
   * @param val_t
   */
  inline void r(val_t v) throw()
  { c_[2] = v; }

  /**
   * ALpha channel
   * @param val_t
   */
  inline void a(val_t v) throw()
  { c_[3] = v; }

  /**
   * Note: This is CV BGR format. color[0] is blue not red.
   * @param unsigned i
   * @return val_t
   */
  inline val_t operator[](unsigned i) const
  { return c_[i]; }

  /**
   * Note: This is CV BGR format. color[0] is blue not red.
   * @param unsigned i
   * @return val_t
   */
  inline val_t& operator[](unsigned i)
  { return c_[i]; }

  /**
   * To cv scalar
   * @return cv::Scalar
   */
  inline operator cv::Scalar () const throw()
  { return cv::Scalar(c_[0], c_[1], c_[2], c_[3]); }

  /**
   * To std::string
   * @return std::string
   */
  inline operator std::string () const
  {
    char c[16];
    ::snprintf(c,16,"#%02x%02x%02x", (int)c_[2]&0xff,(int)c_[1]&0xff,(int)c_[0]&0xff);
    std::string s(c);
    if(c_[3]) { ::snprintf(c, 16,"%02x", (int)c_[3]&0xff); }
    return s;
  }

  /**
   * Returns the maximum quantum value (for 8bit === 255 === 0xff)
   * @return const val_t
   */
  inline static const val_t max_val() throw()
  { return max_quantum; }


  /**
   * Gray scale colour, where values are between 0.0 and 1.0.
   * @param float v
   * @return basic_color
   */
  inline static basic_color grayscale(float_t v) throw()
  {
    v *= max_quantum;
    v = v < 0 ? 0 : (v > max_quantum ? max_quantum : v);
    return basic_color(v,v,v);
  }

  /**
   * RGB color from r,g,b, where values are between 0.0 and 1.0.
   * @param float r
   * @param float g
   * @param float b
   * @param float a=0
   * @return basic_color
   */
  inline static basic_color rgb(float_t r, float_t g, float_t b, float_t a=0) throw()
  {
    return basic_color(
      max_quantum * (r>1 ? 1 : (r<=0 ? 0 : r)), 255 * (g>1 ? 1 : (g<=0 ? 0 : g)),
      max_quantum * (b>1 ? 1 : (b<=0 ? 0 : b)), 255 * (a>1 ? 1 : (a<=0 ? 0 : a))
    );
  }

  /**
   * RGB color from HSV input components, where values are between 0.0 and 1.0.
   * @param float h
   * @param float s
   * @param float v
   * @return basic_color
   */
  inline static basic_color hsv(float_t h, float_t s, float_t v)
  {
    s = s<0? 0 : (s>1? 1:s);
    v = v<0? 0 : (v>1? 1:v);
    h = h-((int)h)+1; h -= ((int)h); h *= 6.0;
    int a = (int)h;
    float f = h - a;
    float p = v * (1. - s);
    float q = v * (1. - s * f);
    float t = v * (1. - s * (1.-f));
    switch(a) {
      case 0: return rgb(v,t,p);
      case 1: return rgb(q,v,p);
      case 2: return rgb(p,v,t);
      case 3: return rgb(p,q,v);
      case 4: return rgb(t,p,v);
      case 5: return rgb(v,p,q);
      default:return rgb(v,t,p);
    }
  }

  /**
   * RGB color from HSL input components, where values are between 0.0 and 1.0.
   * @param float h
   * @param float s
   * @param float l
   * @return basic_color
   */
  inline static basic_color hsl(float_t h, float_t s, float_t l)
  { s*=(l<0.5)?(l):(1-l); return hsv(h, 2.0*s/(l+s), l+s); }

  /**
   * RGB heat map color from "intensity" value between 0.0 and 1.0
   * @parem float v
   * @return basic_color
   */
  inline static basic_color heatmap(float_t v)
  { v = (v<0? 0:(v>1? 1:v)); return hsv(1-(0.2+(v*0.8)), 1, 0.25+(v*0.75)); }

public:

  /**
   * Color constants
   */
  static const basic_color white; /// #000000 alpha=0
  static const basic_color black; /// #ffffff alpha=0
  static const basic_color red;   /// #ff0000 alpha=0
  static const basic_color green; /// #00ff00 alpha=0
  static const basic_color blue;  /// #0000ff alpha=0

protected:

  val_t c_[4];
};

template <typename val_t, typename float_t, val_t max_quantum>
std::ostream& operator<<(std::ostream& o, const basic_color<val_t,float_t,max_quantum>& c)
{ o << (std::string)c; return o; }

#define defcolor(NAME, VAL) \
  template <typename val_t, typename float_t, val_t m> \
  const basic_color<val_t,float_t,m> basic_color<val_t,float_t,m>::NAME(VAL)
  defcolor(black , "000");
  defcolor(white , "fff");
  defcolor(red   , "f00");
  defcolor(green , "0f0");
  defcolor(blue  , "00f");
#undef defcolor
// </editor-fold>

// <editor-fold desc="basic_point2d" defaultstate="collapsed">
/**
 * Point class
 */
template <typename crd_t, typename float_t>
class basic_point2d
{
public:
  inline basic_point2d() throw() : x_(0), y_(0)
  {}

  inline basic_point2d(crd_t x_, crd_t y_) throw() : x_(x_), y_(y_)
  {}

  inline basic_point2d(const basic_point2d &p) throw() : x_(p.x_), y_(p.y_)
  {}

  inline basic_point2d(const cv::Point2i &p) throw() : x_(p.x), y_(p.y)
  {}

  virtual ~basic_point2d()
  { ; }

  inline crd_t x() const throw()
  { return x_; }

  inline crd_t y() const throw()
  { return y_; }

  inline basic_point2d & x(crd_t x__) throw()
  { x_=x__; return *this; }

  inline basic_point2d & y(crd_t y__) throw()
  { y_=y__; return *this; }

  inline basic_point2d& operator=(const basic_point2d& p) throw()
  { x_=p.x_; y_=p.y_; return *this; }

  inline bool operator==(const basic_point2d& p) const throw()
  { return x_==p.x_ && y_==p.y_; }

  inline bool operator>=(const basic_point2d& p) const throw()
  { return x_>=p.x_ && y_>=p.y_; }

  inline bool operator<=(const basic_point2d& p) const throw()
  { return x_<=p.x_ && y_<=p.y_; }

  inline bool operator>(const basic_point2d& p) const throw()
  { return x_>p.x_ && y_>p.y_; }

  inline bool operator<(const basic_point2d& p) const throw()
  { return x_<p.x_ && y_<p.y_; }

  inline operator cv::Point2i () const
  { return cv::Point2i(x_, y_); }

  inline operator cv::Point2f () const
  { return cv::Point2f(x_, y_); }

  inline basic_point2d operator+(const basic_point2d& p) const throw()
  { return basic_point2d(x_+p.x_, y_+p.y_); }

  inline basic_point2d operator-(const basic_point2d& p) const throw()
  { return basic_point2d(x_-p.x_, y_-p.y_); }

  inline basic_point2d& operator+=(const basic_point2d& p) throw()
  { x_+=p.x_; y_+=p.y_; return *this; }

  inline basic_point2d& operator-=(const basic_point2d& p) throw()
  { x_-=p.x(); y_-=p.y(); return *this; }

  inline basic_point2d& operator*=(double d) throw()
  { x_*=d; y_*=d; return *this; }

  inline basic_point2d& operator/=(double d)
  { x_/=d; y_/=d; return *this; }

  inline double norm() const throw()
  { return sqrt((double)x_*(double)x_ + (double)y_*(double)y_); }

protected:

  crd_t x_, y_;
};

template <typename crd_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_point2d<crd_t,float_t> p)
{ os << "{ x:" << p.x() << ", y:" << p.y() << " }"; return os; }

// </editor-fold>

// <editor-fold desc="basic_circle2d" defaultstate="collapsed">
/**
 * Circle class
 */
template <typename crd_t, typename float_t>
class basic_circle2d
{
public:
  typedef basic_point2d<crd_t,float_t> point_t;

  inline basic_circle2d() throw() : x_(0), y_(0), r_(0)
  {}

  inline basic_circle2d(crd_t x, crd_t y, crd_t r) throw() : x_(x),y_(y),r_(r)
  {}

  inline basic_circle2d(point_t p, crd_t r) throw() : x_(p.x()),y_(p.y()),r_(r)
  {}

  inline basic_circle2d(const basic_circle2d &c) throw() : x_(c.x_),y_(c.y_),r_(c.r_)
  {}

  virtual ~basic_circle2d()
  { ; }

  inline basic_circle2d& operator=(const basic_circle2d& c) throw()
  { x_=c.x_; y_=c.y_; r_=c.r_; return *this; }

  inline crd_t x() const throw()
  { return x_; }

  inline crd_t y() const throw()
  { return y_; }

  inline crd_t r() const throw()
  { return r_; }

  inline basic_circle2d & x(crd_t x) throw()
  { x_=x; return *this; }

  inline basic_circle2d & y(crd_t y) throw()
  { y_=y; return *this; }

  inline basic_circle2d & r(crd_t r) throw()
  { r_=r; return *this; }

  inline operator cv::Point () const
  { return cv::Point(x_, y_); }

  inline double area() const throw()
  { return r_ * r_ * M_PI; }

  inline basic_point2d<crd_t, float_t> center() const
  { return basic_point2d<crd_t, float_t>(x_,y_); }

  inline double circumference() const throw()
  { return  r_ * M_PI * 2.0; }

protected:
  crd_t x_, y_, r_;
};

template <typename crd_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_circle2d<crd_t,float_t> c)
{ os << "{ x:" << c.x() << ", y:" << c.y() << ", r:" << c.r() << " }"; return os; }

// </editor-fold>

// <editor-fold desc="basic_line2d" defaultstate="collapsed">
/**
 * Line class
 */
template <typename crd_t, typename float_t>
class basic_line2d
{
public:

  inline basic_line2d() throw() : p0_(), p1_()
  {}

  inline basic_line2d(crd_t x0, crd_t y0, crd_t x1, crd_t y1) throw()
    : p0_(x0,y0), p1_(x1,y1)
  {}

  inline basic_line2d(const basic_line2d &l) throw() : p0_(l.p0_), p1_(l.p1_)
  {}

  inline basic_line2d(const basic_point2d<crd_t, float_t> &p0_,
    const basic_point2d<crd_t, float_t> &p1_) throw() : p0_(p0_), p1_(p1_)
  {}

  virtual ~basic_line2d()
  { ; }

  inline basic_point2d<crd_t, float_t> p0() const throw()
  { return p0_; }

  inline basic_point2d<crd_t, float_t> p1() const throw()
  { return p1_; }

  inline basic_line2d & p0(basic_point2d<crd_t, float_t> p) throw()
  { p0_=p; return *this; }

  inline basic_line2d & p1(basic_point2d<crd_t, float_t> p) throw()
  { p1_=p; return *this; }

  inline basic_line2d& operator=(const basic_line2d& l) throw()
  { p0_=l.p0_; p1_=l.p1_; return *this; }

  inline bool operator==(const basic_line2d& l) const throw()
  { return p0_==l.p0_ && p1_==l.p1_; }

  inline double length() const throw()
  { return (p1_-p0_).norm(); }

protected:
  basic_point2d<crd_t, float_t> p0_, p1_;
};

template <typename crd_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_line2d<crd_t,float_t> l)
{ os << "{ x0:" << l.p0().x() << ", y0:" << l.p0().y()
     << ", x1:" << l.p1().x() << ", y1:" << l.p1().y()
     << " }"; return os; }

// </editor-fold>

// <editor-fold desc="basic_rect2d" defaultstate="collapsed">
/**
 * Rectangle class
 */
template <typename crd_t, typename float_t>
class basic_rect2d
{
  typedef basic_point2d<crd_t, float_t> point_t;

public:

  /**
   * Empty rect
   */
  inline basic_rect2d() throw() : p0_(), sz_()
  { ; }

  /**
   * 2 point component rect
   */
  inline basic_rect2d(crd_t x0, crd_t y0, crd_t x1, crd_t y1) throw()
  { p0_.x(std::min(x0,x1)); p0_.y(std::min(y0,y1)); sz_.x(std::abs(x1-x0)+1);
    sz_.y(std::abs(y1-y0)+1);
  }

  /**
   * Copy rect
   */
  inline basic_rect2d(const basic_rect2d &r) throw() : p0_(r.p0_), sz_(r.sz_)
  { ; }

  /**
   * 2 point rect, top-left, bottom-right
   */
  inline basic_rect2d(const point_t &p0, const point_t &p1) throw()
  { p0_.x(std::min(p0.x(),p1.x())); p0_.y(std::min(p0.y(),p1.y()));
    sz_.x(std::abs(p1.x()-p0.x())); sz_.y(std::abs(p1.y()-p0.y()));
  }

  /**
   * Point (top-left) and size components.
   */
  inline basic_rect2d(const point_t &p, crd_t width, crd_t height) throw() : p0_(p)
  { sz_.x(width<0?0:width); sz_.y(height<0?0:height); }

  /**
   * Copy cv::Rect
   */
  inline basic_rect2d(const cv::Rect r) throw()
  { p0_.x(r.x); p0_.y(r.y); sz_.x(r.width); sz_.y(r.height); }

  virtual ~basic_rect2d()
  { ; }

  inline point_t p() const throw()
  { return p0_; }

  inline point_t size() const throw()
  { return sz_; }

  inline crd_t width() const throw()
  { return sz_.x(); }

  inline void width(crd_t v) throw()
  { sz_.x(v<0?0:v); }

  inline crd_t height() const throw()
  { return sz_.y(); }

  inline void height(crd_t v) throw()
  { sz_.y(v<0?0:v); }

  inline crd_t x() const throw()
  { return p0_.x(); }

  inline void x(crd_t v) throw()
  { p0_.x(v); }

  inline crd_t y() const throw()
  { return p0_.y(); }

  inline void y(crd_t v) throw()
  { p0_.y(v); }

  inline point_t tl() const throw()
  { return p0_; }

  inline point_t br() const throw()
  { return basic_point2d<crd_t, float_t>(p0_.x()+(sz_.x()>0 ? sz_.x()-1 : 0),
              p0_.y()+(sz_.y()>0 ? sz_.y()-1 : 0)); }

  inline point_t center() const throw()
  { return basic_point2d<crd_t, float_t>(p0_.x()+sz_.x()/2, p0_.y()+sz_.y()/2); }

  inline basic_rect2d & p(point_t p) throw()
  { p0_=p; return *this; }

  inline basic_rect2d & size(point_t p) throw()
  { sz_=p; return *this; }

  inline basic_rect2d& operator=(const basic_rect2d& r) throw()
  { p0_=r.p0_; sz_=r.sz_; return *this; }

  inline bool operator==(const basic_rect2d& r) const throw()
  { return p0_ == r.p0_ && sz_ == r.sz_; }

  inline operator cv::Rect () const
  { return cv::Rect((cv::Point2i) p0_, (cv::Size2i(sz_.x(), sz_.y()))); }

  inline double area() const throw()
  { return sz_.x()*sz_.y(); }

  inline double circumference() const throw()
  { return  sz_.x() * 2.0 + sz_.y() * 2.0; }

  inline operator std::vector<cv::Point2i> () const
  {
    std::vector<cv::Point2i> v;
    v.push_back(cv::Point2i(p0_.x(), p0_.y()));
    v.push_back(cv::Point2i(p0_.x()+sz_.x(), p0_.y()));
    v.push_back(cv::Point2i(p0_.x()+sz_.x(), p0_.y()+sz_.y()));
    v.push_back(cv::Point2i(p0_.x(), p0_.y()+sz_.y()));
    v.push_back(cv::Point2i(p0_.x(), p0_.y()));
    return v;
  }

protected:

  basic_point2d<crd_t, float_t> p0_, sz_;
};

template <typename crd_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_rect2d<crd_t,float_t> l)
{ os << "{ x:" << l.x() << ", y:" << l.y() << ", w:" << l.width() << ", h:"
    << l.height() << " }"; return os; }

// </editor-fold>

// <editor-fold desc="basic_rotated_rect2d" defaultstate="collapsed">
/**
 * Rotates rectangle class
 */
template <typename crd_t, typename float_t>
class basic_rotated_rect2d
{
public:

  inline basic_rotated_rect2d() throw() : c_(), sz_(), a_(0)
  {}

  inline basic_rotated_rect2d(crd_t x, crd_t y, crd_t w, crd_t h, float_t a=0)
    throw() : c_(x,y), sz_(w,h), a_(a)
  {}

  inline basic_rotated_rect2d(const basic_rotated_rect2d &r) throw()
    : c_(r.c_), sz_(r.sz_), a_(r.a_)
  {}

  inline basic_rotated_rect2d(const basic_point2d<crd_t, float_t> &c,
    const basic_point2d<crd_t, float_t> &sz, float_t a=0) throw() : c_(c), sz_(sz),
    a_(a)
  {}

  inline basic_rotated_rect2d(const cv::RotatedRect &r) throw() : c_(r.center)
  { sz_.x(r.size.width); sz_.y(r.size.height); a_ = r.angle; }

  virtual ~basic_rotated_rect2d()
  { ; }

  inline basic_point2d<crd_t, float_t> center() const throw()
  { return c_; }

  inline basic_point2d<crd_t, float_t> size() const throw()
  { return sz_; }

  inline crd_t angle() const throw()
  { return a_; }

  inline basic_rotated_rect2d & center(const basic_point2d<crd_t, float_t>& p) throw()
  { c_ = p; return *this; }

  inline basic_rotated_rect2d & size(const basic_point2d<crd_t, float_t>& sz) throw()
  { sz_ = sz; return *this; }

  inline basic_rotated_rect2d & angle(crd_t a) const throw()
  { a_=a; return *this; }

  inline basic_rotated_rect2d& operator=(const basic_rotated_rect2d& r) throw()
  { c_=r.c_; sz_=r.sz_; a_=r.a_; return *this; }

  inline operator cv::Point () const throw()
  { return c_; }

  inline operator cv::RotatedRect () const
  { return cv::RotatedRect(cv::Point2f(c_.x(), c_.y()), cv::Size2f(sz_.x(), sz_.y()), a_); }

  inline double area() const throw()
  { return sz_.x()*sz_.y(); }

  inline double circumference() const throw()
  { return  sz_.x() * 2.0 + sz_.y() * 2.0; }

  inline operator std::vector<cv::Point2i> () const
  {
    std::vector<cv::Point2i> v;
    cv::Point2f rpt[4];
    ((cv::RotatedRect)*this).points(rpt);
    for(int i=0; i<4; ++i) v.push_back(cv::Point2i(rpt[i].x,rpt[i].y));
    return v;
  }

protected:

  basic_point2d<crd_t, float_t> c_, sz_;
  float_t a_;
};

// </editor-fold>

// <editor-fold desc="basic_ellipse2d" defaultstate="collapsed">
/**
 * Ellipse class
 */
template <typename crd_t, typename float_t>
class basic_ellipse2d
{
public:

  inline basic_ellipse2d() throw() : c_(), sz_(), a_(0)
  {}

  inline basic_ellipse2d(crd_t x, crd_t y, crd_t w, crd_t h, float_t a=0) throw()
    : c_(x,y), sz_(w,h), a_(a)
  {}

  inline basic_ellipse2d(const basic_ellipse2d &r) throw() : c_(r.c_), sz_(r.sz_),
    a_(r.a_)
  {}

  inline basic_ellipse2d(const basic_point2d<crd_t, float_t> &c, const
    basic_point2d<crd_t, float_t> &sz, float_t a=0) throw() : c_(c), sz_(sz), a_(a)
  {}

  virtual ~basic_ellipse2d()
  { ; }

  inline basic_point2d<crd_t, float_t> center() const throw()
  { return c_; }

  inline basic_point2d<crd_t, float_t> size() const throw()
  { return sz_; }

  inline crd_t angle() const throw()
  { return a_; }

  inline basic_ellipse2d & center(const basic_point2d<crd_t, float_t>& p) throw()
  { c_ = p; return *this; }

  inline basic_ellipse2d & size(const basic_point2d<crd_t, float_t>& sz) throw()
  { sz_ = sz; return *this; }

  inline basic_ellipse2d & angle(crd_t a) const throw()
  { a_=a; return *this; }

  inline basic_ellipse2d& operator=(const basic_ellipse2d& r) throw()
  { c_=r.c_; sz_=r.sz_; a_=r.a_; return *this; }

  inline operator cv::Point () const
  { return c_; }

  inline operator cv::RotatedRect () const
  { return cv::RotatedRect(cv::Point2f(c_.x(), c_.y()), cv::Size2f(sz_.x(), sz_.y()), a_); }

  inline double area() const throw()
  { return  sz_.x() * sz_.y() * M_PI/4; }

  inline double circumference() const throw()
  { return (double)2.0 * M_PI * sqrt((sz_.x()/2*sz_.x()/2)+(sz_.y()/2)*(sz_.y()/2)); }

protected:

  basic_point2d<crd_t, float_t> c_, sz_;
  float_t a_;
};
// </editor-fold>

// <editor-fold desc="basic_polynom2d" defaultstate="collapsed">
/**
 * Polynom
 */
template <typename crd_t, typename float_t>
class basic_polynom2d : public std::vector<basic_point2d<crd_t,float_t> >
{
public:

  typedef basic_point2d<crd_t,float_t> point_t;
  typedef basic_rect2d<crd_t,float_t> rect_t;
  typedef typename std::vector<point_t> vect_t;

  inline basic_polynom2d() : vect_t()
  { ; }

  inline basic_polynom2d(const point_t &p) : vect_t()
  { vect_t::push_back(p); }

  inline basic_polynom2d(const rect_t &r) : vect_t()
  {
    vect_t::push_back(point_t(0,0));
    vect_t::push_back(point_t(r.br().x(), 0));
    vect_t::push_back(point_t(r.br().x(), r.br().y() ));
    vect_t::push_back(point_t(0, r.br().y()));
    vect_t::push_back(point_t(0,0));
  }

  inline basic_polynom2d(const std::vector<cv::Point2i> &p) : vect_t()
  {
    if(p.empty()) return;
    vect_t::reserve(p.size());
    for(typename std::vector<cv::Point2i>::const_iterator it=p.begin(); it!=p.end(); ++it) {
      vect_t::push_back(point_t(*it));
    }
  }

  virtual ~basic_polynom2d()
  { ; }

  inline basic_polynom2d& add(const point_t &p)
  { vect_t::push_back(p); return *this; }

  inline basic_polynom2d& add(const basic_polynom2d &p)
  {
    if(p.empty()) return *this;
    typename vect_t::const_iterator it=p.begin();
    if(!vect_t::empty() && *(vect_t::rbegin()) == *(p.begin())) ++it;
    while(it!=p.end()) { vect_t::push_back(*it); ++it; }
    return *this;
  }

  inline basic_polynom2d& operator=(const point_t &p)
  { vect_t::clear(); push_back(p); return *this; }

  inline basic_polynom2d& operator+=(const point_t &p)
  { return add(p); }

  inline basic_polynom2d& operator+=(const basic_polynom2d &p)
  { return add(p); }

  /**
   * Returns true if the polynom is closed
   * @return
   */
  virtual bool closed() const
  {
    if(vect_t::empty()) return false;
    point_t p = *vect_t::begin(); p-=*vect_t::rbegin();
    return p.norm() < 1.5;
  }

  /**
   * To cv::Point2i vector
   * @return std::vector<cv::Point2i>
   */
  inline operator std::vector<cv::Point2i> () const
  {
    std::vector<cv::Point2i> v;
    for(typename vect_t::const_iterator it=vect_t::begin(); it!=vect_t::end(); ++it) {
      v.push_back((cv::Point2i)*it);
    }
    return v;
  }

  /**
   * Approximate lower resolution polynom with defined maximum error.
   * @param float_t max_error=1
   * @return basic_polynom2d&
   */
  inline basic_polynom2d& approx(float_t max_error=1)
  {
    std::vector<cv::Point2i> v = *this, v2;
    (*this).clear();
    cv::approxPolyDP(v, v2, max_error, closed());
    return ((*this) = basic_polynom2d(v2));
  }

  /**
   * Returns true if the polynom is convex
   * @return bool
   */
  inline bool convex() const
  { return cv::isContourConvex(cv::Mat((std::vector<cv::Point2i>)(*this))); }

  /**
   * Returns the perimeter / arc length of the polynom
   * @return double
   */
  inline double perimeter() const
  {
    if(vect_t::size() < 2) return 0;
    typename vect_t::const_iterator p0 = vect_t::begin();
    typename vect_t::const_iterator p1 = vect_t::begin()+1;
    double len = closed() ? 0 : (vect_t::front()-vect_t::back()).norm();
    for(; p1 != vect_t::end(); ++p0, ++p1) len += (*p1-*p0).norm();
    return len;
  }

  /**
   * Returns the perimeter / arc length of the polynom
   * @return double
   */
  inline double arc_length() const
  { return perimeter(); }

  /**
   * Returns the arc length of the polynom
   * @return
   */
  inline rect_t bounding_rect() const
  {
    if(vect_t::empty()) return rect_t();
    crd_t x1, y1, x0, y0;
    x1 = x0 = vect_t::front().x();
    y1 = y0 = vect_t::front().y();
    for(typename vect_t::const_iterator it = vect_t::begin(); it != vect_t::end(); ++it) {
      if(it->x() > x1) x1 = it->x(); else if(it->x() < x0) x0 = it->x();
      if(it->y() > y1) y1 = it->y(); else if(it->y() < y0) y0 = it->y();
    }
    return rect_t(x0, y0, x1, y1);
  }

};

template <typename crd_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_polynom2d<crd_t,float_t> p)
{
  if(p.empty()) return os;
  os << '[' << (*p.begin());
  for(typename basic_polynom2d<crd_t,float_t>::const_iterator it=p.begin()+1; it!=p.end(); ++it) {
    os << "," << (*it);
  }
  os << ']';
  return os;
}

// </editor-fold>

// <editor-fold desc="basic_contour2d" defaultstate="collapsed">
/**
 * Contour
 */
template <typename crd_t, typename float_t>
class basic_contour2d : public basic_polynom2d<crd_t,float_t>
{
public:
  typedef basic_point2d<crd_t,float_t> point_t;
  typedef basic_polynom2d<crd_t,float_t> poly_t;
  typedef basic_rect2d<crd_t,float_t> rect_t;
  typedef typename std::vector<point_t> vect_t;

  inline basic_contour2d() : poly_t()
  { }

  inline basic_contour2d(const point_t &p) : poly_t(p)
  { }

  inline basic_contour2d(const rect_t &r) : basic_polynom2d<crd_t,float_t>(r)
  { }

  inline basic_contour2d(const std::vector<cv::Point2i> &p) : basic_polynom2d<crd_t,float_t>(p)
  { }

  virtual ~basic_contour2d()
  { ; }

  /**
   * Returns true if the contour is closed (always true for contours not for the polygon base class)
   * @return bool
   */
  virtual bool closed() const throw()
  { return true; }

  /**
   * Returns the area of the contour
   * @param bool oriented
   * @return crd_t
   */
  inline crd_t area(bool oriented=false) const
  { return cv::contourArea(cv::Mat((std::vector<cv::Point2i>)(*this)), oriented); }

};

template <typename crd_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_contour2d<crd_t,float_t> c)
{ return operator<< (os, (basic_polynom2d<crd_t,float_t>) c); }

// </editor-fold>

// <editor-fold desc="basic_contour_tree" defaultstate="collapsed">
/**
 * Contour tree. As the class uses pointers internally, you cannot add, insert
 * or erase, only clear node contents and children using clear().
 */
template <typename crd_t, typename val_t, typename float_t>
class basic_contour_tree
{
  friend class basic_image<crd_t, val_t, float_t>;

public: // Types

  typedef basic_contour2d<crd_t, float_t> elem_t;
  typedef std::vector<basic_contour_tree> container_t;

public: // Initialisers

  inline basic_contour_tree() : elem_(), children_(), parent_(0), level_(0)
  { ; }

  inline basic_contour_tree(elem_t e) : elem_(e), children_(), parent_(0), level_(0)
  { ; }

  virtual ~basic_contour_tree()
  { clear(); }

public: // Conversions, info

  inline bool empty() const throw()
  { return children_.empty() && elem_.empty(); }

  inline bool has_children() const throw()
  { return !children_.empty(); }

  inline bool has_parent() const throw()
  { return parent_!=0; }

  inline const basic_contour_tree& parent() const
  { return parent_!=0 ? *parent_ : *this; }

  inline const elem_t& val() const throw()
  { return elem_; }

  inline elem_t& val() throw()
  { return elem_; }

  inline operator const elem_t& () const throw()
  { return elem_; }

  inline const container_t& children() const throw()
  { return children_; }

  inline operator const container_t& () const throw()
  { return children_; }

  inline void clear()
  { children_.clear(); elem_.clear(); }

protected:

  elem_t elem_;
  container_t children_;
  basic_contour_tree* parent_;
  int level_;
};
// </editor-fold>

/**
 * Class basic_image
 */
template <typename crd_t, typename val_t, typename float_t>
class basic_image
{
public:

  // <editor-fold desc="types" defaultstate="collapsed">
  typedef typename cv::Mat mat_t; // dev: can be added to template args
  typedef typename std::vector<cv::Mat> mat_vector_t;
  typedef char char_t;  // dev: can be added to template args

  typedef enum { ls_no=8, ls_connected=4, ls_antialiased=CV_AA, ls_filled=CV_FILLED } line_style_type;
  typedef enum { csp_gray=0, csp_rgb=1, csp_hsv=2 } colorspace_t;
  typedef enum { load_as_is=0, load_gray, load_color } load_mode_t;

  typedef basic_image_exception<char_t> exception;
  typedef basic_color<val_t, float_t> color_t;
  typedef basic_point2d<crd_t, float_t> point_t;
  typedef basic_line2d<crd_t, float_t> line_t;
  typedef basic_rect2d<crd_t, float_t> rect_t;
  typedef basic_circle2d<crd_t, float_t> circle_t;
  typedef basic_rotated_rect2d<crd_t, float_t> rotated_rect_t;
  typedef basic_ellipse2d<crd_t, float_t> ellipse_t;
  typedef basic_polynom2d<crd_t, float_t> polynom_t;
  typedef basic_contour2d<crd_t, float_t> contour_t;
  typedef std::vector<line_t> line_list_t;
  typedef std::vector<circle_t> circle_list_t;
  typedef std::vector<contour_t> contour_list_t;
  typedef struct { contour_list_t outer, inner; } contour_ccomp_list_t;
  typedef basic_contour_tree<crd_t,val_t,float_t> contour_nested_list_t;
  typedef basic_histogram<crd_t,val_t,float_t> histogram_t;
  // </editor-fold>

public:

  // <editor-fold desc="constructors / assignment, type conversion" defaultstate="collapsed">

  // Constructors: Cloning ensures continuous cv::Mat's

  /**
   * Standard
   */
  basic_image() throw() : m_(), file_(), is_hsv_(false)
  { ; }

  /**
   * Copy constructor, clones the cv::Mat
   * @param const basic_image &im
   */
  basic_image(const basic_image &im)  : m_(im.m_.clone()), file_(), is_hsv_(false)
  { ; }

  /**
   * Copy constructor from cv::Mat, clones the cv::Mat
   * @param const mat_t &m
   */
  basic_image(const mat_t &m) : m_(m.clone()), file_(), is_hsv_(false)
  { ; }

  /**
   * Constructor by size and color space
   * @param crd_t width
   * @param crd_t height
   * @param colorspace_t colorsp
   */
  inline basic_image(crd_t width, crd_t height, colorspace_t colorsp) : m_(), file_(),
          is_hsv_(colorsp==csp_hsv)
  { m_ = cv::Mat::zeros((int)height, (int)width, colorsp == csp_gray ? CV_8UC1 : CV_8UC3); }

  /**
   * Constructor by size and color space and predefined "background color".
   * @param crd_t width
   * @param crd_t height
   * @param colorspace_t colorsp
   */
  basic_image(crd_t width, crd_t height, colorspace_t colorsp, const color_t & color) :
    m_(), file_(), is_hsv_(colorsp==csp_hsv)
  {
    m_.create((int)height, (int)width, colorsp == csp_gray ? CV_8UC1 : CV_8UC3);
    m_ = (cv::Scalar) color;
  }

  /**
   * Constructor by file
   * @param const std::string& file
   */
  inline basic_image(const std::string& file, load_mode_t mode=load_as_is) : m_(), file_(),
          is_hsv_(false)
  { load(file, mode); }

  virtual ~basic_image()
  { ; }

  inline basic_image& operator=(const basic_image& im)
  { m_ = im.mat().clone(); is_hsv_ = im.is_hsv_; return *this; }

  inline basic_image& operator=(const mat_t& m)
  { m_ = m.clone(); is_hsv_=false; return *this; }

  inline mat_t& mat() throw()
  { return m_; }

  inline const mat_t& mat() const throw()
  { return m_; }

  inline operator mat_t& () throw()
  { return m_; }

  inline operator const mat_t& () throw()
  { return m_; }

  inline basic_image& mat(const mat_t& m)
  { m_ = m.clone(); return *this; }

  // </editor-fold>

public: // Import / export

  // <editor-fold desc="load/save" defaultstate="collapsed">

  /**
   * Load an image
   * @param const std::string& file
   * @param load_mode mode=load_as_is
   * @return basic_image&
   */
  inline basic_image& load(const std::string& file, load_mode_t mode=load_as_is)
  {
    try {
      file_ = "";
      m_ = mat_t();
      m_ = cv::imread(file, (
        mode == load_as_is ? CV_LOAD_IMAGE_UNCHANGED : (
        mode == load_gray ? CV_LOAD_IMAGE_GRAYSCALE :
        CV_LOAD_IMAGE_COLOR
      )));
      file_ = file;
      is_hsv_ = false;
    } catch(...) {
      iex("Failed to load image.");
    }
    return *this;
  }

  /**
   * Returns the file the image was loaded from.
   * @param file const std::string&
   * @return
   */
  inline const std::string& file() const throw()
  { return file_; }

  /**
   * Save the image as RGB or gray scale, depending on the image with the file
   * name that it was loaded with.
   * @return basic_image&
   */
  inline basic_image& save()
  { save(file_); return *this; }

  /**
   * Save the image as RGB or gray scale, depending on the image. Does not modify the
   * loaded file path retrieved with file(), except this path is empty (image is not a loaded
   * file).
   * @param const std::string &file
   * @return basic_image&
   */
  inline basic_image& save(const std::string &file)
  {
    try{
      if(!empty()) { if(is_hsv_) to_rgb(); cv::imwrite(file, m_); }
      if(file_.empty()) file_ = file;
    } catch(...) {
      iex("Failed to save image.");
    }
    return *this;
  }

  // </editor-fold>

public: // Image information

  // <editor-fold desc="general getters" defaultstate="collapsed">
  /**
   * Return maximum quantum value
   * @return double
   */
  inline const val_t max_quant() const throw()
  { return color_t::max_val(); }  // CHECK more common way: _m.elemSize1() or the like

  /**
   * Returns if the image is empty.
   * @return
   */
  inline bool empty() const throw()
  { return !m_.data || m_.cols == 0 || m_.rows == 0; }

  /**
   * Returns the number of image channels.
   * @return unsigned
   */
  inline unsigned channels() const throw()
  { return m_.channels(); }

  /**
   * Image width
   * @return crd_t
   */
  inline crd_t width() const throw()
  { return m_.cols; }

  /**
   * Image height
   * @return crd_t
   */
  inline crd_t height() const throw()
  { return m_.rows; }

  /**
   * Image area
   * @return crd_t
   */
  inline crd_t area() const throw()
  { return m_.rows * m_.cols; }

  /**
   * Image size rectangle
   * @return rect_t
   */
  inline point_t size() const throw()
  { return point_t(width(), height()); }

  /**
   * Returns if the image color space is currently RGB (t.m. matrix is BGR, 3ch)
   * @return bool
   */
  inline bool is_rgb() const throw()
  { return !is_hsv_ && channels()>1; }

  /**
   * Returns if the image color space is currently gray scale (1 channel).
   * @return bool
   */
  inline bool is_gray() const throw()
  { return channels()<2; }

  /**
   * Returns if the image color space is currently HSV (3 channels)
   * @return bool
   */
  inline bool is_hsv() const throw()
  { return is_hsv_ && channels()>1; }

  /**
   * Returns the current color space of the image
   * @return colorspace_t
   */
  inline colorspace_t color_space() const throw()
  { return (channels()<2) ? csp_gray : (is_hsv_ ? csp_hsv : csp_rgb); }

  /**
   * Sets the current color space of the image
   * @param colorspace_t
   */
  inline basic_image& color_space(colorspace_t c)
  {
    if(color_space() == c) return *this;
    else if(c == csp_gray) return to_grayscale();
    else if(c == csp_rgb) return to_rgb();
    else if(c == csp_hsv) return to_hsv();
    return *this;
  }

  // </editor-fold>

  // <editor-fold desc="num_nonzero_pixels, nonzero_pixels" defaultstate="collapsed">
  /**
   * Returns the number of nonzero pixels. If the image is not gray scale,
   * a temporary gray scale copy is created (performance!).
   * @return long
   */
  inline long num_nonzero_pixels() const
  {
    return (channels()<2) ? cv::countNonZero(m_) :
      basic_image(*this).to_grayscale().num_nonzero_pixels();
  }

 /**
   * Return nonzero pixel coordinates.
   * Note: E.g.
   *
   *    struct my_point_type; // MUST have constructor my_point_type(x, y)
   *    std::list<my_point_type> pts;
   *    unsigned n = img.unsigned nonzero_pixels(pts);
   *
   * @param typename point_list_type& points
   * @return unsigned
   */
  template <typename point_list_type>
  inline unsigned nonzero_pixels(point_list_type& points) const
  {
    points.clear();
    basic_image im = *this;
    im.to_grayscale();
    unsigned np = im.num_nonzero_pixels();
    if(!np) return 0;
    register unsigned n = np;
    for(register int iy=im.height()-1; iy >= 0; --iy)  {
      for(register int ix=im.width()-1; ix >= 0; --ix)  {
        if(im.gray_value_at(ix, iy)) {
          points.push_back(typename point_list_type::value_type(ix,iy));
          if(--n <= 0) return np;
        }
      }
    }
    return np;
  }

 /**
   * Return nonzero pixel coordinates. Specialised version with `reserve()` call.
   * Note: E.g.
   *
   *    struct my_point_type; // MUST have constructor my_point_type(x, y)
   *    std::vector<my_point_type> pts;
   *    unsigned n = img.unsigned nonzero_pixels(pts);
   *
   * @param std::vector<point_type>& points
   * @return unsigned
   */
  template <typename point_type>
  inline unsigned nonzero_pixels(std::vector<point_type>& points) const
  {
    points.clear();
    basic_image im = *this;
    im.to_grayscale();
    unsigned np = im.num_nonzero_pixels();
    if(!np) return 0;
    points.reserve(np);
    register unsigned n = np;
    for(register int iy=im.height()-1; iy >= 0; --iy)  {
      for(register int ix=im.width()-1; ix >= 0; --ix)  {
        if(im.gray_value_at(ix, iy)) {
          points.push_back(point_type(ix,iy));
          if(--n <= 0) return np;
        }
      }
    }
    return np;
  }
  // </editor-fold>

  // <editor-fold desc="bounding_rect" defaultstate="collapsed">
  /**
   * Determines the bounding rectangle of the image. T.m. the image content
   * rectangle without the border. Uses x/y Sobel gradient.
   * @param float_t gradient_threshold=0.2
   * @param int gradient_radius=0
   * @return rect_t
   */
  inline rect_t bounding_rect(float_t gradient_threshold, int gradient_radius=0) const
  {
    // +++ check if finding external contours is faster +++
    basic_image im(*this);
    im.gradient(gradient_radius).threshold(gradient_threshold);
    int x0=0, x1=im.width()-1, y0=0, y1=im.height()-1;
    while(y1>0 && cv::sum(im.mat().row(y1))[0]==0) y1--;
    while(y0<y1 && cv::sum(im.mat().row(y0))[0]==0) y0++;
    while(x1>0 && cv::sum(im.mat().col(x1))[0]==0) x1--;
    while(x0<x1 && cv::sum(im.mat().col(x0))[0]==0) x0++;
    return rect_t(x0, y0, x1, y1);
  }

  /**
   * Determines the bounding rectangle of the image. All nonzero pixels count.
   * @return rect_t
   */
  inline rect_t bounding_rect() const
  {
    basic_image im(*this);
    im.threshold(0.1);
    int x0=0, x1=im.width()-1, y0=0, y1=im.height()-1;
    while(y1>0 && cv::sum(im.mat().row(y1))[0]==0) y1--;
    while(y0<y1 && cv::sum(im.mat().row(y0))[0]==0) y0++;
    while(x1>0 && cv::sum(im.mat().col(x1))[0]==0) x1--;
    while(x0<x1 && cv::sum(im.mat().col(x0))[0]==0) x0++;
    return rect_t(x0, y0, x1, y1);
  }
  // </editor-fold>

  // <editor-fold desc="min, max, mean" defaultstate="collapsed">
  /**
   * Returns the minimum value in the image matrix (each channel separately).
   * @return val_t
   */
  inline color_t min() const
  {
    if(is_gray()) {
      double d; cv::minMaxLoc(m_, &d, 0); return color_t((int)d);
    } else {
      color_t c; double d;
      mat_vector_t chs = split_channels();
      for(unsigned i=0; i<3&&i<chs.size(); ++i) {
        cv::minMaxLoc(chs[i], &d, 0); c[i] = d;
      }
      return c;
    }
  }

  /**
   * Returns the minimum value in the image matrix (each channel separately).
   * @return val_t
   */
  inline color_t max() const
  {
    if(is_gray()) {
      double d; cv::minMaxLoc(m_, 0, &d); return color_t((int)d);
    } else {
      color_t c; double d;
      mat_vector_t chs = split_channels();
      for(unsigned i=0; i<3&&i<chs.size(); ++i) {
        cv::minMaxLoc(chs[i], 0, &d); c[i] = d;
      }
      return c;
    }
  }

  /**
   * Returns the mean "colour" of all image pixels. If the image is gray scale,
   * the all R,G,B have the same value.
   * @return color_t
   */
  inline color_t mean() const
  { return is_gray() ? color_t(cv::mean(m_)[0]) : color_t(cv::mean(m_)); }

  /**
   * Returns the mean "colour" of all image pixels. If the image is gray scale,
   * the all R,G,B have the same value.
   * @return color_t
   */
  inline color_t mean(const basic_image& mask) const
  {
    cv::Mat msk;
    if(!mask.empty()) msk = mask.mat();
    return is_gray() ? color_t(cv::mean(m_, msk)[0]) : color_t(cv::mean(m_, msk));
  }
  // </editor-fold>

  // <editor-fold desc="color_at, color_mean" defaultstate="collapsed">
  /**
   * Returns colour at given coordinates.
   * NOTE: This method is NOT suitable for large scale iterations, as it needs
   * to determine the matrix channels/depth every time it is called (switch()).
   * @param crd_t x
   * @param crd_t y
   * @return color_t
   */
  color_t color_at(crd_t x, crd_t y) const
  {
    color_t c;
    if(channels()!=1 && channels() != 3) {
      iex("color_at(): Unsupported number of channels (must be 1 or 3).");
    } else if(x<0 || y<0 || x>=width() || y>=height()) {
      return color_t();
    }
    switch(m_.depth()) {
      case CV_8U:
        if(channels()==1) {
          c = color_t(m_.at<unsigned char>(y,x));
        } else {
          cv::Vec3b v(m_.at<cv::Vec3b>(y,x));
          c = color_t(v[2],v[1],v[0]);
        }
        break;
      case CV_32F:
        if(channels()==1) {
          c = color_t((unsigned)m_.at<float>(y,x));
        } else {
          cv::Vec3f v(m_.at<cv::Vec3f>(y,x));
          c = color_t(v[2],v[1],v[0]);
        }
        break;
      case CV_64F:
        if(channels()==1) {
          c = color_t((unsigned)m_.at<double>(y,x));
        } else {
          cv::Vec3d v(m_.at<cv::Vec3d>(y,x));
          c = color_t(v[2],v[1],v[0]);
        }
        break;
      default:
        iex("color_at(): Unsupported depth.");
    }
    return c;
  }

 /**
   * Sets the color at a given pixel.
   * NOTE: This method is NOT suitable for large scale iterations, as it needs
   * to determine the matrix channels/depth every time it is called (switch()).
   * @param crd_t x
   * @param crd_t y
   * @param const color_t& c
   * @return color_t
   */
  basic_image& color_at(crd_t x, crd_t y, const color_t& c)
  {
    if(channels()!=1 && channels() != 3) {
      iex("color_at(): Unsupported number of channels (must be 1 or 3).");
    } else if(x<0 || y<0 || x>=width() || y>=height()) {
      return *this;
    }
    switch(m_.depth()) {
      case CV_8U:
        if(channels()==1) {
          m_.at<unsigned char>(y,x) = (((int)c[0]) + c[1] + c[2])/3;
        } else {
          m_.at<cv::Vec3b>(y,x) = cv::Vec3b(c[2], c[1], c[0]);
        }
        break;
      case CV_32F:
        if(channels()==1) {
          m_.at<float>(y,x) = (((float)c[0]) + c[1] + c[2])/3;
        } else {
          m_.at<cv::Vec3f>(y,x) = cv::Vec3f(c[2], c[1], c[0]);
        }
        break;
      case CV_64F:
        if(channels()==1) {
          m_.at<double>(y,x) = (((double)c[0]) + c[1] + c[2])/3;
        } else {
          m_.at<cv::Vec3d>(y,x) = cv::Vec3d(c[2], c[1], c[0]);
        }
        break;
      default:
        iex("color_at(): Unsupported depth.");
    }
    return *this;
  }

 /**
   * Sets the gray scale value at a given pixel.
   * NOTE: This method is NOT suitable for large scale iterations, as it needs
   * to determine the matrix channels/depth every time it is called (switch()).
   * @param crd_t x
   * @param crd_t y
   * @param val_t c
   * @return color_t
   */
  basic_image& color_at(crd_t x, crd_t y, val_t c)
  {
    if(channels()!=1 && channels() != 3) {
      iex("color_at(): Unsupported number of channels (must be 1 or 3).");
    } else if(x<0 || y<0 || x>=width() || y>=height()) {
      return *this;
    }
    switch(m_.depth()) {
      case CV_8U:
        if(channels()==1) {
          m_.at<unsigned char>(y,x) = c;
        } else {
          m_.at<cv::Vec3b>(y,x) = cv::Vec3b(c,c,c);
        }
        break;
      case CV_32F:
        if(channels()==1) {
          m_.at<float>(y,x) = c;
        } else {
          m_.at<cv::Vec3f>(y,x) = cv::Vec3f(c,c,c);
        }
        break;
      case CV_64F:
        if(channels()==1) {
          m_.at<double>(y,x) = c;
        } else {
          m_.at<cv::Vec3d>(y,x) = cv::Vec3d(c,c,c);
        }
        break;
      default:
        iex("color_at(): Unsupported depth.");
    }
    return *this;
  }

  /**
   * Returns colour at given point coordinates
   * @param const point_t &p
   * @return color_t
   */
  inline color_t color_at(const point_t &p) const
  { return color_at(p.x(), p.y()); }

 /**
   * Returns the gray value at a given pixel or -1 if the image is not gray.
   * @param crd_t x
   * @param crd_t y
   * @return int
   */
  inline int gray_value_at(crd_t x, crd_t y) const
  {
    #if defined(DEBUG) || defined(_DEBUG) || defined(_DEBUG_)
    return (channels() !=1 || m_.depth() != CV_8U) ? -1 : m_.at<unsigned char>(y,x);
    #else
    return m_.at<unsigned char>(y,x);
    #endif
  }

  // </editor-fold>

public: // Image (matrix) conversions, channel selection

  // <editor-fold desc="color conversion" defaultstate="collapsed">
  /**
   * Converts the image to gray scale (if not yet gray scale)
   * @return basic_image&
   */
  inline basic_image& to_grayscale()
  {
    if(m_.type() != CV_8U) m_.convertTo(m_, CV_8U);
    if(channels() > 1) {
      if(is_hsv_) cv::cvtColor(m_, m_, CV_HSV2BGR);
      cv::cvtColor(m_, m_, CV_BGR2GRAY);
    }
    is_hsv_ = false;
    return *this;
  }

  /**
   * Converts the image to BGR, 8 bits per channel (if not yet gray scale)
   * @return basic_image&
   */
  inline basic_image& to_rgb()
  {
    if(m_.type() != CV_8U) m_.convertTo(m_, CV_8U);
    if(channels() > 1 && is_hsv_) cv::cvtColor(m_, m_, CV_HSV2BGR);
    if(channels() == 1) cv::cvtColor(m_, m_, CV_GRAY2BGR);
    is_hsv_ = false;
    return *this;
  }

  /**
   * Converts image to HSV, 8 bits per channel
   * @return basic_image&
   */
  inline basic_image& to_hsv()
  {
    if(m_.type() != CV_8U) m_.convertTo(m_, CV_8U);
    if(channels() == 1) cv::cvtColor(m_, m_, CV_GRAY2BGR);
    if(is_hsv_) return *this;
    cv::cvtColor(m_, m_, CV_BGR2HSV);
    is_hsv_ = true;
    return *this;
  }
  // </editor-fold>

  // <editor-fold desc="channel extraction / composition" defaultstate="collapsed">
  /**
   * Splits image into the three RGB channels, each one gray scale image.
   * @return mat_vector;
   */
  inline mat_vector_t split_channels() const
  { mat_vector_t v; cv::split(m_, v); return v; }

  /**
   * Merges channels in to this image.
   * @param const mat_vector_t& channels;
   * @return basic_image& (*this)
   */
  inline basic_image& merge_channels(const mat_vector_t& chs)
  { cv::merge(chs, m_); return *this; }

  /**
   * Splits image into the three HSV channels, each one gray scale image.
   * @return mat_vector;
   */
  inline mat_vector_t split_channels_hsv() const
  {
    basic_image im = *this;
    im.to_rgb();
    cv::cvtColor(im.mat(), im.mat(), CV_BGR2HSV);
    mat_vector_t v;
    cv::split(im.mat(), v);
    return v;
  }

  /**
   * Returns RGB B channel as gray scale image
   * @return basic_image
   */
  inline basic_image channel_blue() const
  { return (channels()==1) ? basic_image(m_) : basic_image(split_channels()[0]); }

  /**
   * Returns RGB G channel as gray scale image
   * @return basic_image
   */
  inline basic_image channel_green() const
  { return (channels()==1) ? basic_image(m_) : basic_image(split_channels()[1]); }

  /**
   * Returns RGB R channel as gray scale image
   * @return basic_image
   */
  inline basic_image channel_red() const
  { return (channels()==1) ? basic_image(m_) : basic_image(split_channels()[2]); }

  /**
   * Returns BGR or HSV H channel as gray scale image
   * @return basic_image
   */
  inline basic_image channel_hsv_hue() const
  { return is_hsv_ ? split_channels()[0] : split_channels_hsv()[0]; }

  /**
   * Returns BGR or HSV S channel as gray scale image
   * @return basic_image
   */
  inline basic_image channel_hsv_sat() const
  { return is_hsv_ ?  split_channels()[1] : split_channels_hsv()[1]; }

  /**
   * Returns BGR or HSV V channel as gray scale image
   * @return basic_image
   */
  inline basic_image channel_hsv_val() const
  { return is_hsv_ ?  split_channels()[2] : split_channels_hsv()[2]; }

  // </editor-fold>

public: // Methods returning a copy / new image

  // <editor-fold desc="subimage, in_color_range" defaultstate="collapsed">
  /**
   * Returns the sub image (a COPY) specified by the rectangle `r`.
   * @param const rect_t &r
   * @return basic_image
   */
  inline basic_image subimage(const rect_t &r) const
  { return m_((cv::Rect)r).clone(); }

  /**
   * Creates a mask where all pixels that are in range of the specified colour
   * are nonzero, all other colours are black (0). The mask is a CLONED image.
   * The difference is normed to 0..1 (for 0..max_quant()). Geometric color
   * difference is basically the magnitude of the colour interpreted as 3d vector
   * compared, instead of simple cv::inRange(). For gray scale images,
   * geometric colour difference is not applied.
   * @param const color_t &c
   * @param float_t max_difference=0.1
   * @param bool geometric_color_difference=false
   * @return basic_image
   */
  basic_image in_color_range(const color_t &c, float_t max_difference=0.1,
      bool geometric_color_difference=false)
  {
    double d = max_difference * max_quant();
    if(channels()==1) {
      mat_t m(m_.rows, m_.cols, CV_8U);
      cv::inRange(m_, cv::Scalar(c.r()-d), cv::Scalar(c.r()+d), m);
      return m;
    } else if(!geometric_color_difference) {
      mat_t m(m_.rows, m_.cols, CV_8U);
      cv::inRange(m_,
        cv::Scalar(c.b()-d, c.g()-d, c.r()-d),
        cv::Scalar(c.b()+d, c.g()+d, c.r()+d),
      m);
      return m;
    } else {
      cv::Mat3b m;
      {
        basic_image im(*this);
        im.to_rgb(); // Ensure BGR uint8_t, continuous
        m = im.mat();
      }
      cv::Mat1b r(m.rows, m.cols);
      cv::Mat1b::iterator rit = r.begin();
      for(cv::Mat3b::const_iterator it=m.begin(); it != m.end(); ++it, ++rit) {
        long cb=(*it)[0]-c.b(), cg=(*it)[1]-c.g(), cr=(*it)[2]-c.r();
        *rit = sqrt(cb*cb+cg*cg+cr*cr) <= d ? 255 : 0;
      }
      return r;
    }
  }
  // </editor-fold>

public: // Methods modifying the image (and returning a reference to itself))

  // <editor-fold desc="clear" defaultstate="collapsed">
  /**
   * Clears the contents and resets all instance variables to defaults.
   * @return basic_image&
   */
  inline basic_image& clear()
  { m_ = cv::Mat(); is_hsv_ = false; file_ = ""; return *this; }

  // </editor-fold>

  // <editor-fold desc="negate/invert" defaultstate="collapsed">

  /**
   * Negate the image ("invert")
   * @return basic_image&
   */
  inline basic_image& negate()
  {
    if(m_.channels() == 1) {
      m_ = max_quant() - m_;
    } else {
      mat_vector_t im;
      cv::split(m_, im);
      for(typename mat_vector_t::iterator it=im.begin(); it!=im.end(); ++it) {
        (*it) = max_quant() - (*it);
      }
      cv::merge(im, m_);
      im.clear();
    }
    return *this;
  }

  /**
   * Negate the image ("invert")
   * @return basic_image&
   */
  inline basic_image& invert()
  {
    return negate();
  }

  // </editor-fold>

  // <editor-fold desc="extend, crop" defaultstate="collapsed">

  typedef enum { em_color=0, em_replicate, em_reflect, em_wrap } extend_mode_t;

  /**
   * Extend an image area
   * @param const point_t &top_left
   * @param const point_t &bottom_right
   * @param extend_mode_t how=em_color
   * @param const color_t &c=color_t()
   * @return basic_image&
   */
  inline basic_image& extend(const point_t &top_left, const point_t &bottom_right,
      extend_mode_t how=em_color, const color_t &c=color_t())
  {
    int t;
    switch(how) {
      case em_replicate: t=cv::BORDER_REPLICATE; break;
      case em_reflect: t=cv::BORDER_REFLECT; break;
      case em_wrap: t=cv::BORDER_WRAP; break;
      case em_color: default: t=cv::BORDER_CONSTANT;
    }
    cv::copyMakeBorder(m_, m_, top_left.y(), bottom_right.y(), top_left.x(),
        bottom_right.x(), t, (cv::Scalar) c);
    return *this;
  }

  /**
   * Extend an image area
   * @param crd_t px
   * @param extend_mode_t how=em_color
   * @param const color_t &c=color_t()
   * @return basic_image&
   */
  inline basic_image& extend(crd_t px, extend_mode_t how=em_color,
      const color_t &c=color_t())
  { return extend(point_t(px,px), point_t(px,px), how, c); }


  /**
   * Extend to a specified width and height. If `centered==true`, border at top,
   * left, bottom, right is added; otherwise only bottom and right.
   * @param crd_t width
   * @param crd_t height
   * @param bool centered=false
   * @param extend_mode_t how=em_color
   * @param const color_t &c=color_t()
   * @return basic_image&
   */
  inline basic_image& extend(crd_t width, crd_t height, bool centered=false,
      extend_mode_t how=em_color, const color_t &c=color_t())
  {
    if(width<this->width() && height<this->height()) return *this;
    crd_t dw = width < this->width() ? 0 : width-(this->width());
    crd_t dh = height < this->height() ? 0 : height-(this->height());
    crd_t t, b, l, r;
    if(centered) {
      l = r = dw/2; r += (l+r < dw) ? 1 : 0; // rounding
      t = b = dh/2; b += (t+b < dh) ? 1 : 0;
    } else {
      t = l = 0; r = dw; b = dh;
    }
    dbg2("t=" << t << ", l=" << l << ", b=" << b << ", r=" << r);
    return extend(point_t(t,l), point_t(b,r), how, c);
  }

  /**
   * Crop the image to the specified region
   * @param const rect_t& region
   * @return basic_image&
   */
  inline basic_image& crop(const rect_t& region)
  {
    crd_t x0, y0, x1, y1;
    if(region.width()<1 || region.height()<1) { clear(); return *this; }
    x0 = region.tl().x(); y0 = region.tl().y();
    x1 = region.br().x(); y1 = region.br().y();
    x0 = x0<0 ? 0 : x0;
    y0 = y0<0 ? 0 : y0;
    x1 = x1>=width() ? width()-1 : x1;
    y1 = y1>=height() ? height()-1 : y1;
    rect_t r(x0, y0, x1, y1);
    if(r.width() >= width() && r.height() > height()) return *this;
    m_ = m_((cv::Rect)r).clone(); // Re-assign to ensure continuous mat.
    return *this;
  }

  /**
   * Crop the image to the specified width/height from the the bottom-right
   * corner.
   * @param crd_t width
   * @param crd_t height
   * @return basic_image&
   */
  inline basic_image& crop(crd_t width, crd_t height)
  { return crop(rect_t(0, 0, width, height)); }

  // </editor-fold>

  // <editor-fold desc="resample, stretch, scale, resize" defaultstate="collapsed">

  /**
   * Resamples the image, doubles the width and height.
   * @return basic_image&
   */
  inline basic_image& resample_up()
  { cv::pyrUp(m_, m_); return *this; }

  /**
   * Resamples the image, halves the width and height.
   * @return basic_image&
   */
  inline basic_image& resample_down()
  { cv::pyrDown(m_, m_); return *this; }


  /**
   * Resizing interpolation methods
   * @enum resize_interpolation_t
   */
  typedef enum {
    resz_bilinear = cv::INTER_LINEAR,
    resz_nearest = cv::INTER_NEAREST,
    resz_cubic = cv::INTER_CUBIC,
    resz_area = cv::INTER_AREA,
    resz_laczos4 = cv::INTER_LANCZOS4
  } resize_interpolation_t;

  /**
   * Resizes the image to (width)x(height) without keeping the aspect ratio.
   * @param crd_t width
   * @param crd_t height
   * @param resize_interpolation_t interpolation
   * @return basic_image&
   */
  inline basic_image& resize_stretch(crd_t width, crd_t height, resize_interpolation_t
    interpolation=resz_area)
  {
    if(width<=0 || height<=0) { m_ = cv::Mat(0,0,m_.type()); return *this; }
    cv::resize(m_, m_, cv::Size(width, height), 0, 0, (int)interpolation);
    return *this;
  }

  /**
   * Resizes the image using a scale factor (>0). The factor refers to the edge lengths.
   * @param float_t factor
   * @return basic_image&
   */
  inline basic_image& resize_scale(float_t factor)
  {
    if(factor<=0) return *this;
    resize_stretch(::round(factor*width()), ::round(factor*height()),
        factor<1 ? resz_area : resz_nearest);
    return *this;
  }

  /**
   * Resizes the image to a maximum height/width. Keeps aspect ratio. Scales
   * up and down.
   * @param crd_t max_width
   * @param crd_t max_height
   * @return basic_image&
   */
  inline basic_image& resize(crd_t max_width, crd_t max_height)
  {
    float_t fw=(float_t)max_width/width(), fh=(float_t)max_height/height();
    if(fw < fh) {
      return resize_stretch(max_width, fw*height());
    } else {
      return resize_stretch(fh*width(), max_height);
    }
  }

  // </editor-fold>

  // <editor-fold desc="level, normalize, histogram_equalize" defaultstate="collapsed">
  /**
   * Linear stretch the colour value from min->max ---> 0->max_quant()
   * @return basic_image&
   */
  inline basic_image& level()
  {
    if(m_.channels() == 1) {
      double min, max;
      cv::minMaxLoc(m_, &min, &max);
      if(max > min) { m_ -= min; m_ *= (max_quant() / (max-min)); }
    } else {
      mat_vector_t im;
      cv::split(m_, im);
      for(typename mat_vector_t::iterator it=im.begin(); it!=im.end(); ++it) {
        double min, max;
        cv::minMaxLoc((*it), &min, &max);
        if(max > min) { (*it) -= min; (*it) *= (max_quant() / (max-min)); }
      }
      cv::merge(im, m_);
      im.clear();
    }
    return *this;
  }

  /**
   * Linear stretch the colour value from min->max ---> 0->max_quant(),
   * where the image colour range is cropped to min_0_1, max_0_1.
   * @param double min_0_1
   * @param double max_0_1
   * @return basic_image&
   */
  basic_image& level(double min_0_1, double max_0_1)
  {
    min_0_1 = min_0_1 > 1 ? 1 : min_0_1 < 0 ? 0 : min_0_1;
    max_0_1 = max_0_1 > 1 ? 1 : max_0_1 < 0 ? 0 : max_0_1;
    if(channels() == 1) {
      m_.convertTo(m_, CV_8U);
      cv::threshold(m_, m_, min_0_1*max_quant(), max_quant(), cv::THRESH_TOZERO);
      cv::threshold(m_, m_, max_0_1*max_quant(), max_quant(), cv::THRESH_TRUNC);
      double min, max;
      cv::minMaxLoc(m_, &min, &max);
      if(max > min) { m_ -= min; m_ *= (max_quant() / (max-min)); }
    } else {
      mat_vector_t im;
      cv::split(m_, im);
      for(typename mat_vector_t::iterator it=im.begin(); it!=im.end(); ++it) {
        (*it).convertTo((*it), CV_8U);
        cv::threshold((*it), (*it), min_0_1*max_quant(), max_quant(), cv::THRESH_TOZERO);
        cv::threshold((*it), (*it), max_0_1*max_quant(), max_quant(), cv::THRESH_TRUNC);
        double min, max;
        cv::minMaxLoc((*it), &min, &max);
        if(max > min) { (*it) -= min; (*it) *= (max_quant() / (max-min)); }
      }
      cv::merge(im, m_);
      im.clear();
    }
    return *this;
  }

  /**
   * Normalise image to values between min (0 to 1.0) and max (0 to 1.0).
   * @param double min_0_1
   * @param double max_0_1
   * @return basic_image&
   */
  inline basic_image& normalize(double min_0_1, double max_0_1)
  {
    min_0_1 = min_0_1 > 1 ? 1 : min_0_1 < 0 ? 0 : min_0_1;
    max_0_1 = max_0_1 > 1 ? 1 : max_0_1 < 0 ? 0 : max_0_1;
    cv::normalize(m_, m_, min_0_1*max_quant(), max_0_1*max_quant(), cv::NORM_MINMAX);
    return *this;
  }

  /**
   * Equalises the image histogram (for each channel independently)
   * @return basic_image&
   */
  inline basic_image& histogram_equalize()
  {
    if(channels()==1) {
      cv::equalizeHist(m_, m_);
    } else {
      mat_vector_t chs = split_channels();
      for(unsigned i=0; i<chs.size(); ++i)
        cv::equalizeHist(chs[i], chs[i]);
      cv::merge(chs, m_);
    }
    return *this;
  }

  // </editor-fold>

  // <editor-fold desc="threshold, threshold_adaptive" defaultstate="collapsed">
  /**
   * Sets all pixels with a value smaller than `value_0_1*max_quant()` to 0,
   * all other pixels to max_quant().
   * @param double value_0_1
   * @return basic_image&
   */
  inline basic_image& threshold(double value_0_1=.5)
  {
    value_0_1 = value_0_1 > 1 ? 1 : value_0_1 < 0 ? 0 : value_0_1;
    to_grayscale();
    cv::threshold(m_, m_, value_0_1*max_quant(), max_quant(), cv::THRESH_BINARY);
    return *this;
  }

  /**
   * Sets all pixels with a value smaller than the corresponding component of the colour
   * to black.
   * @param double value_0_1
   * @return basic_image&
   */
  inline basic_image& threshold(const color_t& c)
  {
    if(is_gray() || (c[0]==c[1] && c[0]==c[2])) {
      int cc = c[0];
      if(cc < c[1]) cc = c[1];
      if(cc < c[2]) cc = c[2];
      if(cc >= max_quant()) cc = max_quant()-1;
      cv::threshold(m_, m_, cc, max_quant(), cv::THRESH_BINARY);
    } else {
      mat_vector_t chs = split_channels();
      if(chs.size() > 3) chs.resize(3);
      for(typename mat_vector_t::size_type i=0; i<chs.size(); ++i) {
        int cc = c[i] >= max_quant() ? max_quant()-1 : c[i];
        cv::threshold(chs[i], chs[i], cc, max_quant(), cv::THRESH_BINARY);
      }
      merge_channels(chs);
    }
    return *this;
  }

  /**
   * Adaptive threshold
   * @param crd_t radius=0
   * @param bool gauss=false
   * @return basic_image&
   */
  inline basic_image& threshold_adaptive(crd_t radius=0, bool gauss=false)
  {
    to_grayscale();
    if(radius<2) radius = 2;
    cv::adaptiveThreshold(m_, m_, max_quant(), gauss ? cv::ADAPTIVE_THRESH_GAUSSIAN_C :
      cv::ADAPTIVE_THRESH_MEAN_C, cv::THRESH_BINARY, 2*radius+1, 0);
    return *this;
  }

  // </editor-fold>

  // <editor-fold desc="dilate, erode" defaultstate="collapsed">
  /**
   * Dilates the image
   * @param unsigned radius
   * @return basic_image&
   */
  inline basic_image& dilate(unsigned radius=1)
  {
    cv::dilate(m_, m_, cv::getStructuringElement(cv::MORPH_ELLIPSE,
      cv::Size(2*radius+1, 2*radius+1), cv::Point(radius, radius)
    ));
    return *this;
  }

  /**
   * Erodes the image
   * @param unsigned radius
   * @return basic_image&
   */
  inline basic_image& erode(unsigned radius=1)
  {
    cv::erode(m_, m_, cv::getStructuringElement(cv::MORPH_ELLIPSE,
      cv::Size(2*radius+1, 2*radius+1), cv::Point(radius, radius)
    ));
    return *this;
  }

  // </editor-fold>

  // <editor-fold desc="blur" defaultstate="collapsed">

  /**
   * blur() methods:
   */
  typedef enum { blur_box=0, blur_gauss, blur_median } blur_type;

  /**
   * Blur with given "radius" (kernel size 2*dist+1).
   * @param int kernel_radius=1
   * @return basic_image&
   */
  inline basic_image& blur(int kernel_radius=1, blur_type how=blur_gauss)
  {
    if((kernel_radius=kernel_radius > 32 ? 32 : kernel_radius) > 0) {
      if(how == blur_gauss) {
        cv::GaussianBlur(m_, m_, cv::Size(2*kernel_radius+1, 2*kernel_radius+1), 0, 0 );
      } else if(how == blur_median) {
        cv::medianBlur(m_, m_, 2*kernel_radius+1);
      } else {
        cv::blur(m_, m_, cv::Size(2*kernel_radius+1, 2*kernel_radius+1));
      }
    }
    return *this;
  }

  // </editor-fold>

  // <editor-fold desc="sobel, scharr, laplacian, gradient" defaultstate="collapsed">

  /**
   * Applies Sobel spatial image derivative. Returns the absolute values of
   * each (pixel) location derivative (values are 0 to something >0).
   * The "kernel_radius" expands to kernel size 2*kernel_radius+1.
   * Note: dx, dy are the filter orders in x/y, valid values
   * 0(=component off) to kernel_size. On scale_to_max_quantum scales==true
   * the values are shifted half the quantum range (127 for 8 bit), so that
   * negative gradients are <128, positive gradient values >=128;
   * @param int dx=1
   * @param int dy=1
   * @param int kernel_radius=0
   * @return basic_image
   */
  inline basic_image& sobel(int dx=1, int dy=1, int kernel_radius=0)
  {
    kernel_radius = kernel_radius<0 ? 0:kernel_radius;
    kernel_radius = 2*kernel_radius+1;
    dx=dx<0?0: (dx>kernel_radius? kernel_radius:dx);
    dy=dy<0?0: (dy>kernel_radius? kernel_radius:dy);
    mat_t m;
    cv::Sobel(m_, m, CV_32F, dx, dy, kernel_radius, 1, 0, cv::BORDER_DEFAULT);
    m = abs(m);
    double mi, ma;
    cv::minMaxLoc(m, &mi, &ma);
    m *= (double)max_quant()/ma;
    m_ = m;
    return *this;
  }

  /**
   * Applies Scharr filter. Returns the absolute values of each (pixel) location
   * derivative (values are normed 0 to max_quant()). The data type of each pixel
   * is float32.
   * Note: dx, dy are the filter orders in x/y, valid values
   * 0(=component off) to kernel_size. On scale_to_max_quantum scales==true
   * the values are shifted half the quantum range (127 for 8 bit), so that
   * negative gradients are <128, positive gradient values >=128;
   * @param int dx
   * @param int dy
   * @return basic_image&
   */
  inline basic_image& scharr(int dx, int dy)
  {
    dx=dx<0?0: (dx>1? 1:dx);
    dy=dy<0?0: (dy>1? 1:dy);
    mat_t m;
    cv::Scharr(m_, m, CV_32F, dx, dy, 1, 0, cv::BORDER_DEFAULT);
    m = abs(m);
    double mi, ma;
    cv::minMaxLoc(m, &mi, &ma);
    m *= (double)max_quant()/ma;
    m_ = m;
    return *this;
  }

  /**
   * Laplacian ("second derivative"). The "kernel_radius" expands to kernel size
   * 2*kernel_radius+1.
   * @param int kernel_radius=0
   * @return basic_image&
   */
  inline basic_image& laplacian(int kernel_radius=0)
  {
    kernel_radius=kernel_radius<0 ? 0:kernel_radius;
    mat_t m0; mat().convertTo(m0, CV_32F);
    mat_t m(m0.rows, m0.cols, m0.type());
    cv::Laplacian(m0, m, CV_32F, 2*kernel_radius+1, 1, 0, cv::BORDER_DEFAULT);
    m = abs(m);
    double mi, ma;
    cv::minMaxLoc(m, &mi, &ma);
    m *= (double)max_quant()/ma;
    m_ = m;
    return *this;
  }

  /**
   * Gradient calculation from Sobel in x and y. Returns ("geometric") magnitude
   * of x/y values (normed to max_quant()). The data type of the image object
   * changes to float_t. The "kernel_radius" expands to kernel size
   * 2*kernel_radius+1.
   *
   * @param unsigned kernel_radius=0
   * @return basic_image&
   */
  inline basic_image& gradient(int kernel_radius=0)
  {
    if(channels()>1) to_grayscale();
    cv::Mat1d xm, ym;
    {
      cv::Sobel(m_, xm, CV_64F, 1,0, 2*kernel_radius+1, 1,0, cv::BORDER_DEFAULT);
      cv::Sobel(m_, ym, CV_64F, 0,1, 2*kernel_radius+1, 1,0, cv::BORDER_DEFAULT);
    }
    if(ym.size != xm.size) {
      dbg("Sobel x/y derivative matrices SIZE MISSMATCH");
    } else {
      cv::Mat1d mags(xm.rows, xm.cols);
      cv::magnitude(xm,ym,mags);
      double mi, ma;
      mags = abs(mags);
      cv::minMaxLoc(mags, &mi, &ma);
      if(ma!=0) mags *= (double)max_quant()/ma;
      mags.convertTo(m_, sizeof(float_t)==8 ? CV_64F : CV_32F);
      //dbg("gradient(): max magnitude=" << ma);
    }
    return *this;
  }
  // </editor-fold>

  // <editor-fold desc="edge/corner detection" defaultstate="collapsed">

  /**
   * Canny edge filter. The `kernel_radius` expands to 2*kernel_radius+1.
   * @param double thres1=100
   * @param double thres2=120
   * @param int kernel_radius=1
   * @return basic_image&
   */
  inline basic_image& edges_canny(crd_t thres1=100, crd_t thres2=120, int kernel_radius=1)
  { cv::Canny(m_, m_, thres1, thres2, 2*kernel_radius+1, true); return *this;}

  /**
   * Harris corner detector. The `kernel_radius` expands to 2*kernel_radius+1.
   * @param k=0.01
   * @param block_size=3
   * @param kernel_radius=1
   * @return basic_image&
   */
  inline basic_image& corners_harris(double k=0.01, int block_size=3, int kernel_radius=1)
  {
    if(channels() > 1) to_grayscale();
    mat_t m = mat_t::zeros(m_.size(), CV_32FC1);
    cv::cornerHarris(m_, m, block_size, 2*kernel_radius+1, k, cv::BORDER_DEFAULT );
    cv::normalize(m, m, 0, max_quant(), cv::NORM_MINMAX, CV_32FC1, mat_t());
    cv::convertScaleAbs(m, m_);
    return *this;
  }
  // </editor-fold>

public:  // Methods returning non-image data types

  // <editor-fold desc="contours" defaultstate="collapsed">

  typedef enum {
    cfa_code=0, cfa_none=1, cfa_simple=2, cfa_l1=3, cfa_kcos=4, cfa_link_runs=5
  } contour_chain_approx_t;

  /**
   * Find contours. Every nonzero (nonblack) pixel is used. Colour images are
   * (copy) converted to gray scale.
   * @param bool external_only=false
   * @param contour_chain_approx_t chain_approx=cf_simple
   * @return contour_list_t
   */
  inline contour_list_t contours(bool external_only=false,
          contour_chain_approx_t chain_approx=cfa_simple) const
  {
    contour_list_t contours;
    std::vector<std::vector<cv::Point2i> > cvcontours;
    {
      std::vector<cv::Vec4i> hierarchy;
      basic_image im(m_);
      if(im.channels()>1) im.to_grayscale();
      cv::findContours(im.mat(), cvcontours, hierarchy, external_only ?
        CV_RETR_EXTERNAL : CV_RETR_LIST, chain_approx);
    }
    while(!cvcontours.empty()) {
      contours.push_back((contour_t)cvcontours.back());
      cvcontours.pop_back();
    }
    return contours;
  }

  /**
   * Find contours. Every nonzero (nonblack) pixel is used. Colour images are
   * (copy) converted to gray scale. Contours are organised in external boundaries
   * and holes. "Holes in holes" are treated as external structures again, and
   * so on (even contour levels are external, odd levels are holes).
   * @param contour_chain_approx_t chain_approx=cf_simple
   * @return contour_ccomp_list_t
   */
  inline contour_ccomp_list_t contours_ccomp(
          contour_chain_approx_t chain_approx=cfa_simple) const
  {
    contour_ccomp_list_t contours;
    std::vector<std::vector<cv::Point2i> > cvcontours;
    std::vector<cv::Vec4i> cvhierarchy;
    {
      basic_image im(m_);
      if(im.channels()>1) im.to_grayscale();
      cv::findContours(im.mat(), cvcontours, cvhierarchy, CV_RETR_CCOMP, chain_approx);
    }
    while(!cvcontours.empty()) {
      if(cvhierarchy.back()[3] < 0) {
        contours.outer.push_back((contour_t)cvcontours.back());
      } else {
        contours.inner.push_back((contour_t)cvcontours.back());
      }
      cvcontours.pop_back();
      cvhierarchy.pop_back();
    }
    return contours;
  }

  /**
   * Finds contours and returns them as tree structure, where inner contours
   * the child elements. The root node contour is always the whole image (t.m.
   * a rectangular contour with 4 points: (0,0,image width, image height).
   * @param contour_chain_approx_t chain_approx=cfa_simple
   * @return contour_nested_list_t
   */
  contour_nested_list_t contours_tree(
    contour_chain_approx_t chain_approx=cfa_simple) const
  {
    if(!width() || !height()) return (contour_t(rect_t(point_t(0,0),point_t(0,0))));
    contour_nested_list_t contours(contour_t(rect_t(0,0,width()-1,height()-1)));
    {
      std::vector<contour_nested_list_t> nodes; // Nodes (contents);
      std::vector<int> parent_indices, parents;
      {
        std::vector<std::vector<cv::Point2i> > cvcontours;
        std::vector<cv::Vec4i> cvhierarchy;
        {
          basic_image im(m_);
          if(im.channels()>1) im.to_grayscale();
          cv::findContours(im.mat(), cvcontours, cvhierarchy, CV_RETR_TREE, chain_approx);
        }
        std::map<int,int> avail_parents;
        for(unsigned i=0; i<cvcontours.size(); ++i) {
          unsigned parent = cvhierarchy[i][3];
          nodes.push_back(contour_nested_list_t(cvcontours[i]));
          parent_indices.push_back(parent);
          avail_parents[parent] = 0;
        }
        for(typename std::map<int,int>::iterator it=avail_parents.begin();
                it!=avail_parents.end(); ++it) {
          parents.push_back(it->first);
        }
        std::sort(parents.begin(), parents.end());
        if(parents.front()==-1) parents.erase(parents.begin());
      } while(0);

      // --> Rollup parents from back to front
      while(!parents.empty()) {
        int pi = parents.back();
        parents.pop_back();
        for(unsigned i=0; i<parent_indices.size(); ++i) {
          if(parent_indices[i]==pi) {
            nodes[pi].children_.push_back(nodes[i]); // push copy
            nodes[pi].children_.back().parent_ = &(nodes[pi]);
            nodes[i].clear(); // save mem
          }
        }
      }
      for(unsigned i=0; i<parent_indices.size(); ++i) {
        if(parent_indices[i] < 0) {
          contours.children_.push_back(nodes[i]);
          contours.children_.back().parent_ = &contours;
          nodes[i].clear();
        }
      }
    }
    return contours;
  }

  // </editor-fold>

  // <editor-fold desc="hough lines/circles" defaultstate="collapsed">
  /**
   * Hough line transform
   * @param double acc_vote_threshold=1
   * @param double min_line_length=10
   * @param double max_line_point_gap=2
   * @param double acc_angle_resolution_deg=1.0
   * @param double acc_distance_resolution_px=1
   * @return line_list_t (vector of line_t)
   */
  inline line_list_t hough_lines(
    double acc_vote_threshold=1,
    double min_line_length=10, double max_line_point_gap=2,
    double acc_angle_resolution_deg=1.0, double acc_distance_resolution_px=1
  ){
    std::vector<cv::Vec4i> hlines;
    line_list_t lines;
    cv::HoughLinesP(m_, hlines,
      acc_distance_resolution_px, acc_angle_resolution_deg*M_PI/360,
      acc_vote_threshold, min_line_length, max_line_point_gap
    );
    for(size_t i=0; i<hlines.size(); ++i) {
      const cv::Vec4i &l = hlines[i];
      lines.push_back(line_t(point_t(l[0],l[1]), point_t(l[2],l[3])));
    }
    return lines;
  }

  /**
   * Hough circle transform
   * Note: min_centre_distance==-1 --> max_radius
   * @param double max_radius=10
   * @param double min_radius=5
   * @param double canny_threshold=50
   * @param acc_vote_threshold=50
   * @param min_centre_distance=-1
   * @param relative_acc_resolution=1
   * @return circle_list_t (vector of circle_t)
   */
  inline circle_list_t hough_circles(
    double max_radius=10,
    double min_radius=5,
    double canny_threshold=50,
    double acc_vote_threshold=50,
    double min_centre_distance=-1,
    double relative_acc_resolution=1
  ){
    circle_list_t circles;
    std::vector<cv::Vec3f> hc;
    int res = relative_acc_resolution < 2 ? 1 : relative_acc_resolution;
    if(acc_vote_threshold < 1) acc_vote_threshold = 1;
    if(canny_threshold < 1) canny_threshold = 1;
    double dist = min_centre_distance <= 0 ? max_radius : min_centre_distance;
    cv::HoughCircles(
      m_, hc, CV_HOUGH_GRADIENT, res, dist,
      canny_threshold, acc_vote_threshold, min_radius, max_radius
    );
    for(size_t i=0; i<hc.size(); ++i) {
      const cv::Vec3f &c = hc[i];
      circles.push_back(circle_t(c[0],c[1],c[2]));
    }
    return circles;
  }
  // </editor-fold>

  // <editor-fold desc="histogram" defaultstate="collapsed">
  /**
   * @see basic_histogram<>
   * @param int num_bins=0
   * @param const basic_image& mask=basic_image()
   * @return histogram_t
   */
  histogram_t histogram(int num_bins=0, const basic_image& mask=basic_image()) const
  { return histogram_t(*this, num_bins, mask); }

  // </editor-fold>

public: // Drawing

  // <editor-fold desc="draw, fill, floodfill, mask" defaultstate="collapsed">

  /**
   * Draw a point with defined size.
   * @param const point_t &point
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const point_t &point, const color_t &c=color_t(), double size=1)
  {
    if(size < 2) {
      color_at(point.x(), point.y(), c);
    } else {
      cv::circle(m_, (cv::Point)point, 1, (cv::Scalar) c, size);
    }
    return *this;
  }

  /**
   * Draw a line with defined width.
   * @param const line_t &line
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const line_t &line, const color_t &c=color_t(), double size=1)
  { cv::line(m_, (cv::Point)line.p0(), (cv::Point)line.p1(), (cv::Scalar) c, size);
      return *this; }

  /**
   * Draw a circle outline with defined line width.
   * @param const circle_t &circle
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const circle_t &circle, const color_t &c=color_t(), double size=1)
  {
    if(circle.r() <= 0 || circle.x() <=0 || circle.y() <=0) return *this;
    cv::circle(m_, (cv::Point)circle, circle.r(), (cv::Scalar) c, size); return *this;
  }

  /**
   * Draw an ellipse outline with defined line width.
   * @param const ellipse_t &e
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const ellipse_t &e, const color_t &c=color_t(), double size=1)
  { cv::ellipse(m_, (cv::RotatedRect)e, (cv::Scalar) c, size); return *this; }

  /**
   * Draw an upright rectangle outline
   * @param const rect_t &r
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const rect_t &r, const color_t &c=color_t(), double size=1)
  { cv::rectangle(m_, (cv::Rect)r, (cv::Scalar) c, size); return *this; }

  /**
   * Draw a rotated rectangle outline
   * @param const rotated_rect_t &r
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const rotated_rect_t &r, const color_t& c=color_t(), double size=1)
  {
    if(r.size().x() && r.size().y() && size >= 1) {
      cv::polylines(m_, std::vector<std::vector<cv::Point2i> >(1, r), true,
          (cv::Scalar) c, (int)size);
    }
    return *this;
  }

  /**
   * Draw a polynom outline.
   * @param const polynom_t &p
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const polynom_t &p, const color_t &c=color_t(), double size=1)
  {
    if(!p.empty()) {
      cv::polylines(m_, std::vector<std::vector<cv::Point2i> >(1,p), p.closed(),
          (cv::Scalar) c, (int)size);
    }
    return *this;
  }

  /**
   * Draw a contour outline
   * @param const contour_t &p
   * @param const color_t &c=color_t()
   * @param double size=1
   * @return basic_image& (*this)
   */
  inline basic_image& draw(const contour_t &p, color_t c=color_t(), double size=1)
  {
    if(!p.empty()) {
      cv::polylines(m_, std::vector<std::vector<cv::Point2i> >(1,p), true,
      (cv::Scalar) c, (int)size);
    }
    return *this;
  }

  /**
   * Draws a channel of the histogram
   * @param const histogram_t &hist
   * @param int channel
   * @param const img_t::rect_t &where=img_t::rect_t()
   * @param color_t curve_color = white
   * @param color_t background_color = black
   * @param color_t border_color = gray
   * @param int grid_div_x = 4
   * @param int grid_div_y = 4
   * @param color_t grid_color = gray
   * @param int line_width = 0
   * @return basic_image&
   */
  inline basic_image& draw(const histogram_t &hist, int channel,
      const rect_t &where=rect_t(), color_t curve_color=color_t(255),
      color_t background_color=color_t(), color_t border_color=color_t(80),
      int grid_div_x=4, int grid_div_y=4, color_t grid_color=color_t(80), int line_width=0
      )
  {
    hist.draw_channel(*this, channel, where, curve_color, background_color, border_color,
      grid_div_x, grid_div_y, grid_color, line_width);
    return *this;
  }

  /**
   * Draws an image into another one the image is scaled to the rect size.
   * @param const basic_image& im
   * @param rect_t where
   * @return basic_image&
   */
  inline basic_image& draw(const basic_image& im, rect_t where)
  {
    if(where.width()<1 || where.height()<1 || im.empty()) return *this;
    basic_image iim(im);
    iim.resize_stretch(where.width(), where.height());
    if(where.br().x()>=width() || where.br().y()>=height()) {
      where = rect_t(point_t(where.p()), point_t(width()-1, height()-1));
      iim.crop(where.width(), where.height());
      dbg2("draw(): Modified (rect_t) where=" << where);
    }
    cv::Mat mm = m_((cv::Rect) where);
    iim.mat().copyTo(mm);
    return *this;
  }

  /**
   * Draws an image into another one with mask.
   * @param const basic_image& im
   * @param const basic_image& mask
   * @return basic_image&
   */
  inline basic_image& draw(const basic_image& im, const basic_image& mask)
  {
    if(color_space() != im.color_space()) {
      basic_image im2(im);
      im2.color_space(color_space());
      im2.mat().copyTo(m_, mask.mat());
    } else {
      im.mat().copyTo(m_, mask.mat());
    }
    return *this;
  }

  /**
   * Draws the histogram
   * @param const histogram_t &hist
   * @param const img_t::rect_t &where=img_t::rect_t()
   * @return basic_image&
   */
  inline basic_image& draw(const histogram_t &hist, const rect_t &where)
  { hist.draw(*this, where); return *this; }

  /**
   * Fill the whole image with a defined colour.
   * @return basic_image& (*this)
   */
  inline basic_image& fill(const color_t &c)
  { m_ = (cv::Scalar)c; return *this; }

  /**
   * Fill am area defined by a CV point vector
   * @param const std::vector<cv::Point2i> &p
   * @param const color_t &c=color_t()
   * @return basic_image& (*this)
   */
  inline basic_image& fill(const std::vector<cv::Point2i> &p, const color_t &c=color_t())
  { if(p.empty()) return *this;
    std::vector<std::vector<cv::Point2i> > v; v.push_back((std::vector<cv::Point2i>)p);
    cv::fillPoly(m_, v,(cv::Scalar) c); return *this; }

  /**
   * Fill a circle area
   * @param const circle_t &circle
   * @param const color_t &c=color_t()
   * @return basic_image& (*this)
   */
  inline basic_image& fill(const circle_t &cr, const color_t &c=color_t())
  { if(cr.r()>0) cv::circle(m_, (cv::Point)cr, cr.r(), (cv::Scalar) c, ls_filled);
      return *this; }

  /**
   * Fill an ellipse area
   * @param const ellipse_t &e
   * @param const color_t &c=color_t()
   * @return basic_image& (*this)
   */
  inline basic_image& fill(const ellipse_t &e, const color_t &c=color_t())
  { cv::ellipse(m_, (cv::RotatedRect)e, (cv::Scalar) c, ls_filled); return *this; }

  /**
   * Fill an upright rectangle area
   * @param const rect_t &r
   * @param const color_t &c=color_t()
   * @return basic_image& (*this)
   */
  inline basic_image& fill(const rect_t &r, const color_t &c=color_t())
  { cv::rectangle(m_, (cv::Rect)r, (cv::Scalar) c, ls_filled); return *this; }

  /**
   * Fill the area of a rotated rectangle
   * @param const rotated_rect_t &r
   * @param const color_t &c=color_t()
   * @return basic_image& (*this)
   */
  inline basic_image& fill(const rotated_rect_t &r, const color_t &c=color_t())
  { return fill((std::vector<cv::Point2i>)r, c); }

  /**
   * Fill the area of a polynom
   * @param const polynom_t &p
   * @param const color_t &c=color_t()
   * @return basic_image& (*this)
   */
  inline basic_image& fill(const polynom_t &p, const color_t &c=color_t())
  { return fill((std::vector<cv::Point2i>)p, c); }

  /**
   * Flood fill from seed point with defined colour and tolerance. The tolerance
   * of type color_t can be specified as integer as well and will be implicitly
   * converted to the corresponding gray scale colour.
   * @param const point_t& seed_point
   * @param const color_t& c
   * @param const color_t& tolerance
   * @param bool difference_from_seed_point=false
   * @return basic_image&
   */
  inline basic_image& floodfill(const point_t &seed_point, const color_t &c,
      const color_t &tolerance, bool difference_from_seed_point=false)
  {
    int flags = 4|(255<<8)|(difference_from_seed_point ? cv::FLOODFILL_FIXED_RANGE:0);
    cv::Mat mask;
    cv::floodFill(m_, mask, (cv::Point) seed_point, (cv::Scalar)c , 0,
        (cv::Scalar)tolerance, (cv::Scalar)tolerance, flags);
    return *this;
  }

  /**
   * Flood fill from seed point with defined colour and tolerance. The tolerance
   * of type color_t can be specified as integer as well and will be implicitly
   * converted to the corresponding gray scale colour.
   * @param const point_t& seed_point
   * @param const color_t& c
   * @param const color_t& tolerance
   * @param bool difference_from_seed_point
   * @param basic_image& mask
   * @return basic_image&
   */
  inline basic_image& floodfill(const point_t &seed_point, const color_t &c,
      const color_t &tolerance, bool difference_from_seed_point, basic_image& mask)
  {
    int flags = 4|(255<<8)|(difference_from_seed_point ? cv::FLOODFILL_FIXED_RANGE:0);
    mask.to_grayscale().extend(width()+2, height()+2, true); // ensure type, size
    cv::floodFill(m_, mask.mat(), (cv::Point) seed_point, (cv::Scalar)c , 0,
        (cv::Scalar)tolerance, (cv::Scalar)tolerance, flags);
    mask.crop(1, 1, width(), height());
    return *this;
  }

  /**
   * Flood fill from seed point with the colour of the seed point. The tolerance
   * of type color_t can be specified as integer as well and will be implicitly
   * converted to the corresponding gray scale colour.
   * @param const point_t &seed_point
   * @param const color_t &tolerance=0
   * @param bool difference_from_seed_point=false
   * @return basic_image&
   */
  inline basic_image& floodfill(const point_t &seed_point, const color_t &tolerance,
      bool difference_from_seed_point=false)
  {
    color_t c = color_at(seed_point);
    int flags = 4|(255<<8)|(difference_from_seed_point ? cv::FLOODFILL_FIXED_RANGE:0);
    mat_t mask;
    cv::floodFill(m_, mask, (cv::Point) seed_point, (cv::Scalar)c , 0,
        (cv::Scalar)tolerance, (cv::Scalar)tolerance, flags);
    return *this;
  }

  /**
   * Returns the pixels that would be modified when flood filling from seed
   * point with defined colour and tolerance. The tolerance
   * of type color_t can be specified as integer as well and will be implicitly
   * converted to the corresponding gray scale colour.
   * @param const point_t& seed_point
   * @param const color_t& c
   * @param const color_t& tolerance=0
   * @param bool difference_from_seed_point=false
   * @return basic_image&
   */
  inline basic_image floodfill_mask(const point_t &seed_point,
      const color_t &tolerance=0, bool difference_from_seed_point=false)
  {
    int flags = 4|(255<<8)|cv::FLOODFILL_MASK_ONLY;
    flags |= difference_from_seed_point ? cv::FLOODFILL_FIXED_RANGE : 0;
    mat_t mask = mat_t::zeros(m_.rows+2, m_.cols+2, CV_8U);
    cv::floodFill(m_, mask, (cv::Point) seed_point, 255, 0, (cv::Scalar)tolerance,
        (cv::Scalar)tolerance, flags);
    return mask(cv::Rect(1,1,m_.cols-2,m_.rows-2)).clone();
  }

  /**
   * Masks the image.
   * @param const basic_image & mask
   * @param const color_t& repl_color=color_t::black
   * @return basic_image&
   */
  inline basic_image& mask(const basic_image & mask, const color_t& repl_color=color_t::black)
  {
    basic_image m = mask;
    m.negate();
    mat().setTo((cv::Scalar)repl_color, m.mat());
    return *this;
  }

  // </editor-fold>

protected:

  // <editor-fold desc="instance variables" defaultstate="collapsed">
  mat_t m_;
  std::basic_string<char_t> file_;
  bool is_hsv_;
  // </editor-fold>

};

// <editor-fold desc="basic_image operator <<" defaultstate="collapsed">
/**
 * Output stream dump
 * @param std::ostream&
 * @param basic_image<crd_t,val_t,float_t>
 * @return std::ostream&
 */
template <typename crd_t, typename val_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_image<crd_t,val_t,float_t> im)
{
  os << "{" << std::endl
     << "  width: " << im.width() << "," << std::endl
     << "  height: " << im.height() << "," << std::endl
     << "  type: \"" << (im.is_gray()?"gray":im.is_hsv()?"hsv":"rgb")<< "\""
     << "," << std::endl
     << "  data_type: \"";
  switch(im.mat().depth()) {
    case CV_8U: os << "uint8"; break;
    case CV_8S : os << "int8"; break;
    case CV_16U: os << "uint16"; break;
    case CV_16S: os << "int16"; break;
    case CV_32S: os << "int32"; break;
    case CV_32F: os << "float32"; break;
    case CV_64F: os << "float64"; break;
    default: os << "unknown";
  }
  os << "\"," << std::endl
     << "  channels: " << im.channels() << "," << std::endl
     << "  min: \"" << im.min() << "\"," << std::endl
     << "  max: \"" << im.max() << "\"," << std::endl
     << "  mean: \"" << im.mean() << "\"," << std::endl
     << "  mat_continuous: " << im.mat().isContinuous() << "," << std::endl
     << "  mat_is_submatrix: " << im.mat().isSubmatrix() << "," << std::endl
     << "  mat_dimensions: " << im.mat().dims << "," << std::endl
     << "  mat_refcount: " << (im.mat().refcount ? *im.mat().refcount : 0) << "," << std::endl
     << "  mat_element_size: " << im.mat().elemSize() << "," << std::endl
     << "  mat_element_cell_size: " << im.mat().elemSize1() << "," << std::endl
     << "  mat_buffer_size: " << im.mat().total() << std::endl
      ;

  os << "}" << std::endl;
  return os;
}
// </editor-fold>

// <editor-fold desc="basic_histogram" defaultstate="collapsed">
template <typename crd_t, typename val_t, typename float_t>
class basic_histogram
{
public:

  typedef basic_image<crd_t,val_t,float_t> img_t;
  typedef typename img_t::rect_t rect_t;
  typedef typename img_t::line_t line_t;
  typedef typename img_t::point_t point_t;
  typedef typename img_t::color_t color_t;
  typedef std::vector<cv::Mat> matv_t;

public:

  inline basic_histogram() : histograms_(), is_hsv_(false)
  { ; }

  inline basic_histogram(const basic_histogram& hist) : histograms_(hist.histograms_),
    is_hsv_(hist.is_hsv_)
  { ; }

  /**
   * Directly calculates the histogram. On error, the histogram is empty().
   * @param const img_t &im
   * @param const img_t& mask=img_t()
   * @param unsigned num_bins=0
   */
  inline basic_histogram(const img_t &im, unsigned num_bins=0, const img_t& mask=img_t())
    : histograms_()
  { try{ calculate(im, num_bins, mask); } catch(...) { histograms_.clear(); } }

  virtual ~basic_histogram()
  { ; }

  inline basic_histogram& operator=(const basic_histogram& h)
  { histograms_=h.histograms_; is_hsv_=h.is_hsv_; return *this; }

public:

  /**
   * Returns if the histograms is empty.
   * @return bool
   */
  inline bool empty() const throw()
  { return histograms_.empty(); }

  /**
   * Returns the number of histogram bins (same for all channels)
   * @return unsigned
   */
  inline unsigned num_bins() const throw()
  { return empty() ? 0 : histograms_.front().rows; }

  /**
   * Returns the maximum bin (accumulator) count of a given channel.
   * @param unsigned channel
   * @return int
   */
  inline int max(unsigned channel) const
  {
    if(channel >= histograms_.size()) return 0;
    double mi, ma; cv::minMaxLoc(histograms_[channel], &mi, &ma);
    return ma;
  }

  /**
   * Returns the minimum bin (accumulator) count of a given channel.
   * @param unsigned channel
   * @return int
   */
  inline int min(unsigned channel) const
  {
    if(channel >= histograms_.size()) return 0;
    double mi, ma; cv::minMaxLoc(histograms_[channel], &mi, &ma);
    return mi;
  }

  /**
   * Returns the raw std::vector<cv::Mat> as const reference.
   * @return matv_t
   */
  inline const matv_t & mats() const throw()
  { return histograms_; }


  /**
   * Returns the number of histogram bins (same for all channels). This method
   * is used to operate as end index for gray scale histograms combined with
   * operator [].
   * @return unsigned
   */
  inline unsigned size() const throw()
  { return num_bins(); }


  /**
   * Returns the bin counter of the FIRST CHANNEL at `index` or 0.
   * @param unsigned index
   * @return unsigned
   */
  inline unsigned operator[](unsigned index) const
  { return histograms_.empty() || index >= num_bins() ? 0 :
      (unsigned)histograms_[0].template at<float>(index); }


public:

  /**
   * Calculates the histogram for up to 3 channels. If the specified number of
   * bins is 0, the colour range resolution is assumed (im.max_quant()).
   * @param const img_t &im
   * @param unsigned num_bins=0
   * @param const img_t& mask=img_t()
   * @return basic_histogram&
   */
  basic_histogram& calculate(const img_t &im, unsigned num_bins=0, const img_t& mask=img_t())
  {
    is_hsv_ = im.is_hsv();
    histograms_.clear();
    if(im.empty()) return *this;
    if(num_bins==0) num_bins = im.max_quant()/2+1; // 0xff+1 for 8 bit
    int channels[] = { 0, 1, 2, 3 };
    unsigned n_channels = sizeof(channels) / sizeof(int);
    int n_channel_bins[]  = { (int)num_bins, (int)num_bins, (int)num_bins, (int)num_bins };
    float rng_c0[]={0,is_hsv_?180:256}, rng_c1[]={0,256}, rng_c2[]={0,256}, rng_c3[]={0,256};
    const float* ranges[] = { rng_c0, rng_c1, rng_c2, rng_c3 };
    cv::Mat msk;
    if(!mask.empty()) {
      img_t imsk = mask;
      msk = imsk.to_grayscale().threshold(0.5).mat().clone();
    }
    for(unsigned i=0; i<im.channels() && i<n_channels; ++i) {
      cv::Mat hst;
      cv::calcHist(&im.mat(), 1, channels+i, msk, hst, 1, n_channel_bins+i,
          ranges+i, true, false);
      histograms_.push_back(hst.clone()); // clone: ensure hist not re-referenced
    }
    return *this;
  }

  /**
   * Normalises the histogram (all channels if channel<0, selected channel if
   * >=0) to a range from 0 to max_value. This would be the height of the
   * drawing area if the histogram is plotted.
   * @param int max_value
   * @return int channel=-1
   */
  basic_histogram& normalize(int max_value, int channel=-1)
  {
    if(max_value<=1 || channel>=histograms_.size()) return *this;
    int from = 0, to = histograms_.size()-1;
    if(channel >= 0) from = to = channel;
    for(unsigned i=from; i<=to; ++i) {
      cv::normalize(histograms_[i], histograms_[i], 0, max_value, cv::NORM_MINMAX,
        -1, cv::Mat());
    }
    return *this;
  }

  /**
   * Draws a channel of the histogram
   * @param img_t& im
   * @param int channel
   * @param const img_t::rect_t &where=img_t::rect_t()
   * @param color_t curve_color = white
   * @param color_t background_color = black
   * @param color_t border_color = gray
   * @param int grid_div_x = 4
   * @param int grid_div_y = 4
   * @param color_t grid_color = gray
   * @param int line_width = 0
   * @return basic_histogram&
   */
  const basic_histogram& draw_channel(img_t& im, int channel, const rect_t &where=rect_t(),
      color_t curve_color=color_t(255), color_t background_color=color_t(),
      color_t border_color=color_t(80), int grid_div_x=4, int grid_div_y=4,
      color_t grid_color=color_t(80), int line_width=0
      ) const
  {
    if(channel<0 || channel>=(int)histograms_.size() || where.width()<5 || where.height()<5) {
      return *this;
    }
    if(line_width<1) {
      line_width = (im.height()>im.width()?im.width():im.height())/500;
      line_width = line_width<0 ? 0 : line_width > 5 ? 5 : line_width;
    }

    cv::Mat hst = histograms_[channel].clone();
    int w = where.width()-2; // 2*1px border
    int h = where.height()-2;
    double x_scale = (double)w / num_bins();
    cv::normalize(hst, hst, 0, h, cv::NORM_MINMAX, -1, cv::Mat());
    im.fill(where, background_color);

    if(grid_div_x > 1) {
      for(int j=::round((double)w/grid_div_x); j<=w; j+=::round((double)w/grid_div_x)) {
        im.draw(line_t(
          point_t(where.x()+j, where.y()+1),
          point_t(where.x()+j, where.y()+1+h)
        ), grid_color, line_width);
      }
    }

    if(grid_div_y > 1) {
      for(int j=::round((double)h/grid_div_y); j<=h; j+=::round((double)h/grid_div_y)) {
        im.draw(line_t(
          point_t(where.x()+1, where.y()+j),
          point_t(where.x()+1+w, where.y()+j)
        ), grid_color, line_width);
      }
    }

    im.draw(where, border_color, line_width);

    for(int j=0; j<(int)num_bins(); j++) {
      im.fill(rect_t(
          point_t(where.x()+1+(x_scale*(j+0)), where.y()+1+h-hst.at<float>(j)),
          point_t(where.x()+1+(x_scale*(j+1)), where.y()+1+h)
        ), curve_color);
    }
    return *this;
  }

  /**
   * Draw a histogram with three graphs (channels left to right)
   * @param img_t& im
   * @param const rect_t &where=rect_t()
   * @return basic_histogram&
   */
  const basic_histogram& draw(img_t& im, const rect_t& where) const
  {
    color_t background_color = color_t();
    color_t border_color = color_t(80);
    int grid_div_x=4;
    int grid_div_y=4;
    color_t grid_color = color_t(80);
    int line_width;
    line_width = (im.height()>im.width()?im.width():im.height())/500;
    line_width = line_width<0 ? 0 : line_width > 5 ? 5 : line_width;
    color_t curve_color = color_t(255);
    int n_ch = histograms_.size();
    if(where.width() < (5+(2*line_width))*n_ch || where.height() < (5+(2*line_width))*1) {
      return *this;
    }
    std::vector<color_t> curve_colors;
    curve_colors.push_back("#ff0000"); // R / H / gray
    curve_colors.push_back("#00ff00"); // G / S
    curve_colors.push_back("#0000ff"); // B / V
    curve_colors.push_back("#ffffff"); // unused (alpha)
    float lw = line_width;
    float gap = 5*lw;
    float rx = where.x() + lw;
    float ry = where.y() + lw;
    float hw = ((float)(where.width()-2*lw)/n_ch)-(n_ch>1?(gap)/(n_ch-1):0)-lw;
    float hh = (float)where.height()-2*lw;
    for(int i=0; i<n_ch; ++i) {
      this->draw_channel(im, i, rect_t(point_t(rx+((hw+gap)*i), ry), hw, hh),
        curve_colors[i], background_color, border_color, grid_div_x, grid_div_y,
        grid_color, line_width
      );
    }
    return *this;
  }

protected:

  matv_t histograms_;
  bool is_hsv_;

};

template <typename crd_t, typename val_t, typename float_t>
std::ostream& operator<<(std::ostream& os, basic_histogram<crd_t,val_t,float_t> h)
{
  os << "{" << std::endl
     << "  bins: " << h.num_bins() << "," << std::endl
     << "  channels : " << h.mats().size() << std::endl
     << "}";
  return os;
}

// </editor-fold>

// <editor-fold desc="basic_camera_capture" defaultstate="collapsed">
template <typename Image_Type, typename Channel_Desc_Type=int>
class basic_camera_capture
{
public:

  /**
   * Types
   */
  typedef Image_Type imgage_t;
  typedef Channel_Desc_Type cam_selelct_t;

public:

  /**
   * Constructors, destructor
   */
  inline basic_camera_capture()
  { ; }

  virtual ~basic_camera_capture()
  { close(); }

private:

  inline basic_camera_capture(const basic_camera_capture & c) : cam_()
  { ; }

  inline void operator=(const basic_camera_capture& c)
  { throw "basic_camera_capture: No assignments, please!";  }

public:

  /**
   * Open the camera channel, returns success
   * @param cam_selelct_t
   * @return bool
   */
  bool open(cam_selelct_t channel)
  { close(); try { return cam_.open(channel) && cam_.isOpened(); } catch(...) { return false; } }

  /**
   * Close camera channel
   */
  void close() throw()
  { try { if(cam_.isOpened()) cam_.release(); } catch(...) {;} }

  /**
   * Take a snapshot, return success
   * @param imgage_t & img
   * @return bool
   */
  bool snapshot(imgage_t & img)
  {
    if(!cam_.read((::cv::Mat&)img)) { img.clear(); return false; }
    if(!img.mat().isContinuous()) img.mat() = img.mat().clone(); // No fun if that is not given.
    return true;
  }

private:

  cv::VideoCapture cam_;
};
// </editor-fold>

}}}

// <editor-fold desc="default specialisations" defaultstate="collapsed">
namespace sw { namespace im {
  typedef detail::basic_image<int, unsigned char, double> cvimage;
  typedef detail::basic_camera_capture<cvimage, int> camera_capture;
  typedef cvimage::exception image_exception;
  typedef cvimage::color_t color;
  typedef cvimage::point_t point2d;
  typedef cvimage::line_t line2d;
  typedef cvimage::line_list_t lines2d;
  typedef cvimage::rect_t rect2d;
  typedef cvimage::circle_t circle2d;
  typedef cvimage::circle_list_t circles2d;
  typedef cvimage::ellipse_t ellipse2d;
  typedef cvimage::rotated_rect_t rotated_rect2d;
  typedef cvimage::polynom_t polynom2d;
  typedef cvimage::contour_t contour2d;
  typedef cvimage::contour_list_t contours2d;
  typedef cvimage::contour_nested_list_t nested_contours2d;
}}
// </editor-fold>

// <editor-fold desc="undefs" defaultstate="collapsed">
#undef ex
#undef dbg
#undef dbg2
#endif
// </editor-fold>