roi.hpp

/*
 * roi.hpp
 *
 *     Project: jafar
 *  Created on: Jun 1, 2010
 *      Author: jsola
 */

#ifndef ROI_HPP_
#define ROI_HPP_

#include <cmath>
#include <boost/shared_ptr.hpp>

#include "jmath/jblas.hpp"

#include "image/Image.hpp"
#include "image/imageException.hpp"

#define USE_TEMPLATE_INLINE 0

namespace jafar {
  namespace image {
    using namespace jblas;

//    class RoiAbstract;
//    typedef boost::shared_ptr<RoiAbstract> roi_ptr_t;
/*
Well well... so many solutions...
Do we want a very general ROI, or is convex ROI enough ? with a convex roi some operations can fail (remove, add...)
For a general ROI, do we use a mask (meaning we have to test every pixel), or a slice model ?
For a slice model, do we use a vector of list, just a list, both ?
Do we use std::list, or a custom list which stores its data in an array to optimize memory management ?

Let's start with simple convex roi, while  thinking to a general roi




*/

    class ConvexRoi;
    class SliceRoi;

//#define NEW_ROI

    struct Pixel
    {
      int x, y;
     private:
      int x_end;
      friend class ConvexRoi;
      friend class SliceRoi;
    };
    
    struct Slice
    {
      int y, x1, x2;
     private:
      void* next;
      friend class ConvexRoi;
      friend class SliceRoi;
    };


    class ConvexRoi
    {
      public:
        struct SliceImpl { int x, w; };
        typedef std::vector<SliceImpl> ConvexRoiImpl;
      private:
        cv::Rect rect; // bounding box
        ConvexRoiImpl convex; // convex area, with origin at rect.x,rect.y
        bool hasConvex;
        int count_;
        mutable int y_end;
      private:
        void initConvex()
        {
          convex.resize(rect.height);
          for(ConvexRoiImpl::iterator it = convex.begin(); it != convex.end(); ++it)
          {
            it->x = 0;
            it->w = rect.width;
          }
          hasConvex = true;
        }
      public:
        // constructors
        ConvexRoi(): rect(cv::Rect(0,0,0,0)), hasConvex(false), count_(0) {}
        ConvexRoi(cv::Rect &rect): rect(rect), convex(), hasConvex(false), count_(rect.width*rect.height) {}
        ConvexRoi(int x, int y, int width, int height): rect(x,y,width,height), convex(), hasConvex(false), count_(width*height) {}
        ConvexRoi(jblas::vec2 x, jblas::sym_mat22 P, double scale);
        
        void init(cv::Rect _rect) { rect = _rect; convex.resize(0); hasConvex = false; count_ = rect.width*rect.height; }
        
        // rect interface
        int y(     ) const { return rect.y; }                                  
        int h(     ) const { return rect.height; }                             
        int x(     ) const { return rect.x; }                                  
        int w(     ) const { return rect.width; }                              

        // convex interface
        private:
        int x(int i) const { return rect.x + (hasConvex ? convex[i-rect.y].x : 0); }  
        int w(int i) const { return hasConvex ? convex[i-rect.y].w : rect.width; }    
        public:

        // generic interface, at slice level
        void start_slice(Slice &s) const { s.y = rect.y-1; y_end = y()+h(); }
        bool next_slice(Slice &s) const
        {
          ++s.y;
          if (s.y >= y_end) return false;
          s.x1 = x(s.y);
          s.x2 = x(s.y)+w(s.y)-1;
          return true;
        }

        // most generic interface, at individual pixels level, but harder to use because loses all line consistency
        void start(Pixel &pix) const { pix.y = y(); pix.x = x(pix.y)-1; pix.x_end = x(pix.y)+w(pix.y); y_end = y()+h(); }
        bool next(Pixel &pix) const
        {
          ++pix.x;
          if (pix.x >= pix.x_end)
          {
            ++pix.y;
            if (pix.y >= y_end) return false;
            pix.x = x(pix.y);
            pix.x_end = x(pix.y)+w(pix.y);
          }
          return true;
        }

        // more generic accessors
        int size(  ) const { return w()*h(); }                                 
        int count( ) const { return count_; }                                  
        bool isIn(jblas::vec2 p) const
        {
          int py = std::floor(p(1));
          int px = std::floor(p(0));
          return (py >= y() && py < y()+h() && px >= x(py) && px < x(py)+w(py));
        }

        ConvexRoi  operator+ (const ConvexRoi &roi) const 
        {
          ConvexRoi res(*this);
          // TODO
          JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          return res;
        }
        
        ConvexRoi& operator+=(const ConvexRoi &roi)
        {
          // TODO
          JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          return *this;
        }
        
        ConvexRoi  operator- (const ConvexRoi &roi) const
        {
          ConvexRoi res(*this);
          // TODO
          JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          return res;
        }
        
        ConvexRoi& operator-=(const ConvexRoi &roi)
        {
          // TODO
          JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          return *this;
        }
        
        ConvexRoi  operator* (const ConvexRoi &roi) const
        {
          // first deal with the rects
          cv::Rect newRect;
          newRect.x = std::max(rect.x, roi.rect.x);
          int x2 = std::min(rect.x+rect.width-1, roi.rect.x+roi.rect.width-1);
          newRect.width = std::max(0, x2-newRect.x+1);
          newRect.y = std::max(rect.y, roi.rect.y);
          int y2 = std::min(rect.y+rect.height-1, roi.rect.y+roi.rect.height-1);
          newRect.height = std::max(0, y2-newRect.y+1);
          // now with convex
          if (!hasConvex && !roi.hasConvex)
          {
            return ConvexRoi(newRect);
          } else
          if (hasConvex && roi.hasConvex)
          {
            // TODO
            JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          } else
          {
            ConvexRoi newRoi(newRect);
            newRoi.initConvex();
            newRoi.count_ = 0;
            
            const cv::Rect &oldRect = (hasConvex ? roi.rect : rect);
            const ConvexRoiImpl &oldConvex = (hasConvex ? convex : roi.convex);
            const cv::Rect &oldConvexRect = (hasConvex ? rect : roi.rect);
            
            int shift = std::max(0, oldRect.y - oldConvexRect.y);
            int lastY = -1000000;
            int minX1 = 1000000, maxX2 = -1000000;
//JFR_DEBUG("newRect begin " << newRect);
            for(int y = 0; y < newRect.height; ++y)
            {
              SliceImpl &newSlice = newRoi.convex[y];
              SliceImpl const &oldSlice = oldConvex[y+shift];
              newSlice.x = std::max(newRect.x, oldConvexRect.x+oldSlice.x)-newRect.x;
              int newX2 = std::min(newRect.x+newRect.width-1, oldConvexRect.x+oldSlice.x+oldSlice.w-1);
              newSlice.w = newX2-(newRect.x+newSlice.x)+1;
//JFR_DEBUG("slice y=" << y << " " << newSlice);
              if (newSlice.x >= newRect.width || newSlice.w <= 0)
              {
                if (lastY == -1000000)
                {
                  newRect.y++; newRect.height--;
                  shift++;
                  y--;
                  continue;
                } else
                {
                  newSlice.x = 0;
                  newSlice.w = 0;
                }
//JFR_DEBUG("error branch: newRect " << newRect << ", shift " << shift);
              } else
              {
                lastY = y;
                newRoi.count_ += newSlice.w;
                if (newSlice.x < minX1) minX1 = newSlice.x;
                if (newX2 > maxX2) maxX2 = newX2;
//JFR_DEBUG("no error branch: lastY " << lastY << ", minX1 " << minX1 << ", maxX2 " << maxX2);
              }
            }
            if (lastY == -1000000)
            {
              newRoi.rect.height = 0;
              newRoi.convex.resize(0);
            } else
            {
              newRoi.rect.height = lastY+1;
              newRoi.rect.width = maxX2-minX1-newRect.x+1;
              newRoi.rect.x += minX1;
              newRoi.convex.resize(newRoi.rect.height);
              if (minX1 != 0) for(int y = 0; y < newRect.height; ++y) newRoi.convex[y].x -= minX1;
            }
//JFR_DEBUG("newRect end " << newRect);
//JFR_DEBUG("newRoi.rect end " << newRoi.rect);
            
            newRoi.checkIntegrity();
            return newRoi;
          }
        }
        
        ConvexRoi& operator*=(const ConvexRoi &roi)
        {
          *this = *this * roi;
          return *this;
        }

        ConvexRoi  operator+ (const cv::Point &p) const
        { 
          ConvexRoi res(*this);
          res.rect.x += p.x;
          res.rect.y += p.y;
          return res;
        }
        ConvexRoi  operator- (const cv::Point &p) const
        { 
          ConvexRoi res(*this);
          res.rect.x -= p.x;
          res.rect.y -= p.y;
          return res;
        }
        
        ConvexRoi& operator+=(const cv::Point &p)
        {
          rect.x += p.x;
          rect.y += p.y;
          return *this;
        }
        ConvexRoi& operator-=(const cv::Point &p)
        {
          rect.x -= p.x;
          rect.y -= p.y;
          return *this;
        }
        
        void scaleTo(cv::Rect &newRect)
        {
// JFR_DEBUG("# scaleTo: oldRect " << rect << ", newRect " << newRect);
          cv::Point desiredOr(newRect.x, newRect.y);
          if (hasConvex)
          {
            // init new convex
            ConvexRoi newRoi(newRect);
            newRoi.initConvex();
            ConvexRoiImpl &newConvex = newRoi.convex;
// JFR_DEBUG("convex.size " << newConvex.size());
            double fw = newRect.width / (double)rect.width;
            newRect.x=std::floor(rect.x*fw);
// JFR_DEBUG("width fw " << fw);
            // first scale x
            for(ConvexRoiImpl::iterator it = convex.begin(); it != convex.end(); ++it)
            {
// JFR_DEBUG_BEGIN(); JFR_DEBUG_SEND("oldSlice " << *it);
/*              int x1 = std::floor((rect.x + it->x)*fw);
              int x2 = std::ceil((rect.x + it->x + it->w)*fw - 1);
              it->x = x1 - newRect.x;
              it->w = x2 - x1 + 1;
*/
              it->x = std::floor((rect.x + it->x)*fw) - newRect.x;
              it->w = std::ceil(it->w*fw);
              if (it->x+it->w > newRect.width) it->w = newRect.width-it->x;
// JFR_DEBUG_SEND(", newSlice x " << *it); JFR_DEBUG_END();
            }
            // then scale y
            int i = 0;
            newRoi.count_ = 0;
            double fh = newRect.height / (double)rect.height;
            newRect.y=std::floor(rect.y*fh);
            for(ConvexRoiImpl::iterator it = newConvex.begin(); it != newConvex.end(); ++it, ++i)
            {
              // same FIXME than in ellipse constructor
              //double jj = (i < newRect.height/2 ? (i+1)/fh : i/fh);
              int jj1 = std::floor(i/fh), jj2 = std::ceil((i+1)/fh);
              if (jj2 >= rect.height) jj2 = rect.height-1;
              int x1 = convex[jj1].x, x2 = convex[jj1].x+convex[jj1].w-1;
              for (int j = jj1+1; j <= jj2; ++j)
              {
                if (convex[j].x < x1) x1 = convex[j].x;
                int x2j = convex[j].x+convex[j].w-1;
                if (x2j > x2) x2 = x2j;
              }
              it->x = x1;
              it->w = x2-x1+1;
              newRoi.count_ += it->w;
// JFR_DEBUG("newSlice y " << *it << " by taking extrema between " << jj1 << " and " << jj2);
            }
#ifndef JFR_NDEBUG
            newRoi.rect = newRect;
            newRoi.checkIntegrity();
#endif
            *this = newRoi;
          } else
          {
            count_ = newRect.width*newRect.height;
          }
          rect = newRect;
          if (desiredOr.x >= 0) rect.x = desiredOr.x;
          if (desiredOr.y >= 0) rect.y = desiredOr.y;
          checkIntegrity();
        }
        
        
        ConvexRoi  operator* (double f) const
        {
          ConvexRoi res(*this);
          res *= f;
          return res; 
        }
        
        ConvexRoi& operator*=(double f)
        {
          cv::Rect newRect(-1,-1, std::ceil(rect.width *f), std::ceil(rect.height*f));
          scaleTo(newRect);
          return *this;
        }
        
        ConvexRoi& scale(double f)
        {
          cv::Point center(rect.x+rect.width/2, rect.y+rect.height/2);
          cv::Rect newRect(-1,-1, std::ceil(rect.width *f), std::ceil(rect.height*f));
          scaleTo(newRect);
          cv::Point newcenter(rect.x+rect.width/2, rect.y+rect.height/2);
          cv::Point shift = center-newcenter;
          *this += shift;
          return *this;
        }
        
        
      void draw()
      {
        std::cout << "rect: " << rect << std::endl;
        for(int i = 0; i < rect.height; ++i)
        {
          int x1 = convex[i].x;
          int x2 = convex[i].x+convex[i].w-1;
          for(int j = 0; j < x1; ++j) std::cout << ".";
          for(int j = x1; j <= x2; ++j) std::cout << "#";
          for(int j = x2+1; j < rect.width; ++j) std::cout << ".";
          std::cout << std::endl;
        }
      }
      
      void checkIntegrity()
      {
#ifndef JFR_NDEBUG
        JFR_ASSERT(rect.x > -1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.x < 1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.y > -1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.y < 1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.width >= 0, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.width < 1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.height >= 0, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.height < 1000000, "ConvexRoi::checkIntegrity failed");
        
        if (hasConvex)
        {
          JFR_ASSERT((int)convex.size() == rect.height, "ConvexRoi::checkIntegrity failed");
          int minX1 = 2000000000, maxX2 = -2000000000;
          int recount = 0;
          for(ConvexRoiImpl::iterator it = convex.begin(); it != convex.end(); ++it)
          {
            JFR_ASSERT(std::abs(it->x) < 2000000000 && std::abs(it->w) < 2000000000, "ConvexRoi::checkIntegrity failed")
            if (it->x < minX1) minX1 = it->x;
            int x2 = it->x+it->w-1;
            if (x2 > maxX2) maxX2 = x2;
            recount += it->w;
          }
          if (rect.height > 0 && rect.width > 0)
          {
            JFR_ASSERT(minX1 == 0,"ConvexRoi::checkIntegrity failed");
            JFR_ASSERT(maxX2 == rect.width-1, "ConvexRoi::checkIntegrity failed");
          }
          JFR_ASSERT(recount == count_, "ConvexRoi::checkIntegrity failed");
        }
#endif
      }

      friend std::ostream& operator<<(std::ostream &os, ConvexRoi::SliceImpl const &c);
      friend std::ostream& operator<<(std::ostream& os, ConvexRoi const& roi);
      friend class SliceRoi;

    };
    
    std::ostream& operator<<(std::ostream &os, ConvexRoi::SliceImpl const &c);
    std::ostream& operator<<(std::ostream& os, ConvexRoi const& roi);


#if 1

    class SliceRoi
    {
      public:
        struct SliceImpl
        {
          int x, w; SliceImpl *next;
          SliceImpl(int x, int w): x(x), w(w), next(NULL) {}
          SliceImpl(): x(-1), w(-1), next(NULL) {}
        };
        typedef std::vector<SliceImpl*> SliceRoiImpl;
      private:
        cv::Rect rect; // bounding box
        mutable SliceRoiImpl slices; // area, with origin at rect.x,rect.y
        mutable bool hasSlices;
        int count_;
        mutable int y_end;
      private:
        void initSlices(bool full) const
        {
          slices.resize(rect.height);
          for(SliceRoiImpl::iterator it = slices.begin(); it != slices.end(); ++it)
            *it = full ? new SliceImpl(0,rect.width) : NULL;
          if (full) hasSlices = true;
        }
        void delete_slice_line(SliceImpl ** slice_start)
        {
          SliceImpl *next;
          for(SliceImpl * slice = *slice_start; slice; slice = next)
          {
            next = slice->next;
            delete slice;
          }
          *slice_start = NULL;
        }
        void delete_slices()
        {
          for(SliceRoiImpl::iterator it = slices.begin(); it != slices.end(); ++it) delete_slice_line(&(*it));
          slices.resize(0);
          hasSlices = false;
        }
        static unsigned merge_line(SliceImpl **newSlice, SliceImpl const * slice1, SliceImpl const * slice2)
        {
          JFR_ASSERT(*newSlice == NULL, "*newSlice must be NULL");

          unsigned linecount = 0;
          int right = -1; // first outside at the right extremity
          while (slice1 || slice2)
          {
            // continuity broken, create new slice
            if ((!slice1 || slice1->x > right) && (!slice2 || slice2->x > right))
            {
              if (*newSlice) // if not first time
              {
                (*newSlice)->w = right - (*newSlice)->x;
                linecount += (*newSlice)->w;
                newSlice = &((*newSlice)->next);
              }

              if (slice1 && (!slice2 || slice1->x < slice2->x))
              {
                *newSlice = new SliceImpl(slice1->x, slice1->w);
                right = slice1->x + slice1->w;
                slice1 = slice1->next;
              }
              if (slice2 && (!slice1 || slice2->x <= slice1->x))
              {
                *newSlice = new SliceImpl(slice2->x, slice2->w);
                right = slice2->x + slice2->w;
                slice2 = slice2->next;
              }
            }

            // merge slice1
            if (slice1 && slice1->x <= right)
            {
              if (slice1->x+slice1->w > right) right = slice1->x+slice1->w;
              slice1 = slice1->next;
            }

            // merge slice2
            if (slice2 && slice2->x <= right)
            {
              if (slice2->x+slice2->w > right) right = slice2->x+slice2->w;
              slice2 = slice2->next;
            }

          }
          (*newSlice)->w = right - (*newSlice)->x;
          linecount += (*newSlice)->w;
          return linecount;
        }
      public:
        // constructors
        SliceRoi(): rect(cv::Rect(0,0,0,0)), hasSlices(false), count_(0) {}
        SliceRoi(cv::Rect &rect): rect(rect), slices(), hasSlices(false), count_(rect.width*rect.height) {}
        SliceRoi(int x, int y, int width, int height): rect(x,y,width,height), slices(), hasSlices(false), count_(width*height) {}
        SliceRoi(jblas::vec2 x, jblas::sym_mat22 P, double scale) { *this = ConvexRoi(x, P, scale); }
        SliceRoi(const SliceRoi & roi) { *this = roi; }
        SliceRoi & operator=(const SliceRoi & roi);
        SliceRoi(const ConvexRoi & roi) { *this = roi; }
        SliceRoi & operator=(const ConvexRoi & croi);
        ~SliceRoi() { delete_slices(); }

        void init(cv::Rect _rect) { rect = _rect; slices.resize(0); hasSlices = false; count_ = rect.width*rect.height; }

        // rect interface
        int y(     ) const { return rect.y; }                                  
        int h(     ) const { return rect.height; }                             
        int x(     ) const { return rect.x; }                                  
        int w(     ) const { return rect.width; }                              

        // convex interface
        private:
        int x(int i) const { return rect.x + (hasSlices ? slices[i-rect.y]->x : 0); }  
        int w(int i) const 
        {
          if (!hasSlices) return rect.width;
          SliceImpl *slice = slices[i-rect.y];
          while (slice->next != NULL) slice = slice->next;
          return slice->x + slice->w - slices[i-rect.y]->x;
        }
        public:

        // generic interface, at slice level
        void start_slice(Slice &s) const { s.y = rect.y-1; s.next = NULL; y_end = y()+h(); }
        bool next_slice(Slice &s) const
        {
          SliceImpl *slice = (SliceImpl*)s.next;
          if (slice == NULL) // go to next line
          {
            ++s.y;
            if (s.y >= y_end) return false;
            slice = slices[s.y-rect.y];
          }
          s.x1 = rect.x+slice->x;
          s.x2 = s.x1 + slice->w - 1;
          s.next = slice->next;
          return true;
        }

        // most generic interface, at individual pixels level, but harder to use because loses all line consistency
/*        void start(Pixel &pix) const { pix.y = y(); pix.x = x(pix.y)-1; pix.x_end = x(pix.y)+w(pix.y); y_end = y()+h(); }
        bool next(Pixel &pix) const
        {
          ++pix.x;
          if (pix.x >= pix.x_end)
          {
            ++pix.y;
            if (pix.y >= y_end) return false;
            pix.x = x(pix.y);
            pix.x_end = x(pix.y)+w(pix.y);
          }
          return true;
        }
*/
        // more generic accessors
        int size(  ) const { return w()*h(); }                                 
        int count( ) const { return count_; }                                  
        bool isIn(jblas::vec2 p) const
        {
          int py = std::floor(p(1))-rect.y;
          int px = std::floor(p(0))-rect.x;
          if (py < y() || py >= y()+h()) return false;
          for(SliceImpl *slice = slices[py]; slice != NULL; slice = slice->next)
            if (px >= slice->x && px < slice->x+slice->w) return true;
          return false;
        }


        SliceRoi  operator+ (const SliceRoi &roi) const
        {
          cv::Rect newRect;
          newRect.x = std::min(rect.x, roi.rect.x);
          newRect.y = std::min(rect.y, roi.rect.y);
          newRect.width = std::max(rect.x+rect.width-1, roi.rect.x+roi.rect.width-1) - newRect.x + 1;
          newRect.height = std::max(rect.y+rect.height-1, roi.rect.y+roi.rect.height-1) - newRect.y + 1;

          SliceRoi newRoi(newRect);
          newRoi.initSlices(false);
          newRoi.count_ = 0;

          if (!hasSlices) initSlices(true);
          if (!roi.hasSlices) roi.initSlices(true);

          for(int i = 0; i < newRect.height; ++i)
          {
            SliceImpl **newSlice = &(newRoi.slices[i]);
            SliceImpl const * slice1 = ((i < rect.y-newRect.y || i > rect.y+rect.height-1-newRect.y) ? NULL : slices[i-rect.y-newRect.y]);
            SliceImpl const * slice2 = ((i < roi.rect.y-newRect.y || i > roi.rect.y+roi.rect.height-1-newRect.y) ? NULL : roi.slices[i]);

            newRoi.count_ += merge_line(newSlice, slice1, slice2);
          }

          return newRoi;
        }

        SliceRoi& operator+=(const SliceRoi &roi)
        {
          *this = *this + roi;
          return *this;
        }

        SliceRoi  operator- (const SliceRoi &roi) const
        {
          SliceRoi res(*this);
          // TODO
          JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          return res;
        }

        SliceRoi& operator-=(const SliceRoi &roi)
        {
          // TODO
          JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          return *this;
        }

        SliceRoi  operator* (const SliceRoi &roi) const
        {
          // first deal with the rects
          cv::Rect newRect;
          newRect.x = std::max(rect.x, roi.rect.x);
          int x2 = std::min(rect.x+rect.width-1, roi.rect.x+roi.rect.width-1);
          newRect.width = std::max(0, x2-newRect.x+1);
          newRect.y = std::max(rect.y, roi.rect.y);
          int y2 = std::min(rect.y+rect.height-1, roi.rect.y+roi.rect.height-1);
          newRect.height = std::max(0, y2-newRect.y+1);

          // now with slices
          if (!hasSlices && !roi.hasSlices)
          {
            return SliceRoi(newRect);
          } else
          if (hasSlices && roi.hasSlices)
          {
            // TODO
            JFR_ERROR(ImageException, ImageException::NOT_IMPLEMENTED, "this operator hasn't been implemented yet");
          } else
          {
            SliceRoi newRoi(newRect);
            newRoi.initSlices(false);
            newRoi.count_ = 0;

            const cv::Rect &oldRect = (hasSlices ? roi.rect : rect);
            const SliceRoiImpl &oldSlices = (hasSlices ? slices : roi.slices);
            const cv::Rect &oldSlicesRect = (hasSlices ? rect : roi.rect);

            int shiftY = std::max(0, oldRect.y - oldSlicesRect.y);
            int shiftX = std::max(0, oldRect.x - oldSlicesRect.x);
            int minX1 = 1000000, maxX2 = -1000000;
//JFR_DEBUG("newRect begin " << newRect);

            //-- copy slices
            for(int y = 0; y < newRect.height; ++y)
            {
              SliceImpl **newSlice = &(newRoi.slices[y]);
              SliceImpl const * oldSlice = oldSlices[y+shiftY];
              for(; oldSlice != NULL; oldSlice = oldSlice->next)
              {
                if (oldSlice->x < shiftX) // start before of new rect
                {
                  if (oldSlice->x+oldSlice->w-1 >= shiftX); // ends inside or after new rect
                  { // copy and truncate beginning
                    (*newSlice) = new SliceImpl(0, std::min(newRect.width, oldSlice->w - shiftX));
                    if (0 < minX1) minX1 = 0;
                    int x2 = (*newSlice)->x+(*newSlice)->w-1; if (x2 > maxX2) maxX2 = x2;
                    newRoi.count_ += (*newSlice)->w;
                    newSlice = &((*newSlice)->next);
                  } // otherwise do nothing, throw it away it's out
                }
                else
                {
                  if (oldSlice->x-shiftX < newRect.width) // start inside new rect
                  { // copy and truncate end
                    (*newSlice) = new SliceImpl(oldSlice->x-shiftX, std::min(newRect.width-(oldSlice->x-shiftX), oldSlice->w));
                    if ((*newSlice)->x < minX1) minX1 = (*newSlice)->x;
                    int x2 = (*newSlice)->x+(*newSlice)->w-1; if (x2 > maxX2) maxX2 = x2;
                    newRoi.count_ += (*newSlice)->w;
                    newSlice = &((*newSlice)->next);
                  } // otherwise do nothing, throw it away it's out
                }
              }
            }

            //-- update rect
            int nremove;
            // top
            for(nremove = 0; newRoi.slices[nremove] == NULL && nremove < newRect.height; nremove++);
            newRect.y += nremove; newRect.height -= nremove;
            for(int i = 0; i < newRect.height; ++i) newRoi.slices[i] = newRoi.slices[i+nremove]; newRoi.slices.resize(newRect.height);
            // bottom
            for(nremove = 0; newRoi.slices[newRect.height-1-nremove] == NULL && nremove < newRect.height; nremove++);
            newRect.height -= nremove;
            newRoi.slices.resize(newRect.height);
            // left
            if (minX1 > 0)
            {
              newRect.x += minX1; newRect.width -= minX1;
              for(int i = 0; i < newRect.height; ++i) for(SliceImpl *slice = newRoi.slices[i]; slice; slice = slice->next) slice->x -= minX1;
              maxX2 -= minX1;
            }
            // right
            maxX2 = newRect.width-1-maxX2; // this is now the space to remove at the end
            if (maxX2 > 0) newRect.width -= maxX2;

//JFR_DEBUG("newRect end " << newRect);
//JFR_DEBUG("newRoi.rect end " << newRoi.rect);

            newRoi.checkIntegrity();
            return newRoi;
          }
        }

        SliceRoi& operator*=(const SliceRoi &roi)
        {
          *this = *this * roi;
          return *this;
        }

        SliceRoi  operator+ (const cv::Point &p) const
        {
          SliceRoi res(*this);
          res.rect.x += p.x;
          res.rect.y += p.y;
          return res;
        }
        SliceRoi  operator- (const cv::Point &p) const
        {
          SliceRoi res(*this);
          res.rect.x -= p.x;
          res.rect.y -= p.y;
          return res;
        }

        SliceRoi& operator+=(const cv::Point &p)
        {
          rect.x += p.x;
          rect.y += p.y;
          return *this;
        }
        SliceRoi& operator-=(const cv::Point &p)
        {
          rect.x -= p.x;
          rect.y -= p.y;
          return *this;
        }

        // FIXME apparently only works for reducing...
        void scaleTo(cv::Rect &newRect)
        {
// JFR_DEBUG("# scaleTo: oldRect " << rect << ", newRect " << newRect);
          cv::Point desiredOr(newRect.x, newRect.y); // desired origin ; FIXME do we really need this ?
          if (hasSlices)
          {
            // init new slices
            SliceRoi newRoi(newRect);
            newRoi.initSlices(false);
            SliceRoiImpl &newSlices = newRoi.slices;
// JFR_DEBUG("convex.size " << newSlices.size());
            double fw = newRect.width / (double)rect.width; // width factor
            newRect.x = std::floor(rect.x*fw);
// JFR_DEBUG("width fw " << fw);

            // first scale x, easy
            for(SliceRoiImpl::iterator it = slices.begin(); it != slices.end(); ++it)
            {
              for(SliceImpl *slice = *it; slice != NULL; slice = slice->next)
              {
// JFR_DEBUG_BEGIN(); JFR_DEBUG_SEND("oldSlice " << slice);
                slice->x = std::floor((rect.x + slice->x)*fw - newRect.x);
                slice->w = std::ceil(slice->w*fw);
                if (slice->x+slice->w > newRect.width) slice->w = newRect.width-slice->x;
// JFR_DEBUG_SEND(", newSlice x " << slice); JFR_DEBUG_END();
              }
            }

            // then scale y, harder with interpolation
            // the idea when shrinking is to keep the widest area, so we merge the different lines
            int i = 0;
            newRoi.count_ = 0;
            double fh = newRect.height / (double)rect.height; // height factor
            newRect.y = std::floor(rect.y*fh);
            for(SliceRoiImpl::iterator it = newSlices.begin(); it != newSlices.end(); ++it, ++i)
            {
              // same FIXME than in ellipse constructor
              int jj1 = std::floor(i/fh), jj2 = std::ceil((i+1)/fh); // factor of how many lines to add in between each pair of lines
              if (jj2 >= rect.height) jj2 = rect.height-1;
              SliceImpl *merged_line[2] = { NULL, NULL};
              int midx = 0;
              for (int j = jj1+1; j <= jj2; ++j)
              {
                merge_line(&(merged_line[1-midx]), merged_line[midx], slices[j]);
                delete_slice_line(&(merged_line[midx]));
                midx = 1-midx;
              }
              *it =  merged_line[midx];
              for(SliceImpl *slice = *it; slice; slice = slice->next) newRoi.count_ += slice->w;
// JFR_DEBUG("newSlice y " << *it << " by taking extrema between " << jj1 << " and " << jj2);
            }
#ifndef JFR_NDEBUG
            newRoi.rect = newRect;
            newRoi.checkIntegrity();
#endif
            *this = newRoi;
          } else
          {
            count_ = newRect.width*newRect.height;
          }
          rect = newRect;
          if (desiredOr.x >= 0) rect.x = desiredOr.x;
          if (desiredOr.y >= 0) rect.y = desiredOr.y;
          checkIntegrity();
        }


        SliceRoi  operator* (double f) const
        {
          SliceRoi res(*this);
          res *= f;
          return res;
        }

        SliceRoi& operator*=(double f)
        {
          cv::Rect newRect(-1,-1, std::ceil(rect.width *f), std::ceil(rect.height*f));
          scaleTo(newRect);
          return *this;
        }

        SliceRoi& scale(double f)
        {
          cv::Point center(rect.x+rect.width/2, rect.y+rect.height/2);
          cv::Rect newRect(-1,-1, std::ceil(rect.width *f), std::ceil(rect.height*f));
          scaleTo(newRect);
          cv::Point newcenter(rect.x+rect.width/2, rect.y+rect.height/2);
          cv::Point shift = center-newcenter;
          *this += shift;
          return *this;
        }


      void draw()
      {
        std::cout << "rect: " << rect << std::endl;
        for(int i = 0; i < rect.height; ++i)
        {
          int x0 = 0;
          for(SliceImpl *slice = slices[i]; slice != NULL; slice = slice->next)
          {
            int x1 = slice->x;
            int x2 = slice->x+slice->w-1;
            for(int j = x0; j < x1; ++j) std::cout << ".";
            for(int j = x1; j <= x2; ++j) std::cout << "#";
            x0 = x2+1;
          }
          for(int j = x0+1; j < rect.width; ++j) std::cout << ".";
          std::cout << std::endl;
        }
      }

      void checkIntegrity()
      {
#ifndef JFR_NDEBUG
        JFR_ASSERT(rect.x > -1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.x < 1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.y > -1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.y < 1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.width >= 0, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.width < 1000000, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.height >= 0, "ConvexRoi::checkIntegrity failed");
        JFR_ASSERT(rect.height < 1000000, "ConvexRoi::checkIntegrity failed");

        if (hasSlices)
        {
          JFR_ASSERT((int)slices.size() == rect.height, "ConvexRoi::checkIntegrity failed");
          int minX1 = 2000000000, maxX2 = -2000000000;
          int recount = 0;
          for(SliceRoiImpl::iterator it = slices.begin(); it != slices.end(); ++it)
          for(SliceImpl *slice = *it; slice; slice = slice->next)
          {
            JFR_ASSERT(std::abs(slice->x) < 2000000000 && std::abs(slice->w) < 2000000000, "ConvexRoi::checkIntegrity failed")
            if (slice->x < minX1) minX1 = slice->x;
            int x2 = slice->x+slice->w-1;
            if (x2 > maxX2) maxX2 = x2;
            recount += slice->w;
          }
          if (rect.height > 0 && rect.width > 0)
          {
            JFR_ASSERT(minX1 == 0,"ConvexRoi::checkIntegrity failed");
            JFR_ASSERT(maxX2 == rect.width-1, "ConvexRoi::checkIntegrity failed");
          }
          JFR_ASSERT(recount == count_, "ConvexRoi::checkIntegrity failed");
        }
#endif
      }

      friend std::ostream& operator<<(std::ostream &os, SliceRoi::SliceImpl const &c);
      friend std::ostream& operator<<(std::ostream& os, SliceRoi const& roi);
    };

    std::ostream& operator<<(std::ostream &os, SliceRoi::SliceImpl const &c);
    std::ostream& operator<<(std::ostream& os, SliceRoi const& roi);

#endif

#if 0
    
    struct ConvexRange { int y, x, w; };

    class GeneralRoi: public std::list<ConvexRange>
    {
      private:
        
        typedef std::vector<iterator> RoiLines; // size +1
      
      private:
        
        RoiLines lines;
        
      public:
        GeneralRoi() {}
      public:
        // constructors
        Roi(Roi const & roi)
        {
          rect = roi.rect;
          if (roi.convex) { convex = new ConvexRoi(); *convex = *(roi.convex); } else convex = NULL;
          if (roi.mask) { mask = new image::Image(); *mask = *(roi.mask); } else mask = NULL;
        }
        Roi& operator=(Roi const & roi)
        {
          rect = roi.rect;
          delete convex; if (roi.convex) { convex = new ConvexRoi(); *convex = *(roi.convex); } else convex = NULL;
          delete mask; if (roi.mask) { mask = new image::Image(); *mask = *(roi.mask); } else mask = NULL;
        }
        Roi(cv::Rect &rect): rect(rect), convex(NULL), mask(NULL) {}
        Roi(int x, int y, int width, int height): rect(x,y,width,height), convex(NULL), mask(NULL) {}
        Roi(jblas::vec2 &x, jblas::sym_mat22 &P, double factor)
        {
          // TODO
        }
        Roi(int x, int y, image::Image &mask_): rect(x,y,mask_.width(),mask_.height()), convex(new ConvexRoi()), mask(new image::Image(mask_)
        {
          convex->resize(rect.height);
          for(int i = 0; i < rect.height; i++)
          {
            int j;
//            for(j = 0; j < rect.width; j++) if (mask(
          // TODO
          }
        }
        ~Roi() { delete convex; delete mask; }
        
        iterator begin(int i) { return lines[i]; }
        iterator end(int i) { return lines[i+1]; }
        
        // accessors
        int top   () const { return rect.y; }
        int bottom() const { return rect.y+rect.height-1; }
        int left () const { return rect.x; }
        int right() const { return rect.x+rect.width-1; }
        int left (int i) const { return convex ? rect[i].first  : left() ; }
        int right(int i) const { return convex ? rect[i].second : right(); }
        bool test(int x, int y) const 
        {
          if (mask) return mask.at<uchar>(y-top(),x-left());
          if (convex) return y>=top() && y <= bottom() && x>=left(y-top()) && x<=right(y-top()); 
                 else return y>=top() && y <= bottom() && x>=left()        && x<=right(); 
        
        }
        double getRatio() { return 1; }
        
        // introspection
        bool isRect() const { return convex == NULL && mask == NULL; }
        bool isConvex() const { return mask == NULL; }
        
        // operations
        void scale(double factor) 
        { 
          rect.x=(int)(rect.x*factor);
          rect.y=(int)(rect.y*factor);
          rect.width =(int)(rect.width *factor+0.999);
          rect.height=(int)(rect.height*factor+0.999);
        }
        void shift(int x, int y) 
        { 
          rect.x += x; 
          rect.y += y; 
        }
        void add(Roi &roi) 
        {
          // TODO
        }
        void remove(Roi &roi) 
        {
          // TODO
        }
        void intersect(Roi &roi) 
        {
          // TODO
        }
    };
#endif


#if 0

    class RoiAbstract {
      
      public:
        virtual int top() const = 0;
        virtual int bottom() const = 0;
        virtual int left() const = 0; 
        virtual int right() const = 0; 
        virtual int left(int i) const = 0;
        virtual int right(int i) const = 0;
        virtual int width() const { return right()-left()+1; }
        virtual int height() const { return bottom()-top()+1; }
        virtual bool test(int x, int y) const = 0;
        virtual bool inTest(int x, int y) const = 0;
        virtual bool isConvex() const = 0;
        
        virtual RoiAbstract* clone() = 0;
        virtual void scale(double factor) = 0;
        virtual void shift(int x, int y) = 0;
//        virtual void remove(RoiAbstract &roi) = 0;

        std::ostream& operator<<(std::ostream& os)
        {
          os << "(" << left() << "," << top() << ";" << width() << "," << height() << ")";
          return os;
        }
    };
    

    class RoiRect
#if !USE_TEMPLATE_INLINE
      : public RoiAbstract
#endif
    {
      private:
        cv::Rect roi;
      public:
        RoiRect(cv::Rect &roi): roi(roi) {}
        RoiRect(int x, int y, int width, int height): roi(x,y,width,height) {}
        int top() const { return roi.y; }
        int bottom() const { return roi.y+roi.height-1; }
        int left() const { return roi.x; }
        int right() const { return roi.x+roi.width-1; }
        int left(int i) const { return left(); }
        int right(int i) const { return right(); }
        bool test(int x, int y) const { return y>=top() && y <= bottom() && x>=left() && x<=right(); }
        bool inTest(int x, int y) const { return true; }
        bool isConvex() const { return true; }
        RoiAbstract* clone() { return new RoiRect(roi); }
        void scale(double factor) { roi.x=(int)(roi.x*factor); roi.y=(int)(roi.y*factor); roi.width=(int)(roi.width*factor+0.999);  roi.height=(int)(roi.height*factor+0.999); }
        void shift(int x, int y) { roi.x += x; roi.y += y; }
    };
    
    class RoiEllipse
#if !USE_TEMPLATE_INLINE
      : public RoiAbstract
#endif
    {
      private:
        jblas::vec2 x;
        jblas::sym_mat22 P;
      public:
        RoiEllipse(jblas::vec2 &x, jblas::sym_mat22 &P) { /* TODO */ }
        int top() const { return 0;/* TODO */ }
        int bottom() const { return 0;/* TODO */ }
        int left() const { return 0;/* TODO */ }
        int right() const { return 0;/* TODO */ }
        int left(int i) const { return 0;/* TODO */ }
        int right(int i) const { return 0;/* TODO */ }
        bool test(int x, int y) const { return 0;/* TODO */ }
        bool inTest(int x, int y) const { return true; }
        bool isConvex() const { return true; }
    };
    
    class RoiConvex
#if !USE_TEMPLATE_INLINE
      : public RoiAbstract
#endif
    {
      private:
        int top_;
        int bottom_;
        std::vector<std::pair<int,int> > roi;
      public:
        RoiConvex(int maxHeight) { /* TODO */ }
        RoiConvex(jblas::vec2 &x, jblas::sym_mat22 &P) { /* TODO */ }
        RoiConvex(int top, int bottom, std::vector<int>) { /* TODO */ }
        void intFromEllipse(jblas::vec2 &x, jblas::sym_mat22 &P) { /* TODO */ }
        int top() const { return top_; }
        int bottom() const { return bottom_; }
        int left() const { return 0;/* TODO */ }
        int right() const { return 0;/* TODO */ }
        int left(int i) const { return roi[i].first; }
        int right(int i) const { return roi[i].second; }
        bool test(int x, int y) const { return y>=top_ && y<=bottom_ && x>=roi[y].first && x<=roi[y].second; }
        bool inTest(int x, int y) const { return true; }
        bool isConvex() const { return true; }
    };
    
    
    class RoiMask
#if !USE_TEMPLATE_INLINE
      : public RoiAbstract
#endif
    {
      private:
        cv::Rect roi;
        image::Image mask;
      public:
        RoiMask(int top, int left, image::Image &mask) { /* TODO */ }
        RoiMask(cv::Rect &roi, image::Image &mask) { /* TODO */ }
        int top() const { return roi.y; }
        int bottom() const { return roi.y+roi.height-1; }
        int left() const { return 0;/* TODO */ }
        int right() const { return 0;/* TODO */ }
        int left(int i) const { return roi.x; }
        int right(int i) const { return roi.x+roi.width-1; }
        bool test(int x, int y) const { return y>=top() && y <= bottom() && x>=left(y) && y<=right(y); }
        bool inTest(int x, int y) const { return mask.data(y)[x]; }
        bool isConvex() const { return false; }
    };

#endif 

    
    
    class ROI;
    typedef boost::shared_ptr<ROI> roi_ptr_t;
    
    class ROI {
      public:
        int x, y;
        int width, height;
        ROI(const int ulx = 0, const int uly = 0, const int sx = 0, const int sy = 0) {
          x = ulx;
          y = uly;
          width = sx;
          height = sy;
        }
        veci2 downright() {
          veci2 dr;
          dr(0) = x + width - 1;
          dr(1) = y + height - 1;
          return dr;
        }
        veci2 upleft() {
          veci2 ul;
          ul(0) = x;
          ul(1) = y;
          return ul;
        }
        veci2 size() {
          veci2 sz;
          sz(0) = width;
          sz(1) = height;
          return sz;
        }
        void upleft(const veci & ul){
          x = ul(0);
          y = ul(1);
        }
        void downright(const veci & dr){
          width = dr(0) - x + 1;
          height = dr(1) - y + 1;
        }
        void size(const veci & sz){
          width = sz(0);
          height = sz(1);
        }

        friend std::ostream& operator <<(std::ostream & s, ROI & roi) {
          s << "ROI: " << std::endl;
          s << " .upper left  = " << roi.upleft() << std::endl;
          s << " .lower right = " << roi.downright() << std::endl;
          s << " .size        = " << roi.size() << std::endl;
          return s;
        }

    };
    

  }
}


#endif /* ROI_HPP_ */