lsgrmGraphOperations.txx 33.4 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
//#include <unistd.h>
5
#include <fstream>
6
#include <cstdio>
7
#include <sys/stat.h>
8
9
10
11

namespace lsgrm
{

12
13
14
15
16
17
bool file_exists(const std::string& fileName)
{
    std::ifstream infile(fileName.c_str());
    return infile.good();
}

18
19
20
21
22
23
24
time_t last_mod_time(const std::string& fileName)
{
	struct stat result;
	if (!stat(fileName.c_str(), &result)) return result.st_mtime;
	else return -1;
}

remicres's avatar
remicres committed
25
26
27
28
29
template<class TSegmenter>
typename TSegmenter::ImageType::Pointer ReadImageRegion(
    typename TSegmenter::ImageType * inputPtr,
    typename TSegmenter::ImageType::RegionType region)
{
30

remicres's avatar
remicres committed
31
32
33
34
35
36
37
  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();
38

remicres's avatar
remicres committed
39
40
  return filter->GetOutput();
}
41

remicres's avatar
remicres committed
42
template<class TSegmenter>
43
typename TSegmenter::GraphType
remicres's avatar
remicres committed
44
MergeAllGraphsAndAchieveSegmentation(
45
    const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
46
47
48
49
50
51
    const float& threshold,
    std::vector<ProcessingTile>& tiles,
    const unsigned int nbTilesX,
    const unsigned int nbTilesY,
    const unsigned int imageWidth,
    const unsigned int imageHeight,
52
53
    const unsigned int imageBands,
    unsigned int numberOfIterations)
remicres's avatar
remicres committed
54
{
55

56
57
  std::cout << "--- Graph aggregation...\n" << std::endl;

remicres's avatar
remicres committed
58
  TSegmenter segmenter;
59

remicres's avatar
remicres committed
60
61
  std::cout << "Reading graphs" << std::endl;

remicres's avatar
remicres committed
62
63
64
65
  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; col++)
      {
remicres's avatar
remicres committed
66
67
      std::cout << "\tImporting graph of tile " << (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) << std::endl;

68
      InsertNodesFromTile<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], false);
remicres's avatar
remicres committed
69
70
      }
    }
71

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

remicres's avatar
remicres committed
74
75
76
77
78
  // 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
79
      std::cout << "Cleaning nodes of tile " << (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) << std::endl;
80

81
      std::unordered_map<std::size_t,
remicres's avatar
remicres committed
82
      std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
83

remicres's avatar
remicres committed
84
      std::cout << "\tBuildBorderPixelMap..." << std::endl;
85

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

remicres's avatar
remicres committed
89
      std::cout << "\tRemoveDuplicatedNodes..." << std::endl;
90

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

remicres's avatar
remicres committed
93
      std::cout << "\tUpdateNeighborsOfNoneDuplicatedNodes..." << std::endl;
94

remicres's avatar
remicres committed
95
      UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap, imageWidth, imageHeight);
remicres's avatar
remicres committed
96
97
98
      }
    }

remicres's avatar
remicres committed
99
  std::cout << "Achieve segmentation process" ;
remicres's avatar
remicres committed
100
101
102
103
104
105
106

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

111
  // patch raf (sort using node id directly --> seems ok)
112
  // Sort the nodes top-->down and left-->right
113
  /*
114
115
116
117
118
119
  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);
120

121
122
123
124
	// patch raf (were unsigned int)
    std::size_t offA = 0;
    std::size_t offB = 0;
	// end patch
125
126
127
    for (auto& pix: borderPixelsA)
      if (pix>offA)
        offA=pix;
128

129
130
131
    for (auto& pix: borderPixelsB)
      if (pix>offB)
        offB=pix;
132

133
134
      return offA > offB;
  });
135
136
137
138
139
140
141
  */
  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
142

143
  return segmenter.m_Graph;
144
}
145

remicres's avatar
remicres committed
146
template<class TSegmenter>
147
long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
148
149
150
151
152
153
154
155
    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,
156
157
    bool& isFusion,
	bool resume,
158
	unsigned int& nextTile)
remicres's avatar
remicres committed
159
160
161
162
163
164
165
166
167
168
169
170
171
{
  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++)
      {
        {
172
173
		// Get the current tile
        std::cout << "Processing tile " << row << ", " << col << std::endl;
174
175
        unsigned int currTile = row*nbTilesX + col;
        ProcessingTile currentTile = tiles[currTile];
176
		
177
		// If resume mode, get accumulated memory on previous tiles
178
179
180
181
182
183
184
185
186
		if (resume && currTile < nextTile ) {
			// Load the graph
			std::cout << "\tResuming graph..." << std::endl;
			TSegmenter segmenter;
			ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
			// Retrieve the amount of memory to store this graph
			std::cout << "\tGet graph memory..." << std::endl;
			accumulatedMemory += segmenter.GetGraphMemory();
			continue;
187
188
		}
		
remicres's avatar
remicres committed
189
        // Load the graph
190
191
        std::cout << "\tLoad graph..." << std::endl;
        TSegmenter segmenter;
192
        ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
remicres's avatar
remicres committed
193
194
195
196

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

          std::cout << "\tBuild border pixel map..." << std::endl;
201
          std::unordered_map<std::size_t, std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
202
          BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, currentTile, row, col,
remicres's avatar
remicres committed
203
204
205
206
207
208
209
210
211
212
213
              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;
214
          RemoveUselessNodes<TSegmenter>(currentTile, segmenter.m_Graph,
215
              imageWidth, numberOfNeighborLayers);
remicres's avatar
remicres committed
216
217
218
219
220
221
222
223
224

        }

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

228
        std::cout << "\tPartial segmentation";
229
        auto merge = grm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
remicres's avatar
remicres committed
230
231
232
233
234

        if(merge == true)
          isFusion = true;

        // Remove unstable segments
235
        std::cout << "\tRemove unstable segments..." << std::endl;
236
        RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
remicres's avatar
remicres committed
237
238

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

242
        // Write graph to temporay directory
243
        std::cout << "\tWrite graph..." << std::endl;
244
        WriteGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
245
		
remicres's avatar
remicres committed
246
247
248
        }
      }
    }
249
250
251
    
    if (resume)
		nextTile = 0;
remicres's avatar
remicres committed
252

253

remicres's avatar
remicres committed
254
255
256
257
258
259
260
261
  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++)
      {
        {
262
263
        // Get current tile
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
remicres's avatar
remicres committed
264
265

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

269
        // Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight-1
remicres's avatar
remicres committed
270
271
272
273
        // and write them to the stability margin
        {
          std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;

274
          DetectBorderNodes<TSegmenter>(graph, currentTile,
remicres's avatar
remicres committed
275
276
277
278
              borderNodeMap, imageWidth, imageHeight);

          ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);

279
          WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
remicres's avatar
remicres committed
280
281
282
283
284
285
        }
        std::cout << "Process finished on tile " << (row*nbTilesX + col) << std::endl;
        }
      }
    }
  std::cout << std::endl;
286
  
remicres's avatar
remicres committed
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
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
  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;
    }

363
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
364
365
366
}

template<class TSegmenter>
367
void UpdateNeighborsOfNoneDuplicatedNodes(std::unordered_map<std::size_t,
remicres's avatar
remicres committed
368
369
370
371
372
373
374
375
376
    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)
    {
377
    NeighIDType neighborhood[4];
378
    grm::FOURNeighborhood(neighborhood, pn.first, imageWidth, imageHeight);
remicres's avatar
remicres committed
379
380
381
382
383
384
385
386
387
388
389
390
391

    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)
            {
392
            auto toNeigh = grm::GraphOperations<TSegmenter>::FindEdge(currNode, neigh);
remicres's avatar
remicres committed
393
394
395
396
397
398
399
400
401
402
403
404

            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())
                  {
405
                  NeighIDType pixNeighborhood[4];
406
                  grm::FOURNeighborhood(pixNeighborhood, pix, imageWidth, imageHeight);
remicres's avatar
remicres committed
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431

                  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>
432
void RemoveDuplicatedNodes(std::unordered_map<std::size_t,
remicres's avatar
remicres committed
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
    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();

456
          auto toNit = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, *nit);
remicres's avatar
remicres committed
457
458
459
460
461
462
463

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

          boundary = edg.m_Boundary;

          neighNit->m_Edges.erase(toNit);

464
          auto toRefNode = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, refNode);
remicres's avatar
remicres committed
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491

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

492
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
493
494
495
496
497
498
499
500
501
502
}

// 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,
503
    std::unordered_map<std::size_t,
remicres's avatar
remicres committed
504
    std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
505
506
    const unsigned int imageWidth,
	bool merging)
remicres's avatar
remicres committed
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
{
  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;
      }
523
524
525
526
527
528
529
530
531
	// 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
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
    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;
        }
      }
    }
}

565
566
template<class TSegmenter>
void InsertNodesFromTile(typename TSegmenter::GraphType& graph,
567
    ProcessingTile& tile, bool margin)
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
{
  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
583
584
template<class TSegmenter>
void AddStabilityMargin(typename TSegmenter::GraphType& graph,
585
    std::vector<ProcessingTile>& tiles,
remicres's avatar
remicres committed
586
587
588
    const unsigned int row,
    const unsigned int col,
    const unsigned int nbTilesX,
remicres's avatar
remicres committed
589
    const unsigned int nbTilesY)
remicres's avatar
remicres committed
590
591
592
593
594
{

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

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

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

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

  // Margin to retrieve at top right
  if(row > 0 && col < nbTilesX - 1)
    {
619
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + (col+1)]);
remicres's avatar
remicres committed
620
621
622
623
624
    }

  // Margin to retrieve at bottom right
  if(row < nbTilesY - 1 && col < nbTilesX - 1)
    {
625
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + (col+1)]);
remicres's avatar
remicres committed
626
627
628
629
630
    }

  // Margin to retrieve at bottom left
  if(row < nbTilesY - 1 && col > 0)
    {
631
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + (col-1)]);
remicres's avatar
remicres committed
632
633
634
635
636
    }

  // Margin to retrieve at top left
  if(row > 0 && col > 0)
    {
637
    InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + (col-1)]);
remicres's avatar
remicres committed
638
639
640
641
642
643
    }
}

template<class TSegmenter>
long long unsigned int RunFirstPartialSegmentation(
    typename TSegmenter::ImageType * inputPtr,
644
    const typename TSegmenter::ParamType& params,
remicres's avatar
remicres committed
645
646
647
648
649
650
651
652
    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,
653
    bool& isFusion,
654
655
    bool resume,
    unsigned int& nextTile)
remicres's avatar
remicres committed
656
657
658
{

  using ImageType = typename TSegmenter::ImageType;
659
660
  const unsigned int imageWidth = inputPtr->GetLargestPossibleRegion().GetSize()[0];
  const unsigned int imageHeight = inputPtr->GetLargestPossibleRegion().GetSize()[1];
remicres's avatar
remicres committed
661
662
663

  long long unsigned int accumulatedMemory = 0;
  isFusion = false;
664
  bool accomplished = true;
remicres's avatar
remicres committed
665
666
667
668
669

  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;

670
671
  time_t lastTime = -1;

Gaetano Raffaele's avatar
Gaetano Raffaele committed
672
  for(unsigned int row = 0; row < nbTilesY; ++row)
remicres's avatar
remicres committed
673
    {
Gaetano Raffaele's avatar
Gaetano Raffaele committed
674
    for(unsigned int col = 0; col < nbTilesX; col++)
remicres's avatar
remicres committed
675
676
677
678
      {
        {
        // Reading images
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
679
680
681
        
        if (resume && file_exists(currentTile.nodeFileName) && file_exists(currentTile.edgeFileName))
		{
682
683
684
685
686
			int lastMod = last_mod_time(currentTile.nodeFileName);
			if (lastMod > lastTime) {
				lastTime = lastMod;
				nextTile = row*nbTilesX + col + 1;
			}
687
688
689
690
691
692
693
694
695
696
			std::cout << "\tResuming graph..." << std::endl;
			TSegmenter segmenter;
			ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
			// Retrieve the amount of memory to store this graph
			std::cout << "\tGet graph memory..." << std::endl;
			accumulatedMemory += segmenter.GetGraphMemory();
			if(segmenter.GetComplete() == false)
				isFusion = true;
			continue;
		}
697
698
699
		
		accomplished = false;
		
700
        std::cout << "Processing tile " <<  (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) <<
remicres's avatar
remicres committed
701
702
703
704
705
706
707
            " (" << 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
708
        std::cout << "\tSegmenting";
remicres's avatar
remicres committed
709
710
711
        TSegmenter segmenter;
        segmenter.SetParam(params);
        segmenter.SetThreshold(threshold);
712
        segmenter.SetDoFastSegmentation(false);
remicres's avatar
remicres committed
713
714
715
        segmenter.SetNumberOfIterations(niter);
        segmenter.SetInput(imageTile);
        segmenter.Update();
716
		
remicres's avatar
remicres committed
717
718
719
720
        if(segmenter.GetComplete() == false)
          isFusion = true;

        // Rescale the graph to be in the reference of the image
721
        std::cout << "\tRescaling graph..." << std::endl;
remicres's avatar
remicres committed
722
        RescaleGraph<TSegmenter>(segmenter.m_Graph,
723
            currentTile,
remicres's avatar
remicres committed
724
725
726
727
728
            row,
            col,
            tileWidth,
            tileHeight,
            imageWidth);
729
			
remicres's avatar
remicres committed
730
        // Remove unstable segments
731
        std::cout << "\tRemoving unstable segments..." << std::endl;
732
        RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
remicres's avatar
remicres committed
733
734

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

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

742
        // Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight -1
remicres's avatar
remicres committed
743
        // and write them to the stability margin
744
        std::cout << "\tComputing stability margin..." << std::endl;
remicres's avatar
remicres committed
745
        {
remicres's avatar
remicres committed
746
        std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;
remicres's avatar
remicres committed
747

remicres's avatar
remicres committed
748
749
        DetectBorderNodes<TSegmenter>(segmenter.m_Graph, currentTile,
            borderNodeMap, imageWidth, imageHeight);
remicres's avatar
remicres committed
750

remicres's avatar
remicres committed
751
        ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
remicres's avatar
remicres committed
752

remicres's avatar
remicres committed
753
        WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
754
        
remicres's avatar
remicres committed
755
756
757
758
        }
        }
      } // for each col
    } // for each row
759
760
761
762
763
764
765
    
    if (resume && !accomplished) {
		std::cout << "\tSegmentation recovered during first partial segmentation." << std::endl;
		nextTile = nbTilesX*nbTilesY;
	}
	if (resume && accomplished && nextTile == nbTilesX*nbTilesY)
		nextTile = 0;
remicres's avatar
remicres committed
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
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872

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

873
874

// Generic
remicres's avatar
remicres committed
875
template<class TSegmenter>
876
void ReadGraph(TSegmenter& segmenter,
remicres's avatar
remicres committed
877
878
879
    const std::string& nodesPath,
    const std::string& edgesPath)
{
880
  FILE * nodeStream = fopen(nodesPath.c_str(), "rb");
remicres's avatar
remicres committed
881
  assert(nodeStream != NULL);
882
  FILE * edgeStream = fopen(edgesPath.c_str(), "rb");
remicres's avatar
remicres committed
883
884
  assert(edgeStream != NULL);

885
  segmenter.ReadGraph(nodeStream, edgeStream);
remicres's avatar
remicres committed
886

887
  fclose(nodeStream);
remicres's avatar
remicres committed
888
889
890
  fclose(edgeStream);
}

891
892
893
894
895
896
897
898
899
900
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
901
902
903
904
905
906
907
908
909
910
911
912
913
// 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);

914
915
916
  // Write number of nodes
  std::size_t size = stabilityMargin.size();
  fwrite(&size, sizeof(size), 1, nodeStream);
remicres's avatar
remicres committed
917

918
  TSegmenter seg;
remicres's avatar
remicres committed
919
920
921
  for(auto& kv : stabilityMargin)
    {
    auto node = kv.first;
922
    seg.WriteNode(node, nodeStream);
remicres's avatar
remicres committed
923
924
925
926
927
928
929
930
931
932
933
934
935
936

    // 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())
        {
937
        seg.WriteEdge(edg, edgeStream);
remicres's avatar
remicres committed
938
939
940
941
942
943
944
945
        }
      }
    }

  fclose(nodeStream);
  fclose(edgeStream);
}

946
// Write the graph
remicres's avatar
remicres committed
947
948
template<class TSegmenter>
void WriteGraph(typename TSegmenter::GraphType& graph,
949
950
    const std::string& nodeFile,
    const std::string& edgeFile)
remicres's avatar
remicres committed
951
952
953
954
955
956
957
{

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

958
959
960
  // Write number of nodes
  std::size_t size = graph.m_Nodes.size();
  fwrite(&size, sizeof(size), 1, nodeStream);
remicres's avatar
remicres committed
961

962
  TSegmenter seg;
remicres's avatar
remicres committed
963
964
  for(auto& node : graph.m_Nodes)
    {
965
    seg.WriteNode(node, nodeStream);
remicres's avatar
remicres committed
966
967
968
969
970
971
972

    // 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)
      {
973
      seg.WriteEdge(edg, edgeStream);
remicres's avatar
remicres committed
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
      }
    }

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

1034
  grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
remicres's avatar
remicres committed
1035
1036
1037
1038
1039
1040
1041
1042
}

template<class TSegmenter>
void RemoveEdgeToUnstableNode(typename TSegmenter::NodePointerType nodePtr)
{
  for(auto& edg : nodePtr->m_Edges)
    {
    auto nodeNeighbor = edg.GetRegion();
Gaetano Raffaele's avatar
Revert.    
Gaetano Raffaele committed
1043
1044
1045
    auto EdgeToNode = grm::GraphOperations<TSegmenter>::FindEdge(nodeNeighbor, nodePtr);
    assert(EdgeToNode != nodeNeighbor->m_Edges.end());
    nodeNeighbor->m_Edges.erase(EdgeToNode);
remicres's avatar
remicres committed
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
    }
}

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)
    {
1063
    // Start pixel index of the node (in the tile)
1064
1065
    rowNodeTile = node->m_Id / tile.region.GetSize()[0];
    colNodeTile = node->m_Id % tile.region.GetSize()[0];
1066
1067

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

    // Change also its bounding box
1073
1074
    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
1075
    }
1076

remicres's avatar
remicres committed
1077
1078
}
}
1079
#endif