triskele/include/ArrayTree/GraphWalker.tpp

#ifndef _OTB_TRISKELE_ARRAY_TREE_GRAPHWALKER_TPP
#define _OTB_TRISKELE_ARRAY_TREE_GRAPHWALKER_TPP

// ========================================
inline const Size &
GraphWalker::getSize () const {
  return border.getSize ();
}

// ========================================
inline
GraphWalker::GraphWalker (const Border &border, const Connectivity &connectivity)
  : border (border), connectivity (connectivity) {
  DEF_LOG ("GraphWalker::GraphWalker",  "size: " << border.getSize ());
  BOOST_ASSERT (border.getSize ().width  >= DimImg (0));
  BOOST_ASSERT (border.getSize ().height >= DimImg (0));
}

// ----------------------------------------
inline DimImg
GraphWalker::vertexMaxCount () const {
  return vertexMaxCount (border.getSize ());
}

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

// ----------------------------------------
inline DimEdge
GraphWalker::edgeMaxCount () const {
  return edgeMaxCount (border.getSize ());
}

// ----------------------------------------
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 &rect,
		       vector<Rect> &tiles, vector<Rect> &boundaries, vector<bool> &verticalBoundaries) {
  BOOST_ASSERT (coreCount);
  DEF_LOG ("GraphWalker::setTiles", "coreCount:" << coreCount << " rect:" << rect);
  if (coreCount < 2 || max (rect.width, rect.height) < 4 || min (rect.width, rect.height) < 3) {
    LOG ("rect:" << rect);
    tiles.push_back (rect);
    return;
  }
  bool vertical = rect.width > rect.height;
  bool odd = coreCount & 0x1;
  DimImg thin = (vertical ? rect.width : rect.height) / (odd ? coreCount : 2);
  Rect tileA (rect);
  Rect tileB (rect);
  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);
}

inline void
GraphWalker::setMaxBounds (const vector<Rect> &tiles, vector<DimImg> &vertexMaxBounds, vector<DimImg> &edgesMaxBounds) const {
  unsigned int tileCount = tiles.size ();
  vertexMaxBounds.resize (tileCount+1);
  edgesMaxBounds.resize (tileCount+1);
  DimImg vertexBase = 0, edgesBase = 0;
  vertexMaxBounds [0] = edgesMaxBounds [0] = 0;
  for (unsigned int i = 0; i < tileCount; ++i) {
    Size zoneSize (tiles[i].width, tiles[i].height);
    vertexMaxBounds [i+1] = vertexBase += vertexMaxCount (zoneSize); /* -1, but prety LOG */
    edgesMaxBounds [i+1] = edgesBase += edgeMaxCount (zoneSize);
  }
}

// ========================================
template<typename Funct>
inline DimImg
GraphWalker::forEachVertexIdx (const Funct &lambda /*lambda (DimImg idx)*/) const {
  DimImg maxCount = vertexMaxCount ();
  DimImg vertexCount = 0;
  for (DimImg 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 (DimImg idx)*/) const {
  const Size &size (border.getSize ());
  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 (Point p)*/) const {
  const Size &size (border.getSize ());
  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 (Point p0, Point 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;
  vector<DimImg> indices (dim+1, 0);
  DimImg *histogram = &indices[1];
  forEachVertexIdx ([&histogram, &ef] (const DimImg &idx) {
      ++histogram[(size_t) ef.getValue (idx)];
    });
  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 {
  SMART_DEF_LOG ("GraphWalker::getSortedEdges", "");
  DimEdge edgeCount = getEdges (rect, tileItem, edges, ef);
  // SMART_LOG (endl << printEdges (edges, border.getSize (), edgeCount));
  ef.sort (edges, edgeCount);
  SMART_LOG (endl << printEdges (edges, border.getSize (), 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 {
  SMART_DEF_LOG ("GraphWalker::getCountingSortedEdges", "");
  int bits = 8*sizeof (WeightT);
  BOOST_ASSERT (bits < 17);
  DimEdge dim = 1UL << bits;
  vector<DimImg> counters (dim+1, 0);
  DimImg *indices = &counters [0];
  DimImg *histogram = &indices [1];
  DimImg sum = 0;
  forEachEdgePt (rect, tileItem,
		 [&histogram, &ef] (Point const &a, Point const &b) {
		   ++histogram[(size_t) ef.getWeight (a, b)];
		 });
  
  // get indices by prefix sum
  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
#ifdef ENABLE_SMART_LOG
  const Size &size (border.getSize ());
#endif
  forEachEdgePt (rect, tileItem,
		 [
#ifdef ENABLE_SMART_LOG
		  &size,
#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;
		   SMART_LOG_EXPR (cerr << printEdge (edge, size) << endl);
		 });
  SMART_LOG (endl << printEdges (edges, size, sum));
  return sum;
}

#endif // _OTB_TRISKELE_ARRAY_TREE_GRAPHWALKER_TPP