#ifndef _OTB_TRISKELE_ARRAY_TREE_GRAPHWALKER_TPP
#define _OTB_TRISKELE_ARRAY_TREE_GRAPHWALKER_TPP

// ========================================
inline DimImg
GraphWalker::vertexMaxCount () const {
  return vertexMaxCount (size);
}

// ----------------------------------------
inline DimImg
GraphWalker::vertexMaxCount (const Size &tileSize) const {
  return DimImg (tileSize.width)*DimImg (tileSize.height);
}

// ----------------------------------------
inline DimEdge
GraphWalker::edgeMaxCount () const {
  return edgeMaxCount (size);
}

// ----------------------------------------
inline DimEdge
GraphWalker::edgeMaxCount (const Size &tileSize) const {
  if (!tileSize.width || !tileSize.height)
    return 0;
  DimEdge miniSquare = DimImg (tileSize.width)*DimImg (tileSize.height-1) + DimImg (tileSize.width-1)*DimImg (tileSize.height);
  switch (connectivity) {
  case Connectivity::C4 :
    return miniSquare;
  case Connectivity::C6N :
  case Connectivity::C6P :
    return miniSquare + DimImg (tileSize.width-1)*DimImg (tileSize.height-1);
  default :
    BOOST_ASSERT (connectivity == Connectivity::C8);
    return miniSquare + 2*DimImg (tileSize.width-1)*DimImg (tileSize.height-1);
  };
}
// ----------------------------------------
inline DimEdge
GraphWalker::edgeBoundaryMaxCount (const DimSideImg &side) const {
  if (!side)
    return 0;
  switch (connectivity) {
  case Connectivity::C4 :
    return side;
  case Connectivity::C6N :
  case Connectivity::C6P :
    return side+side-1;
  default :
    BOOST_ASSERT (connectivity == Connectivity::C8);
    return 3*side-2;
  };
}

// ========================================
inline void
GraphWalker::setTiles (const unsigned int &coreCount, const Rect &tile,
		       std::vector<Rect> &tiles, std::vector<Rect> &boundaries, std::vector<bool> &verticalBoundaries) const {
  BOOST_ASSERT (coreCount);
  DEF_LOG ("GraphWalker::setTiles", "coreCount:" << coreCount << " tile:" << tile);
  if (coreCount < 2 || max (tile.width, tile.height) < 4 || min (tile.width, tile.height) < 3) {
    LOG ("tile:" << tile);
    tiles.push_back (tile);
    return;
  }
  bool vertical = tile.width > tile.height;
  bool odd = coreCount & 0x1;
  DimImg thin = (vertical ? tile.width : tile.height) / (odd ? coreCount : 2);
  Rect tileA (tile);
  Rect tileB (tile);
  if (vertical) {
    tileA.width = thin;
    tileB.x += thin;
    tileB.width -= thin;
  } else {
    tileA.height = thin;
    tileB.y += thin;
    tileB.height -= thin;
  }
  setTiles ((odd ? 1 : coreCount/2), tileA, tiles, boundaries, verticalBoundaries);
  boundaries.push_back (tileB);
  verticalBoundaries.push_back (vertical);
  setTiles ((odd ? coreCount-1 : coreCount/2), tileB, tiles, boundaries, verticalBoundaries);
}

// ========================================
template<typename Funct>
inline DimImg
GraphWalker::forEachVertexIdx (const Funct &lambda /*lambda (idx)*/) const {
  DimImg maxCount = vertexMaxCount ();
  DimImg vertexCount = 0;
  for (DimSideImg idx = 0; idx < maxCount; ++idx)
    GRAPHWALKER_CALL_LAMBDA_IDX (idx, vertexCount, lambda (idx));
  return vertexCount;
}

template<typename Funct>
inline DimImg
GraphWalker::forEachVertexIdx (const Rect &rect, const Funct &lambda /*lambda (idx)*/) const {
  DimImg vertexCount = 0;
  DimSideImg const lastX = rect.x+rect.width;
  DimSideImg const lastY = rect.y+rect.height;
  for (DimSideImg y = rect.y; y < lastY; ++y)
    for (DimSideImg x = rect.x; x < lastX; ++x) {
      const Point p (x, y);
      DimImg idx = point2idx (size, p);
      GRAPHWALKER_CALL_LAMBDA_IDX (idx, vertexCount, lambda (idx));
    }
  return vertexCount;
}

template<typename Funct>
inline DimImg
GraphWalker::forEachVertexPt (const Rect &rect, const Funct &lambda /*lambda (p)*/) const {
  DimImg vertexCount = 0;
  DimSideImg const lastX = rect.x+rect.width;
  DimSideImg const lastY = rect.y+rect.height;
  for (DimSideImg y = rect.y; y < lastY; ++y)
    for (DimSideImg x = rect.x; x < lastX; ++x) {
      const Point p (x, y);
      DimImg idx = point2idx (size, p);
      GRAPHWALKER_CALL_LAMBDA_IDX (idx, vertexCount, lambda (p));
    }
  return vertexCount;
}

// ========================================
template<typename Funct>
inline DimEdge
GraphWalker::forEachEdgePt (const Rect &rect, TileItem tileItem, const Funct &lambda /*lambda (p0, p1)*/) const {
  DimImg edgeCount = 0;
  if (!rect.width || !rect.height)
    return edgeCount;
  const DimSideImg firstX = rect.x, endX = rect.x+rect.width;
  const DimSideImg firstY = rect.y, endY = rect.y+rect.height;
  switch (tileItem) {

  case Surface:
    {
      for (DimSideImg x = firstX+1; x < endX; ++x) {
	Point pha1 (x-1, firstY), pha2 (x, firstY);
	GRAPHWALKER_CALL_LAMBDA_PTS (pha1, pha2, edgeCount, lambda (pha1, pha2));
      }
      for (DimSideImg y = firstY+1; y < endY; ++y) {
	Point pva1 (firstX, y-1), pva2 (firstX, y);
	GRAPHWALKER_CALL_LAMBDA_PTS (pva1, pva2, edgeCount, lambda (pva1, pva2));
	for (DimSideImg x = firstX+1; x < endX; ++x) {
	  Point ph1 (x-1, y), pv1 (x, y-1), p2 (x, y);
	  GRAPHWALKER_CALL_LAMBDA_PTS (ph1, p2, edgeCount, lambda (ph1, p2));
	  GRAPHWALKER_CALL_LAMBDA_PTS (pv1, p2, edgeCount, lambda (pv1, p2));
	}
      }
      if (connectivity & Connectivity::C6P)
	for (DimSideImg y = firstY+1; y < endY; ++y)
	  for (DimSideImg x = firstX+1; x < endX; ++x) {
	    Point p1 (x-1, y-1), p2 (x, y);
	    GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
	  }
      if (connectivity & Connectivity::C6N)
	for (DimSideImg y = firstY+1; y < endY; ++y)
	  for (DimSideImg x = firstX+1; x < endX; ++x) {
	    Point p1 (x-1, y), p2 (x, y-1);
	    GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
	  }
    }
    break;

  case Horizontal:
    {
      const DimSideImg y1 = rect.y-1, y2 = rect.y;
      for (DimSideImg x = firstX; x < endX; ++x) {
	Point p1 (x, y1), p2 (x, y2);
	GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
      }
      if (connectivity & Connectivity::C6P)
	for (DimSideImg x = firstX+1; x < endX; ++x) {
	  Point p1 (x-1, y1), p2 (x, y2);
	  GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
	}
      if (connectivity & Connectivity::C6N)
	for (DimSideImg x = firstX+1; x < endX; ++x) {
	  Point p1 (x, y1), p2 (x-1, y2);
	  GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
	}
    }
    break;

  case Vertical:
    {
      const DimSideImg x1 = rect.x-1, x2 = rect.x;
      for (DimSideImg y = firstY; y < endY; ++y) {
	Point p1 (x1, y), p2 (x2, y);
	GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
      }
      if (connectivity & Connectivity::C6P)
	for (DimSideImg y = firstY+1; y < endY; ++y) {
	  Point p1 (x1, y-1), p2 (x2, y);
	  GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
	}
      if (connectivity & Connectivity::C6N)
	for (DimSideImg y = firstY+1; y < endY; ++y)  {
	  Point p1 (x1, y), p2 (x2, y-1);
	  GRAPHWALKER_CALL_LAMBDA_PTS (p1, p2, edgeCount, lambda (p1, p2));
	}
    }
    break;

  default:
    BOOST_ASSERT (false);
  }
  return edgeCount;
}

// ========================================
template <typename WeightT, typename EdgeWeightFunction>
inline WeightT
GraphWalker::getMedian (const EdgeWeightFunction &ef /*ef.getValue (idx)*/) const {
  int bits = 8*sizeof (WeightT);
  BOOST_ASSERT (bits < 17);
  DimEdge dim = 1UL << bits;
  std::vector<DimImg> indices (dim+1, 0);
  DimImg *histogram = &indices[1];
  forEachVertexIdx ([&histogram, &ef] (const DimImg &idx) {
      ++histogram[(size_t) ef.getValue (idx)];
    });
  std::partial_sum (histogram, histogram+dim, histogram);
  return lower_bound (indices.begin (), indices.end (), indices[dim] >> 1) - indices.begin ();
}

// ========================================
template <typename WeightT, typename EdgeWeightFunction>
inline DimEdge
GraphWalker::getEdges (const Rect &rect, TileItem tileItem, Edge<WeightT> *edges, const EdgeWeightFunction &ef /*ef.getWeight (p0, p1)*/) const {
  BOOST_ASSERT (rect.width);
  BOOST_ASSERT (rect.height);
  BOOST_ASSERT (edges);
  DimEdge edgeCount = 0;
  forEachEdgePt (rect, tileItem,
		 [&edges, &edgeCount, &ef] (Point const &a, Point const &b) {
		   Edge<WeightT> &edge = edges[edgeCount++];
		   edge.points[0] = a;
		   edge.points[1] = b;
		   edge.weight = ef.getWeight (a, b);
		 });
  return edgeCount;
}

// ----------------------------------------
template <typename WeightT, typename EdgeWeightFunction>
inline DimEdge
GraphWalker::getSortedEdges (const Rect &rect, TileItem tileItem, Edge<WeightT> *edges, const EdgeWeightFunction &ef /*ef.getWeight (p0, p1) ef.sort ()*/) const {
  DimEdge edgeCount = getEdges (rect, tileItem, edges, ef);
  ef.sort (edges, edgeCount);
  return edgeCount;
}

template <typename WeightT, typename EdgeWeightFunction>
// template DimImg
inline DimEdge
GraphWalker::getCountingSortedEdges (const Rect &rect, TileItem tileItem, Edge<WeightT> *edges, const EdgeWeightFunction &ef /*ef.getWeight (p0, p1)*/) const {
  int bits = 8*sizeof (WeightT);
  BOOST_ASSERT (bits < 17);
  DimEdge dim = 1UL << bits;
  DimImg *indices = new DimImg [dim+1];
  DimImg *histogram = indices+1;
  DimImg sum = 0;
  ::std::fill_n (histogram, dim, 0);
  forEachEdgePt (rect, tileItem,
		 [&histogram, &ef] (Point const &a, Point const &b) {
		   ++histogram[(size_t) ef.getWeight (a, b)];
		 });
  
  // get indices by prefix sum
  std::partial_sum (histogram, histogram+dim, histogram);
  sum = indices[dim];
  if (ef.getDecr ()) {
    for (size_t i = 0; i < dim; ++i)
      histogram[i] = sum - histogram[i];
    indices++;
  } else
    indices[0] = 0;
  // extract edges
  forEachEdgePt (rect, tileItem,
		 [
#ifdef SMART_LOG
		  this,
#endif
		  &ef, &edges, &indices] (Point const &a, Point const &b) {
		   WeightT weight = ef.getWeight (a, b);
		   Edge<WeightT> &edge = edges[indices[(size_t) weight]++];
		   edge.points[0] = a;
		   edge.points[1] = b;
		   edge.weight = weight;
#ifdef SMART_LOG
		   printEdge (cerr, edge, size) << endl;
#endif
		 });
#ifdef SMART_LOG
  cerr << endl;
  for (int i = 0; i < sum; ++i)
    printEdge (cerr, edges[i], size) << endl;
#endif
  delete [] (histogram-1);
  return sum;
}

#endif // _OTB_TRISKELE_ARRAY_TREE_GRAPHWALKER_TPP