lsgrmGraphOperations.txx 31.7 KB
Newer Older
1 2
#ifndef __LSGRM_GRAPH_OPERATIONS_TXX
#define __LSGRM_GRAPH_OPERATIONS_TXX
remicres's avatar
remicres committed
3
//#include "lsgrmGraphOperations.h"
4 5
//#include <unistd.h>
#include <cstdio>
6 7 8 9

namespace lsgrm
{

remicres's avatar
remicres committed
10 11 12 13 14
template<class TSegmenter>
typename TSegmenter::ImageType::Pointer ReadImageRegion(
    typename TSegmenter::ImageType * inputPtr,
    typename TSegmenter::ImageType::RegionType region)
{
15

remicres's avatar
remicres committed
16 17 18 19 20 21 22
  typedef typename otb::MultiChannelExtractROI<typename TSegmenter::ImageType::InternalPixelType,
      typename TSegmenter::ImageType::InternalPixelType> ExtractROIFilterType;
  typename ExtractROIFilterType::Pointer filter = ExtractROIFilterType::New();
  filter->SetInput(inputPtr);
  filter->SetExtractionRegion(region);
  filter->SetReleaseDataFlag(true);
  filter->Update();
23

remicres's avatar
remicres committed
24 25
  return filter->GetOutput();
}
26

remicres's avatar
remicres committed
27
template<class TSegmenter>
28
typename TSegmenter::GraphType
remicres's avatar
remicres committed
29
MergeAllGraphsAndAchieveSegmentation(
30
    const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
31 32 33 34 35 36
    const float& threshold,
    std::vector<ProcessingTile>& tiles,
    const unsigned int nbTilesX,
    const unsigned int nbTilesY,
    const unsigned int imageWidth,
    const unsigned int imageHeight,
37 38
    const unsigned int imageBands,
    unsigned int numberOfIterations)
remicres's avatar
remicres committed
39
{
40

41 42
  std::cout << "--- Graph aggregation...\n" << std::endl;

remicres's avatar
remicres committed
43
  TSegmenter segmenter;
44

remicres's avatar
remicres committed
45 46
  std::cout << "Reading graphs" << std::endl;

remicres's avatar
remicres committed
47 48 49 50
  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; col++)
      {
remicres's avatar
remicres committed
51 52
      std::cout << "\tImporting graph of tile " << (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) << std::endl;

53
      InsertNodesFromTile<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], false);
remicres's avatar
remicres committed
54 55
      }
    }
56

remicres's avatar
remicres committed
57
  std::cout << "Removing duplicated nodes and updating neighbors" << std::endl;
58

remicres's avatar
remicres committed
59 60 61 62 63
  // TODO this might be parallelized...
  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; col++)
      {
remicres's avatar
remicres committed
64
      std::cout << "Cleaning nodes of tile " << (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) << std::endl;
65

66
      std::unordered_map<std::size_t,
remicres's avatar
remicres committed
67
      std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
68

remicres's avatar
remicres committed
69
      std::cout << "\tBuildBorderPixelMap..." << std::endl;
70

remicres's avatar
remicres committed
71
      BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], row, col,
72
          nbTilesX, nbTilesY, borderPixelMap, imageWidth, true);
73

remicres's avatar
remicres committed
74
      std::cout << "\tRemoveDuplicatedNodes..." << std::endl;
75

remicres's avatar
remicres committed
76
      RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
77

remicres's avatar
remicres committed
78
      std::cout << "\tUpdateNeighborsOfNoneDuplicatedNodes..." << std::endl;
79

remicres's avatar
remicres committed
80
      UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap, imageWidth, imageHeight);
remicres's avatar
remicres committed
81 82 83
      }
    }

remicres's avatar
remicres committed
84
  std::cout << "Achieve segmentation process" ;
remicres's avatar
remicres committed
85 86 87 88 89 90 91

  // Segmentation of the graph
  segmenter.SetImageWidth(imageWidth);
  segmenter.SetImageHeight(imageHeight);
  segmenter.SetNumberOfComponentsPerPixel(imageBands);
  segmenter.SetParam(params);
  segmenter.SetThreshold(threshold);
92
  segmenter.SetDoFastSegmentation(false); // was true
93
  segmenter.SetNumberOfIterations(numberOfIterations);
94
  grm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
remicres's avatar
remicres committed
95

96
  // patch raf (sort using node id directly --> seems ok)
97
  // Sort the nodes top-->down and left-->right
98
  /*
99 100 101 102 103 104
  std::sort(segmenter.m_Graph.m_Nodes.begin(), segmenter.m_Graph.m_Nodes.end(),
      [imageWidth](const auto & a, const auto & b) -> bool
  {
    lp::CellLists borderPixelsA, borderPixelsB;
    lp::ContourOperations::GenerateBorderCells(borderPixelsA, a->m_Contour, a->m_Id, imageWidth);
    lp::ContourOperations::GenerateBorderCells(borderPixelsB, b->m_Contour, b->m_Id, imageWidth);
105

106 107 108 109
	// patch raf (were unsigned int)
    std::size_t offA = 0;
    std::size_t offB = 0;
	// end patch
110 111 112
    for (auto& pix: borderPixelsA)
      if (pix>offA)
        offA=pix;
113

114 115 116
    for (auto& pix: borderPixelsB)
      if (pix>offB)
        offB=pix;
117

118 119
      return offA > offB;
  });
120 121 122 123 124 125 126
  */
  std::sort(segmenter.m_Graph.m_Nodes.begin(), segmenter.m_Graph.m_Nodes.end(),
      [](const auto & a, const auto & b) -> bool
  {
      return a->m_Id < b->m_Id;
  });
  // end patch
127

128
  return segmenter.m_Graph;
129
}
130

remicres's avatar
remicres committed
131
template<class TSegmenter>
132
long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
    const float& threshold,
    const unsigned int niter,
    std::vector<ProcessingTile>& tiles,
    const unsigned int nbTilesX,
    const unsigned int nbTilesY,
    const unsigned int imageWidth,
    const unsigned int imageHeight,
    const unsigned int imageBands,
    bool& isFusion)
{
  long long unsigned int accumulatedMemory = 0;
  isFusion = false;

  const unsigned int numberOfNeighborLayers = static_cast<unsigned int>(pow(2, niter + 1) - 2);

  std::cout << "--- Running partial segmentations...\nNumber of neighbor layers " << numberOfNeighborLayers << std::endl;

  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; col++)
      {
154
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
155
      if (MyTurn(row*nbTilesX + col))
156
#endif
remicres's avatar
remicres committed
157
        {
158
        // Get the current tile
159
        std::cout << "Processing tile " << row << ", " << col << std::endl;
160
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
161

remicres's avatar
remicres committed
162
        // Load the graph
163 164
        std::cout << "\tLoad graph..." << std::endl;
        TSegmenter segmenter;
165
        ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
remicres's avatar
remicres committed
166 167 168 169

        // Add stability margin to the graph
        {
          std::cout << "\tAdd stability margin..." << std::endl;
170
          AddStabilityMargin<TSegmenter>(segmenter.m_Graph, tiles,
remicres's avatar
remicres committed
171
              row, col, nbTilesX, nbTilesY);
remicres's avatar
remicres committed
172 173

          std::cout << "\tBuild border pixel map..." << std::endl;
174
          std::unordered_map<std::size_t, std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
175
          BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, currentTile, row, col,
remicres's avatar
remicres committed
176 177 178 179 180 181 182 183 184 185 186
              nbTilesX, nbTilesY, borderPixelMap, imageWidth);

          std::cout << "\tRemove duplicated nodes..." << std::endl;
          RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);

          std::cout << "\tUpdate neighbors.." << std::endl;
          UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap,
              imageWidth,
              imageHeight);

          std::cout << "\tRemove useless.." << std::endl;
187
          RemoveUselessNodes<TSegmenter>(currentTile, segmenter.m_Graph,
188
              imageWidth, numberOfNeighborLayers);
remicres's avatar
remicres committed
189 190 191 192 193 194 195 196 197

        }

        // Segmentation of the graph
        segmenter.SetImageWidth(imageWidth);
        segmenter.SetImageHeight(imageHeight);
        segmenter.SetNumberOfComponentsPerPixel(imageBands);
        segmenter.SetParam(params);
        segmenter.SetThreshold(threshold);
198
        segmenter.SetDoFastSegmentation(false);
remicres's avatar
remicres committed
199 200
        segmenter.SetNumberOfIterations(niter);

201
        std::cout << "\tPartial segmentation";
202
        auto merge = grm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
remicres's avatar
remicres committed
203 204 205 206 207

        if(merge == true)
          isFusion = true;

        // Remove unstable segments
208
        std::cout << "\tRemove unstable segments..." << std::endl;
209
        RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
remicres's avatar
remicres committed
210 211

        // Retrieve the amount of memory to store this graph
212
        std::cout << "\tGet graph memory..." << std::endl;
213
        accumulatedMemory += segmenter.GetGraphMemory();
remicres's avatar
remicres committed
214

215
        // Write graph to temporay directory
216
        std::cout << "\tWrite graph..." << std::endl;
217
        WriteGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
remicres's avatar
remicres committed
218 219 220 221 222 223
        }
      }
    }

#ifdef OTB_USE_MPI
  otb::MPIConfig::Instance()->barrier();
224 225
#endif

remicres's avatar
remicres committed
226 227 228 229 230 231 232
  std::cout << "Add stability margins to graph for the next round..."<< std::endl;

  // During this step we extract the stability margin for the next round
  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; col++)
      {
233
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
234
      if (MyTurn(row*nbTilesX + col))
235
#endif
remicres's avatar
remicres committed
236
        {
237 238
        // Get current tile
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
remicres's avatar
remicres committed
239 240

        // Load the graph
241 242
        typename TSegmenter::GraphType graph;
        ReadGraph<TSegmenter>(graph, currentTile.nodeFileName, currentTile.edgeFileName);
remicres's avatar
remicres committed
243

244
        // Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight-1
remicres's avatar
remicres committed
245 246 247 248
        // and write them to the stability margin
        {
          std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;

249
          DetectBorderNodes<TSegmenter>(graph, currentTile,
remicres's avatar
remicres committed
250 251 252 253
              borderNodeMap, imageWidth, imageHeight);

          ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);

254
          WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
remicres's avatar
remicres committed
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
        }
        std::cout << "Process finished on tile " << (row*nbTilesX + col) << std::endl;
        }
      }
    }
  std::cout << std::endl;

  return accumulatedMemory;
}

template<class TSegmenter>
void RemoveUselessNodes(ProcessingTile& tile,
    typename TSegmenter::GraphType& graph,
    const unsigned int imageWidth,
    const unsigned int numberOfLayers)
{
  using NPtr = typename TSegmenter::NodePointerType;
  using UI = unsigned int;
  unsigned int rowPixel, colPixel;
  std::unordered_map<NPtr, UI> marginNodes;

  for(auto& node : graph.m_Nodes)
    {
    if(node->m_Bbox.m_UX > tile.columns[0] &&
        node->m_Bbox.m_UY > tile.rows[0] &&
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
      continue;
    else if(node->m_Bbox.m_UX > tile.columns[1] ||
        node->m_Bbox.m_UY > tile.rows[1] ||
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[0] ||
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[0])
      continue;
    else
      {
      lp::CellLists borderpixels;
      lp::ContourOperations::GenerateBorderCells(borderpixels, node->m_Contour, node->m_Id, imageWidth);

      for(auto& pix : borderpixels)
        {
        rowPixel = pix / imageWidth;
        colPixel = pix % imageWidth;

        if(rowPixel == tile.rows[0] || rowPixel == tile.rows[1])
          {
          if(colPixel >= tile.columns[0] && colPixel <= tile.columns[1])
            {
            marginNodes.insert(std::pair<NPtr, UI>(node, 0));
            break;
            }
          }
        else if(colPixel == tile.columns[0] || colPixel == tile.columns[1])
          {
          if(rowPixel >= tile.rows[0] && rowPixel <= tile.rows[1])
            {
            marginNodes.insert(std::pair<NPtr, UI>(node, 0));
            break;
            }
          }
        else
          continue;
        }
      }
    }

  ExtractStabilityMargin<TSegmenter>(marginNodes, numberOfLayers);

  for(auto& node : graph.m_Nodes)
    {
    if(node->m_Bbox.m_UX > tile.columns[0] &&
        node->m_Bbox.m_UY > tile.rows[0] &&
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
      continue;
    else if(marginNodes.find(node) == marginNodes.end())
      {
      RemoveEdgeToUnstableNode<TSegmenter>(node);
      node->m_Expired = true;
      }
    else
      continue;
    }

338
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
339 340 341
}

template<class TSegmenter>
342
void UpdateNeighborsOfNoneDuplicatedNodes(std::unordered_map<std::size_t,
remicres's avatar
remicres committed
343 344 345 346 347 348 349 350 351
    std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
    const unsigned int imageWidth,
    const unsigned int imageHeight)
{
  using EdgeType = typename TSegmenter::EdgeType;
  unsigned int boundary;

  for(auto& pn : borderPixelMap)
    {
352
    NeighIDType neighborhood[4];
353
    grm::FOURNeighborhood(neighborhood, pn.first, imageWidth, imageHeight);
remicres's avatar
remicres committed
354 355 356 357 358 359 360 361 362 363 364 365 366

    for(short j = 0; j < 4; ++j)
      {
      if(neighborhood[j] > -1)
        {
        auto isNeigh = borderPixelMap.find(neighborhood[j]);
        if(isNeigh != borderPixelMap.end())
          {
          auto currNode = pn.second.front();
          auto neigh = isNeigh->second.front();

          if(currNode != neigh)
            {
367
            auto toNeigh = grm::GraphOperations<TSegmenter>::FindEdge(currNode, neigh);
remicres's avatar
remicres committed
368 369 370 371 372 373 374 375 376 377 378 379

            if(toNeigh == currNode->m_Edges.end())
              {
              boundary = 0;

              lp::CellLists borderPixels;
              lp::ContourOperations::GenerateBorderCells(borderPixels, currNode->m_Contour, currNode->m_Id, imageWidth);

              for(auto pix : borderPixels)
                {
                if(borderPixelMap.find(pix) != borderPixelMap.end())
                  {
380
                  NeighIDType pixNeighborhood[4];
381
                  grm::FOURNeighborhood(pixNeighborhood, pix, imageWidth, imageHeight);
remicres's avatar
remicres committed
382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406

                  for(short k = 0; k < 4; k++)
                    {
                    if(pixNeighborhood[k] > -1)
                      {
                      auto isNeighPix = borderPixelMap.find(pixNeighborhood[k]);

                      if(isNeighPix != borderPixelMap.end() && isNeighPix->second.front() == neigh)
                        boundary++;
                      }
                    }
                  }
                }

              currNode->m_Edges.push_back(EdgeType(neigh, 0, boundary));
              neigh->m_Edges.push_back(EdgeType(currNode, 0, boundary));
              }
            }
          }
        }
      }
    }
}

template<class TSegmenter>
407
void RemoveDuplicatedNodes(std::unordered_map<std::size_t,
remicres's avatar
remicres committed
408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430
    std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
    typename TSegmenter::GraphType& graph,
    const unsigned int imageWidth)
{
  using EdgeType = typename TSegmenter::EdgeType;
  unsigned int boundary;

  // Explore the border nodes
  for(auto& pn : borderPixelMap)
    {
    // no valid nodes are those who have already been considered.
    if(pn.second.size() > 1)
      {
      auto refNode = pn.second.front();

      // Explore duplicated nodes
      for(auto nit = pn.second.begin() + 1; nit != pn.second.end(); ++nit)
        {
        // Explore their edges
        for(auto& edg : (*nit)->m_Edges)
          {
          auto neighNit = edg.GetRegion();

431
          auto toNit = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, *nit);
remicres's avatar
remicres committed
432 433 434 435 436 437 438

          assert(toNit != edg.GetRegion()->m_Edges.end());

          boundary = edg.m_Boundary;

          neighNit->m_Edges.erase(toNit);

439
          auto toRefNode = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, refNode);
remicres's avatar
remicres committed
440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466

          if(toRefNode == neighNit->m_Edges.end())
            {
            // Create an edge neighNit -> refNode
            neighNit->m_Edges.push_back(EdgeType(refNode, 0, boundary));
            // Create an edge refNode -> neighNit
            refNode->m_Edges.push_back(EdgeType(neighNit, 0, boundary));
            }
          }

        (*nit)->m_Expired = true;
        }

      lp::CellLists borderPixels;
      lp::ContourOperations::GenerateBorderCells(borderPixels, refNode->m_Contour, refNode->m_Id, imageWidth);

      for(auto& pix : borderPixels)
        {
        if(borderPixelMap.find(pix) != borderPixelMap.end())
          {
          borderPixelMap[pix].clear();
          borderPixelMap[pix].push_back(refNode);
          }
        }
      }
    }

467
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
468 469 470 471 472 473 474 475 476 477
}

// Build the map assigning for each border pixel on the borders the node it belongs to.
template<class TSegmenter>
void BuildBorderPixelMap(typename TSegmenter::GraphType& graph,
    ProcessingTile& tile,
    const unsigned int rowTile,
    const unsigned int colTile,
    const unsigned int nbTilesX,
    const unsigned int nbTilesY,
478
    std::unordered_map<std::size_t,
remicres's avatar
remicres committed
479
    std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
480 481
    const unsigned int imageWidth,
	bool merging)
remicres's avatar
remicres committed
482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
{
  unsigned int rowPixel, colPixel;
  unsigned int rowMin = (tile.rows[0] > 0) ? tile.rows[0] - 1 : tile.rows[0];
  unsigned int rowMax = tile.rows[1] + 1;
  unsigned int colMin = (tile.columns[0] > 0) ? tile.columns[0] - 1 : tile.columns[0];
  unsigned int colMax = tile.columns[1] + 1;

  for(auto& node : graph.m_Nodes)
    {
    if(node->m_Bbox.m_UX > tile.columns[0] &&
        node->m_Bbox.m_UY > tile.rows[0] &&
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
      {
      continue;
      }
498 499 500 501 502 503 504 505 506
	// patch raf (exclude nodes all outside tiles - on merging stage only)
	else if(merging && (node->m_Bbox.m_UX > tile.columns[1] ||
        node->m_Bbox.m_UY > tile.rows[1] ||
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[0] ||
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[0]))
      {
      continue;
      }  
	// end patch
remicres's avatar
remicres committed
507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539
    else
      {
      lp::CellLists borderPixels;
      lp::ContourOperations::GenerateBorderCells(borderPixels, node->m_Contour, node->m_Id, imageWidth);

      for(auto& pix : borderPixels)
        {
        rowPixel = pix / imageWidth;
        colPixel = pix % imageWidth;

        if(rowTile > 0 && (rowPixel == tile.rows[0] || rowPixel == rowMin))
          {
          borderPixelMap[pix].push_back(node);
          }
        else if(colTile < nbTilesX - 1 && ( colPixel == tile.columns[1] || colPixel == colMax))
          {
          borderPixelMap[pix].push_back(node);
          }
        else if(rowTile < nbTilesY - 1 && (rowPixel == tile.rows[1] || rowPixel == rowMax))
          {
          borderPixelMap[pix].push_back(node);
          }
        else if(colTile > 0 && ( colPixel == tile.columns[0] || colPixel == colMin))
          {
          borderPixelMap[pix].push_back(node);
          }
        else
          continue;
        }
      }
    }
}

540 541
template<class TSegmenter>
void InsertNodesFromTile(typename TSegmenter::GraphType& graph,
542
    ProcessingTile& tile, bool margin)
543 544 545 546 547 548 549 550 551 552 553 554 555 556 557
{
  typename TSegmenter::GraphType subgraph;
  if (margin)
    {
    ReadGraph<TSegmenter>(subgraph, tile.nodeMarginFileName, tile.edgeMarginFileName);
    }
  else
    {
    ReadGraph<TSegmenter>(subgraph, tile.nodeFileName, tile.edgeFileName);
    }

  graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());

}

remicres's avatar
remicres committed
558 559
template<class TSegmenter>
void AddStabilityMargin(typename TSegmenter::GraphType& graph,
560
    std::vector<ProcessingTile>& tiles,
remicres's avatar
remicres committed
561 562 563
    const unsigned int row,
    const unsigned int col,
    const unsigned int nbTilesX,
remicres's avatar
remicres committed
564
    const unsigned int nbTilesY)
remicres's avatar
remicres committed
565 566 567 568 569
{

  // Margin to retrieve at top
  if(row > 0)
    {
570
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + col]);
remicres's avatar
remicres committed
571 572 573 574 575
    }

  // Margin to retrieve at right
  if(col < nbTilesX - 1)
    {
576
    InsertNodesFromTile<TSegmenter>(graph, tiles[row * nbTilesX + (col+1)]);
remicres's avatar
remicres committed
577 578 579 580 581
    }

  // Margin to retrieve at bottom
  if(row < nbTilesY - 1)
    {
582
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + col]);
remicres's avatar
remicres committed
583
    }
584

remicres's avatar
remicres committed
585 586 587
  // Margin to retrieve at left
  if(col > 0)
    {
588
    InsertNodesFromTile<TSegmenter>(graph, tiles[row * nbTilesX + (col-1)]);
remicres's avatar
remicres committed
589 590 591 592 593
    }

  // Margin to retrieve at top right
  if(row > 0 && col < nbTilesX - 1)
    {
594
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + (col+1)]);
remicres's avatar
remicres committed
595 596 597 598 599
    }

  // Margin to retrieve at bottom right
  if(row < nbTilesY - 1 && col < nbTilesX - 1)
    {
600
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + (col+1)]);
remicres's avatar
remicres committed
601 602 603 604 605
    }

  // Margin to retrieve at bottom left
  if(row < nbTilesY - 1 && col > 0)
    {
606
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + (col-1)]);
remicres's avatar
remicres committed
607 608 609 610 611
    }

  // Margin to retrieve at top left
  if(row > 0 && col > 0)
    {
612
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + (col-1)]);
remicres's avatar
remicres committed
613 614 615 616 617 618
    }
}

template<class TSegmenter>
long long unsigned int RunFirstPartialSegmentation(
    typename TSegmenter::ImageType * inputPtr,
619
    const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
620 621 622 623 624 625 626 627 628 629 630 631
    const float& threshold,
    const unsigned int niter,
    const unsigned int niter2,
    std::vector<ProcessingTile>& tiles,
    const unsigned int nbTilesX,
    const unsigned int nbTilesY,
    const unsigned int tileWidth,
    const unsigned int tileHeight,
    bool& isFusion)
{

  using ImageType = typename TSegmenter::ImageType;
632 633
  const unsigned int imageWidth = inputPtr->GetLargestPossibleRegion().GetSize()[0];
  const unsigned int imageHeight = inputPtr->GetLargestPossibleRegion().GetSize()[1];
remicres's avatar
remicres committed
634 635 636 637 638 639 640 641

  long long unsigned int accumulatedMemory = 0;
  isFusion = false;

  const unsigned int numberOfNeighborLayers = static_cast<unsigned int>(pow(2, niter2 + 1) - 2);

  std::cout << "--- Running fist partial segmentation...\nNumber of neighbor layers " << numberOfNeighborLayers << std::endl;

Gaetano Raffaele's avatar
Gaetano Raffaele committed
642
  for(unsigned int row = 0; row < nbTilesY; ++row)
remicres's avatar
remicres committed
643
    {
Gaetano Raffaele's avatar
Gaetano Raffaele committed
644
    for(unsigned int col = 0; col < nbTilesX; col++)
remicres's avatar
remicres committed
645
      {
646
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
647
      if (MyTurn(row*nbTilesX + col))
648
#endif
remicres's avatar
remicres committed
649 650 651
        {
        // Reading images
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
652

653
        std::cout << "Processing tile " <<  (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) <<
remicres's avatar
remicres committed
654 655 656 657 658 659 660
            " (" << col << ", " << row << ")" <<
            " start: [" << currentTile.region.GetIndex()[0] << ", " << currentTile.region.GetIndex()[1] <<
            "] size: [" << currentTile.region.GetSize()[0] << ", " << currentTile.region.GetSize()[1] << "]" << std::endl;

        typename ImageType::Pointer imageTile = ReadImageRegion<TSegmenter>(inputPtr, currentTile.region);

        // Segmenting image
remicres's avatar
remicres committed
661
        std::cout << "\tSegmenting";
remicres's avatar
remicres committed
662 663 664
        TSegmenter segmenter;
        segmenter.SetParam(params);
        segmenter.SetThreshold(threshold);
665
        segmenter.SetDoFastSegmentation(false);
remicres's avatar
remicres committed
666 667 668
        segmenter.SetNumberOfIterations(niter);
        segmenter.SetInput(imageTile);
        segmenter.Update();
669
		
remicres's avatar
remicres committed
670 671 672 673
        if(segmenter.GetComplete() == false)
          isFusion = true;

        // Rescale the graph to be in the reference of the image
674
        std::cout << "\tRescaling graph..." << std::endl;
remicres's avatar
remicres committed
675
        RescaleGraph<TSegmenter>(segmenter.m_Graph,
676
            currentTile,
remicres's avatar
remicres committed
677 678 679 680 681
            row,
            col,
            tileWidth,
            tileHeight,
            imageWidth);
682
			
remicres's avatar
remicres committed
683
        // Remove unstable segments
684
        std::cout << "\tRemoving unstable segments..." << std::endl;
685
        RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
remicres's avatar
remicres committed
686 687

        // Retrieve the amount of memory to store this graph
688
        std::cout << "\tRetrieving graph memory..." << std::endl;
689
        accumulatedMemory += segmenter.GetGraphMemory();
remicres's avatar
remicres committed
690

691
        // Write graph to temporay directory
692
        std::cout << "\tWriting graph..." << std::endl;
693
        WriteGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
remicres's avatar
remicres committed
694

695
        // Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight -1
remicres's avatar
remicres committed
696
        // and write them to the stability margin
697
        std::cout << "\tComputing stability margin..." << std::endl;
remicres's avatar
remicres committed
698
        {
remicres's avatar
remicres committed
699
        std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;
remicres's avatar
remicres committed
700

remicres's avatar
remicres committed
701 702
        DetectBorderNodes<TSegmenter>(segmenter.m_Graph, currentTile,
            borderNodeMap, imageWidth, imageHeight);
remicres's avatar
remicres committed
703

remicres's avatar
remicres committed
704
        ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
remicres's avatar
remicres committed
705

remicres's avatar
remicres committed
706
        WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
remicres's avatar
remicres committed
707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817
        }
        }
      } // for each col
    } // for each row

  return accumulatedMemory;
}

template<class TSegmenter>
void ExtractStabilityMargin(std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& nodeMap,
    const unsigned int pmax)
{
  std::vector<typename TSegmenter::NodePointerType> startingNodes;
  startingNodes.reserve(nodeMap.size());
  for(auto& kv: nodeMap)
    startingNodes.push_back(kv.first);

  for(auto& n : startingNodes)
    {
    ExploreDFS<TSegmenter>(n, 0, nodeMap, pmax);
    }
}

template<class TSegmenter>
void ExploreDFS(typename TSegmenter::NodePointerType s,
    const unsigned int p,
    std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& Cb,
    const unsigned int pmax)
{
  if(p > pmax)
    return;
  else
    {
    if(Cb.find(s) != Cb.end())
      {
      if(p <= Cb[s])
        {
        Cb[s] = p;
        for(auto edg : s->m_Edges)
          {
          ExploreDFS<TSegmenter>(edg.GetRegion(), p + 1, Cb, pmax);
          }
        }
      else
        return;
      }
    else
      {
      Cb[s] = p;
      for(auto edg : s->m_Edges)
        {
        ExploreDFS<TSegmenter>(edg.GetRegion(), p + 1, Cb, pmax);
        }
      }
    }
}

template<class TSegmenter>
void DetectBorderNodes(typename TSegmenter::GraphType& graph,
    ProcessingTile& tile,
    std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& borderNodeMap,
    const unsigned int imageWidth,
    const unsigned int imageHeight)
{
  using NP = typename TSegmenter::NodePointerType;
  using Uint = unsigned int;
  unsigned int rowPixel, colPixel;

  for(auto& node : graph.m_Nodes)
    {
    if(node->m_Bbox.m_UX > tile.columns[0] &&
        node->m_Bbox.m_UY > tile.rows[0] &&
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
      continue;
    else
      {
      lp::CellLists borderpixels;
      lp::ContourOperations::GenerateBorderCells(borderpixels, node->m_Contour, node->m_Id, imageWidth);

      for(auto& pix : borderpixels)
        {
        rowPixel = pix / imageWidth;
        colPixel = pix % imageWidth;

        if(tile.rows[0] > 0 && rowPixel == tile.rows[0])
          {
          borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
          break;
          }
        else if(tile.columns[1] < imageWidth - 1 && colPixel == tile.columns[1])
          {
          borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
          break;
          }
        else if(tile.rows[1] < imageHeight - 1 && rowPixel == tile.rows[1])
          {
          borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
          break;
          }
        else if(tile.columns[0] > 0 && colPixel == tile.columns[0])
          {
          borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
          break;
          }
        else continue;
        }
      }
    }
}

818 819

// Generic
remicres's avatar
remicres committed
820
template<class TSegmenter>
821
void ReadGraph(TSegmenter& segmenter,
remicres's avatar
remicres committed
822 823 824
    const std::string& nodesPath,
    const std::string& edgesPath)
{
825
  FILE * nodeStream = fopen(nodesPath.c_str(), "rb");
remicres's avatar
remicres committed
826
  assert(nodeStream != NULL);
827
  FILE * edgeStream = fopen(edgesPath.c_str(), "rb");
remicres's avatar
remicres committed
828 829
  assert(edgeStream != NULL);

830
  segmenter.ReadGraph(nodeStream, edgeStream);
remicres's avatar
remicres committed
831

832
  fclose(nodeStream);
remicres's avatar
remicres committed
833 834 835
  fclose(edgeStream);
}

836 837 838 839 840 841 842 843 844 845
template<class TSegmenter>
void ReadGraph(typename TSegmenter::GraphType& graph,
    const std::string& nodesPath,
    const std::string& edgesPath)
{
  TSegmenter segmenter;
  ReadGraph<TSegmenter>(segmenter, nodesPath, edgesPath);
  graph = segmenter.m_Graph;
}

remicres's avatar
remicres committed
846 847 848 849 850 851 852 853 854 855 856 857 858
// Write stability margin
template<class TSegmenter>
void WriteStabilityMargin(std::unordered_map<
    typename TSegmenter::NodePointerType,
    unsigned int>& stabilityMargin,
    const std::string& nodesPath,
    const std::string& edgesPath)
{
  FILE * nodeStream = fopen(nodesPath.c_str(), "wb");
  assert(nodeStream != NULL);
  FILE * edgeStream = fopen(edgesPath.c_str(), "wb");
  assert(edgeStream != NULL);

859 860 861
  // Write number of nodes
  std::size_t size = stabilityMargin.size();
  fwrite(&size, sizeof(size), 1, nodeStream);
remicres's avatar
remicres committed
862

863
  TSegmenter seg;
remicres's avatar
remicres committed
864 865 866
  for(auto& kv : stabilityMargin)
    {
    auto node = kv.first;
867
    seg.WriteNode(node, nodeStream);
remicres's avatar
remicres committed
868 869 870 871 872 873 874 875 876 877 878 879 880 881

    // Write only edges pointing to nodes which are in the stability margin.
    fwrite(&(node->m_Id), sizeof(node->m_Id), 1, edgeStream);
    std::size_t edgeSize = node->m_Edges.size();
    for(auto& edg : node->m_Edges)
      {
      if(stabilityMargin.find(edg.GetRegion()) == stabilityMargin.end())
        edgeSize--;
      }
    fwrite(&(edgeSize), sizeof(edgeSize), 1, edgeStream);
    for(auto& edg : node->m_Edges)
      {
      if(stabilityMargin.find(edg.GetRegion()) != stabilityMargin.end())
        {
882
        seg.WriteEdge(edg, edgeStream);
remicres's avatar
remicres committed
883 884 885 886 887 888 889 890
        }
      }
    }

  fclose(nodeStream);
  fclose(edgeStream);
}

891
// Write the graph
remicres's avatar
remicres committed
892 893
template<class TSegmenter>
void WriteGraph(typename TSegmenter::GraphType& graph,
894 895
    const std::string& nodeFile,
    const std::string& edgeFile)
remicres's avatar
remicres committed
896 897 898 899 900 901 902
{

  FILE * nodeStream = fopen(nodeFile.c_str(), "wb");
  assert(nodeStream != NULL);
  FILE * edgeStream = fopen(edgeFile.c_str(), "wb");
  assert(edgeStream != NULL);

903 904 905
  // Write number of nodes
  std::size_t size = graph.m_Nodes.size();
  fwrite(&size, sizeof(size), 1, nodeStream);
remicres's avatar
remicres committed
906

907
  TSegmenter seg;
remicres's avatar
remicres committed
908 909
  for(auto& node : graph.m_Nodes)
    {
910
    seg.WriteNode(node, nodeStream);
remicres's avatar
remicres committed
911 912 913 914 915 916 917

    // Write edges
    fwrite(&(node->m_Id), sizeof(node->m_Id), 1, edgeStream);
    std::size_t edgeSize = node->m_Edges.size();
    fwrite(&(edgeSize), sizeof(edgeSize), 1, edgeStream);
    for(auto& edg : node->m_Edges)
      {
918
      seg.WriteEdge(edg, edgeStream);
remicres's avatar
remicres committed
919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978
      }
    }

  fclose(nodeStream);
  fclose(edgeStream);
}

template<class TSegmenter>
void RemoveUnstableSegments(typename TSegmenter::GraphType& graph,
    ProcessingTile& tile,
    const unsigned int imageWidth)
{

  unsigned int rowPixel, colPixel;
  bool stable;

  for(auto& node : graph.m_Nodes)
    {
    if(node->m_Bbox.m_UX >= tile.columns[0] &&
        node->m_Bbox.m_UY >= tile.rows[0] &&
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 <= tile.columns[1] &&
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 <= tile.rows[1])
      continue;
    else if(node->m_Bbox.m_UX > tile.columns[1] ||
        node->m_Bbox.m_UY > tile.rows[1] ||
        node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[0] ||
        node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[0])
      {
      node->m_Expired = true;
      RemoveEdgeToUnstableNode<TSegmenter>(node);
      }
    else
      {
      lp::CellLists borderpixels;
      lp::ContourOperations::GenerateBorderCells(borderpixels, node->m_Contour, node->m_Id, imageWidth);
      stable = false;

      for(auto& pix : borderpixels)
        {
        rowPixel = pix / imageWidth;
        colPixel = pix % imageWidth;

        if(rowPixel >= tile.rows[0] &&
            rowPixel <= tile.rows[1] &&
            colPixel >= tile.columns[0] &&
            colPixel <= tile.columns[1])
          {
          stable = true;
          break;
          }
        }

      if(!stable)
        {
        node->m_Expired = true;
        RemoveEdgeToUnstableNode<TSegmenter>(node);
        }
      }
    }

979
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
980 981 982 983 984 985 986 987
}

template<class TSegmenter>
void RemoveEdgeToUnstableNode(typename TSegmenter::NodePointerType nodePtr)
{
  for(auto& edg : nodePtr->m_Edges)
    {
    auto nodeNeighbor = edg.GetRegion();
988
    auto EdgeToNode = grm::GraphOperations<TSegmenter>::FindEdge(nodeNeighbor, nodePtr);
remicres's avatar
remicres committed
989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007
    assert(EdgeToNode != nodeNeighbor->m_Edges.end());
    nodeNeighbor->m_Edges.erase(EdgeToNode);
    }
}

template<class TSegmenter>
void RescaleGraph(typename TSegmenter::GraphType& graph,
    ProcessingTile& tile,
    const unsigned int rowTile,
    const unsigned int colTile,
    const unsigned int tileWidth,
    const unsigned int tileHeight,
    const unsigned int imageWidth)
{
  unsigned int rowNodeTile, colNodeTile;
  unsigned int rowNodeImg, colNodeImg;

  for(auto& node : graph.m_Nodes)
    {
1008
    // Start pixel index of the node (in the tile)
1009 1010
    rowNodeTile = node->m_Id / tile.region.GetSize()[0];
    colNodeTile = node->m_Id % tile.region.GetSize()[0];
1011 1012

    // Start pixel index of the node (in the image)
1013 1014
    rowNodeImg = rowTile * tileHeight + rowNodeTile - tile.margin[0];
    colNodeImg = colTile * tileWidth + colNodeTile - tile.margin[3];
1015
    node->m_Id = (std::size_t)rowNodeImg * (std::size_t)imageWidth + (std::size_t)colNodeImg;
remicres's avatar
remicres committed
1016 1017

    // Change also its bounding box
1018 1019
    node->m_Bbox.m_UX = colTile * tileWidth + node->m_Bbox.m_UX - tile.margin[3];
    node->m_Bbox.m_UY = rowTile * tileHeight + node->m_Bbox.m_UY - tile.margin[0];
remicres's avatar
remicres committed
1020
    }
1021

remicres's avatar
remicres committed
1022 1023
}
}
1024
#endif