lsgrmGraphOperations.txx 31 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
72
      BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], row, col,
          nbTilesX, nbTilesY, borderPixelMap, imageWidth);
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
97
98
99
100
101
102
  // Sort the nodes top-->down and left-->right
  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);
103

104
105
106
107
108
    unsigned int offA = 0;
    unsigned int offB = 0;
    for (auto& pix: borderPixelsA)
      if (pix>offA)
        offA=pix;
109

110
111
112
    for (auto& pix: borderPixelsB)
      if (pix>offB)
        offB=pix;
113

114
115
      return offA > offB;
  });
116

117
  return segmenter.m_Graph;
118
}
119

remicres's avatar
remicres committed
120
template<class TSegmenter>
121
long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    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++)
      {
143
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
144
      if (MyTurn(row*nbTilesX + col))
145
#endif
remicres's avatar
remicres committed
146
        {
147
        // Get the current tile
148
        std::cout << "Processing tile " << row << ", " << col << std::endl;
149
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
150

remicres's avatar
remicres committed
151
        // Load the graph
152
153
        std::cout << "\tLoad graph..." << std::endl;
        TSegmenter segmenter;
154
        ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
remicres's avatar
remicres committed
155
156
157
158

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

          std::cout << "\tBuild border pixel map..." << std::endl;
163
          std::unordered_map<std::size_t, std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
164
          BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, currentTile, row, col,
remicres's avatar
remicres committed
165
166
167
168
169
170
171
172
173
174
175
              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;
176
          RemoveUselessNodes<TSegmenter>(currentTile, segmenter.m_Graph,
177
              imageWidth, numberOfNeighborLayers);
remicres's avatar
remicres committed
178
179
180
181
182
183
184
185
186

        }

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

190
        std::cout << "\tPartial segmentation";
191
        auto merge = grm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
remicres's avatar
remicres committed
192
193
194
195
196

        if(merge == true)
          isFusion = true;

        // Remove unstable segments
197
        std::cout << "\tRemove unstable segments..." << std::endl;
198
        RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
remicres's avatar
remicres committed
199
200

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

204
        // Write graph to temporay directory
205
        std::cout << "\tWrite graph..." << std::endl;
206
        WriteGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
remicres's avatar
remicres committed
207
208
209
210
211
212
        }
      }
    }

#ifdef OTB_USE_MPI
  otb::MPIConfig::Instance()->barrier();
213
214
#endif

remicres's avatar
remicres committed
215
216
217
218
219
220
221
  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++)
      {
222
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
223
      if (MyTurn(row*nbTilesX + col))
224
#endif
remicres's avatar
remicres committed
225
        {
226
227
        // Get current tile
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
remicres's avatar
remicres committed
228
229

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

233
        // Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight-1
remicres's avatar
remicres committed
234
235
236
237
        // and write them to the stability margin
        {
          std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;

238
          DetectBorderNodes<TSegmenter>(graph, currentTile,
remicres's avatar
remicres committed
239
240
241
242
              borderNodeMap, imageWidth, imageHeight);

          ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);

243
          WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
remicres's avatar
remicres committed
244
245
246
247
248
249
250
251
252
253
254
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
        }
        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;
    }

327
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
328
329
330
}

template<class TSegmenter>
331
void UpdateNeighborsOfNoneDuplicatedNodes(std::unordered_map<std::size_t,
remicres's avatar
remicres committed
332
333
334
335
336
337
338
339
340
    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)
    {
341
    NeighIDType neighborhood[4];
342
    grm::FOURNeighborhood(neighborhood, pn.first, imageWidth, imageHeight);
remicres's avatar
remicres committed
343
344
345
346
347
348
349
350
351
352
353
354
355

    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)
            {
356
            auto toNeigh = grm::GraphOperations<TSegmenter>::FindEdge(currNode, neigh);
remicres's avatar
remicres committed
357
358
359
360
361
362
363
364
365
366
367
368

            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())
                  {
369
                  NeighIDType pixNeighborhood[4];
370
                  grm::FOURNeighborhood(pixNeighborhood, pix, imageWidth, imageHeight);
remicres's avatar
remicres committed
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395

                  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>
396
void RemoveDuplicatedNodes(std::unordered_map<std::size_t,
remicres's avatar
remicres committed
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
    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();

420
          auto toNit = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, *nit);
remicres's avatar
remicres committed
421
422
423
424
425
426
427

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

          boundary = edg.m_Boundary;

          neighNit->m_Edges.erase(toNit);

428
          auto toRefNode = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, refNode);
remicres's avatar
remicres committed
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455

          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);
          }
        }
      }
    }

456
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
457
458
459
460
461
462
463
464
465
466
}

// 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,
467
    std::unordered_map<std::size_t,
remicres's avatar
remicres committed
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
    std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
    const unsigned int imageWidth)
{
  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;
      }
    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;
        }
      }
    }
}

519
520
template<class TSegmenter>
void InsertNodesFromTile(typename TSegmenter::GraphType& graph,
521
    ProcessingTile& tile, bool margin)
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
{
  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
537
538
template<class TSegmenter>
void AddStabilityMargin(typename TSegmenter::GraphType& graph,
539
    std::vector<ProcessingTile>& tiles,
remicres's avatar
remicres committed
540
541
542
    const unsigned int row,
    const unsigned int col,
    const unsigned int nbTilesX,
remicres's avatar
remicres committed
543
    const unsigned int nbTilesY)
remicres's avatar
remicres committed
544
545
546
547
548
{

  // Margin to retrieve at top
  if(row > 0)
    {
549
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + col]);
remicres's avatar
remicres committed
550
551
552
553
554
    }

  // Margin to retrieve at right
  if(col < nbTilesX - 1)
    {
555
    InsertNodesFromTile<TSegmenter>(graph, tiles[row * nbTilesX + (col+1)]);
remicres's avatar
remicres committed
556
557
558
559
560
    }

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

remicres's avatar
remicres committed
564
565
566
  // Margin to retrieve at left
  if(col > 0)
    {
567
    InsertNodesFromTile<TSegmenter>(graph, tiles[row * nbTilesX + (col-1)]);
remicres's avatar
remicres committed
568
569
570
571
572
    }

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

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

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

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

template<class TSegmenter>
long long unsigned int RunFirstPartialSegmentation(
    typename TSegmenter::ImageType * inputPtr,
598
    const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
599
600
601
602
603
604
605
606
607
608
609
610
    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;
611
612
  const unsigned int imageWidth = inputPtr->GetLargestPossibleRegion().GetSize()[0];
  const unsigned int imageHeight = inputPtr->GetLargestPossibleRegion().GetSize()[1];
remicres's avatar
remicres committed
613
614
615
616
617
618
619
620
621
622
623
624

  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;

  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; col++)
      {
625
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
626
      if (MyTurn(row*nbTilesX + col))
627
#endif
remicres's avatar
remicres committed
628
629
630
        {
        // Reading images
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
631

632
        std::cout << "Processing tile " <<  (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) <<
remicres's avatar
remicres committed
633
634
635
636
637
638
639
            " (" << 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
640
        std::cout << "\tSegmenting";
remicres's avatar
remicres committed
641
642
643
        TSegmenter segmenter;
        segmenter.SetParam(params);
        segmenter.SetThreshold(threshold);
644
        segmenter.SetDoFastSegmentation(false);
remicres's avatar
remicres committed
645
646
647
648
649
650
651
652
        segmenter.SetNumberOfIterations(niter);
        segmenter.SetInput(imageTile);
        segmenter.Update();

        if(segmenter.GetComplete() == false)
          isFusion = true;

        // Rescale the graph to be in the reference of the image
653
        std::cout << "\tRescaling graph..." << std::endl;
remicres's avatar
remicres committed
654
        RescaleGraph<TSegmenter>(segmenter.m_Graph,
655
            currentTile,
remicres's avatar
remicres committed
656
657
658
659
660
661
662
            row,
            col,
            tileWidth,
            tileHeight,
            imageWidth);

        // Remove unstable segments
663
        std::cout << "\tRemoving unstable segments..." << std::endl;
664
        RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
remicres's avatar
remicres committed
665
666

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

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

674
        // Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight -1
remicres's avatar
remicres committed
675
        // and write them to the stability margin
676
        std::cout << "\tComputing stability margin..." << std::endl;
remicres's avatar
remicres committed
677
        {
remicres's avatar
remicres committed
678
        std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;
remicres's avatar
remicres committed
679

remicres's avatar
remicres committed
680
681
        DetectBorderNodes<TSegmenter>(segmenter.m_Graph, currentTile,
            borderNodeMap, imageWidth, imageHeight);
remicres's avatar
remicres committed
682

remicres's avatar
remicres committed
683
        ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
remicres's avatar
remicres committed
684

remicres's avatar
remicres committed
685
        WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
remicres's avatar
remicres committed
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
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
        }
        }
      } // 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;
        }
      }
    }
}

797
798

// Generic
remicres's avatar
remicres committed
799
template<class TSegmenter>
800
void ReadGraph(TSegmenter& segmenter,
remicres's avatar
remicres committed
801
802
803
    const std::string& nodesPath,
    const std::string& edgesPath)
{
804
  FILE * nodeStream = fopen(nodesPath.c_str(), "rb");
remicres's avatar
remicres committed
805
  assert(nodeStream != NULL);
806
  FILE * edgeStream = fopen(edgesPath.c_str(), "rb");
remicres's avatar
remicres committed
807
808
  assert(edgeStream != NULL);

809
  segmenter.ReadGraph(nodeStream, edgeStream);
remicres's avatar
remicres committed
810

811
  fclose(nodeStream);
remicres's avatar
remicres committed
812
813
814
  fclose(edgeStream);
}

815
816
817
818
819
820
821
822
823
824
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
825
826
827
828
829
830
831
832
833
834
835
836
837
// 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);

838
839
840
  // Write number of nodes
  std::size_t size = stabilityMargin.size();
  fwrite(&size, sizeof(size), 1, nodeStream);
remicres's avatar
remicres committed
841

842
  TSegmenter seg;
remicres's avatar
remicres committed
843
844
845
  for(auto& kv : stabilityMargin)
    {
    auto node = kv.first;
846
    seg.WriteNode(node, nodeStream);
remicres's avatar
remicres committed
847
848
849
850
851
852
853
854
855
856
857
858
859
860

    // 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())
        {
861
        seg.WriteEdge(edg, edgeStream);
remicres's avatar
remicres committed
862
863
864
865
866
867
868
869
        }
      }
    }

  fclose(nodeStream);
  fclose(edgeStream);
}

870
// Write the graph
remicres's avatar
remicres committed
871
872
template<class TSegmenter>
void WriteGraph(typename TSegmenter::GraphType& graph,
873
874
    const std::string& nodeFile,
    const std::string& edgeFile)
remicres's avatar
remicres committed
875
876
877
878
879
880
881
{

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

882
883
884
  // Write number of nodes
  std::size_t size = graph.m_Nodes.size();
  fwrite(&size, sizeof(size), 1, nodeStream);
remicres's avatar
remicres committed
885

886
  TSegmenter seg;
remicres's avatar
remicres committed
887
888
  for(auto& node : graph.m_Nodes)
    {
889
    seg.WriteNode(node, nodeStream);
remicres's avatar
remicres committed
890
891
892
893
894
895
896

    // 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)
      {
897
      seg.WriteEdge(edg, edgeStream);
remicres's avatar
remicres committed
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
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
      }
    }

  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);
        }
      }
    }

958
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
959
960
961
962
963
964
965
966
}

template<class TSegmenter>
void RemoveEdgeToUnstableNode(typename TSegmenter::NodePointerType nodePtr)
{
  for(auto& edg : nodePtr->m_Edges)
    {
    auto nodeNeighbor = edg.GetRegion();
967
    auto EdgeToNode = grm::GraphOperations<TSegmenter>::FindEdge(nodeNeighbor, nodePtr);
remicres's avatar
remicres committed
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
    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)
    {
987
    // Start pixel index of the node (in the tile)
988
989
    rowNodeTile = node->m_Id / tile.region.GetSize()[0];
    colNodeTile = node->m_Id % tile.region.GetSize()[0];
990
991

    // Start pixel index of the node (in the image)
992
993
    rowNodeImg = rowTile * tileHeight + rowNodeTile - tile.margin[0];
    colNodeImg = colTile * tileWidth + colNodeTile - tile.margin[3];
remicres's avatar
remicres committed
994
995
996
    node->m_Id = rowNodeImg * imageWidth + colNodeImg;

    // Change also its bounding box
997
998
    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
999
    }
1000

remicres's avatar
remicres committed
1001
1002
}
}
1003
#endif