lsgrmController.txx 16.3 KB
Newer Older
1
2
#ifndef __LSGRM_CONTROLLER_TXX
#define __LSGRM_CONTROLLER_TXX
remicres's avatar
remicres committed
3
#include "lsgrmController.h"
4
5
6
7

namespace lsgrm
{

remicres's avatar
remicres committed
8
9
10
template<class TSegmenter>
Controller<TSegmenter>::Controller()
{
11
  m_TilingMode = LSGRM_TILING_AUTO;
remicres's avatar
remicres committed
12
13
14
15
16
17
18
19
20
  m_Margin = 0;
  m_NumberOfIterations = 0;
  m_NumberOfFirstIterations = 0;
  m_TileHeight = 0;
  m_TileWidth = 0;
  m_NbTilesX = 0;
  m_NbTilesY = 0;
  m_Threshold = 75;
  m_Memory = 0;
21
  m_Resuming = false;
22

remicres's avatar
remicres committed
23
24
25
26
27
28
29
}

template<class TSegmenter>
Controller<TSegmenter>::~Controller()
{
}

30
31
32
33
34
35
36
template<class TSegmenter>
void Controller<TSegmenter>::Modified()
{
  Superclass::Modified();
  m_Tiles.clear();
}

37
38
39
40
/*
 * Run the segmentation
 * TODO: compute the correct number of iterations !
 */
remicres's avatar
remicres committed
41
42
43
template<class TSegmenter>
void Controller<TSegmenter>::RunSegmentation()
{
remicres's avatar
remicres committed
44
  itkDebugMacro(<< "Entering RunSegmentation()");
remicres's avatar
remicres committed
45

46
47
  CheckMemorySize();

48
  if (m_TilingMode == LSGRM_TILING_AUTO || m_TilingMode == LSGRM_TILING_USER)
remicres's avatar
remicres committed
49
    {
50
51
52
53
    if(m_TilingMode == LSGRM_TILING_AUTO)
      {
      this->GetAutomaticConfiguration();
      }
54
    else // m_TilingMode is LSGRM_TILING_USER
55
      {
remicres's avatar
remicres committed
56
57
      m_NbTilesX = std::floor(m_InputImage->GetLargestPossibleRegion().GetSize()[0] / m_TileWidth);
      m_NbTilesY = std::floor(m_InputImage->GetLargestPossibleRegion().GetSize()[1] / m_TileHeight);
58
59
60
61
      m_Margin = static_cast<unsigned int>(pow(2, m_NumberOfFirstIterations + 1) - 2);
      }

    std::cout <<
remicres's avatar
remicres committed
62
63
64
65
66
67
68
        "--- Configuration: " <<
        "\n\tAvailable RAM: " << m_Memory <<
        "\n\tInput image dimensions: " << m_InputImage->GetLargestPossibleRegion().GetSize() <<
        "\n\tNumber of first iterations: " << m_NumberOfFirstIterations <<
        "\n\tStability margin: " << m_Margin <<
        "\n\tRegular tile size: " << m_TileWidth << " x " << m_TileHeight <<
        "\n\tTiling layout: " << m_NbTilesX << " x " << m_NbTilesY << std::endl;
remicres's avatar
remicres committed
69
70

    // Compute the splitting scheme
71
72
    m_Tiles = SplitOTBImage<ImageType>(m_InputImage, m_TileWidth, m_TileHeight, m_Margin,
        m_NbTilesX, m_NbTilesY, m_TemporaryFilesPrefix);
remicres's avatar
remicres committed
73

74
    // If there is only one tile, then fallback to LSGRM_TILING_NONE case
75
76
77
78
79
    if (m_Tiles.size() == 1)
      {
      std::cout << "Only one tile is needed. Fallback to tiling=none." << std::endl;
      SetTilingModeNone();
      }
80
81
82
83
    }

  if (m_TilingMode == LSGRM_TILING_AUTO || m_TilingMode == LSGRM_TILING_USER)
    {
84
    const unsigned int numberOfIterationsForPartialSegmentations = 3; // TODO: find a smart value
85
86
    unsigned int numberOfIterationsRemaining = m_NumberOfIterations;

remicres's avatar
remicres committed
87
88
89
90
91
    // Boolean indicating if there are remaining fusions
    bool isFusion = false;

    // Run first partial segmentation
    boost::timer t; t.restart();
92
	
93
94
	// temp. patch, maybe calculate real current memory after resuming graphs.
	long long unsigned int accumulatedMemory = 2 * m_Memory;
95
	unsigned int nextTile = m_NbTilesX*m_NbTilesY;
96
97
	
	accumulatedMemory = RunFirstPartialSegmentation<TSegmenter>(
remicres's avatar
remicres committed
98
99
100
101
        m_InputImage,
        m_SpecificParameters,
        m_Threshold,
        m_NumberOfFirstIterations,
102
        numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
103
104
105
106
107
        m_Tiles,
        m_NbTilesX,
        m_NbTilesY,
        m_TileWidth,
        m_TileHeight,
108
        isFusion,
109
110
111
        m_Resuming,
        nextTile);
        
112
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
113
    GatherUsefulVariables(accumulatedMemory, isFusion);
114
#endif
remicres's avatar
remicres committed
115
116
117
118

    // Time monitoring
    ShowTime(t);

119
120
121
122
123
124
125
126
127
128
	if (m_Resuming & nextTile < m_NbTilesX*m_NbTilesY) {
		std::cout << "Detected resumable process from tile #" << nextTile << std::endl;
		std::cout << "Forcing regular pass on the whole graph to update stability margins." << std::endl;
	}
	else {
		std::cout << "Process resumed at first stage." << std::endl;
		StopResumingMode();
	}

    while(m_Resuming || (accumulatedMemory > m_Memory && isFusion))
remicres's avatar
remicres committed
129
      {
remicres's avatar
remicres committed
130
131
      isFusion = false;
      accumulatedMemory = RunPartialSegmentation<TSegmenter>(
remicres's avatar
remicres committed
132
133
          m_SpecificParameters,
          m_Threshold,
134
          numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
135
136
137
138
139
140
          m_Tiles,
          m_NbTilesX,
          m_NbTilesY,
          m_InputImage->GetLargestPossibleRegion().GetSize()[0],
          m_InputImage->GetLargestPossibleRegion().GetSize()[1],
          m_InputImage->GetNumberOfComponentsPerPixel(),
141
142
          isFusion,
		  m_Resuming,
143
		  nextTile);
144
145
      
	  if (m_Resuming) StopResumingMode();
remicres's avatar
remicres committed
146

147
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
148
      GatherUsefulVariables(accumulatedMemory, isFusion);
149
#endif
150

remicres's avatar
remicres committed
151
152
      // Time monitoring
      ShowTime(t);
153

154
155
156
      // Update number of remaining iterations
      if (numberOfIterationsRemaining < numberOfIterationsForPartialSegmentations)
        {
157
        break;
158
159
160
161
162
        }
      else
        {
        numberOfIterationsRemaining -= numberOfIterationsForPartialSegmentations;
        }
remicres's avatar
remicres committed
163
      }
164
	
remicres's avatar
remicres committed
165
#ifdef OTB_USE_MPI
166
    // Only the master process is doing the next part
remicres's avatar
remicres committed
167
    // TODO: Use the MPI process wich has the largest amount of memory
remicres's avatar
remicres committed
168
169
170
171
172
173
174
    if (otb::MPIConfig::Instance()->GetMyRank() != 0)
      return;
#endif

    if(accumulatedMemory <= m_Memory)
      {
      // Merge all the graphs
175
        m_OutputGraph = MergeAllGraphsAndAchieveSegmentation<TSegmenter>(
remicres's avatar
remicres committed
176
177
178
179
180
181
182
          m_SpecificParameters,
          m_Threshold,
          m_Tiles,
          m_NbTilesX,
          m_NbTilesY,
          m_InputImage->GetLargestPossibleRegion().GetSize()[0],
          m_InputImage->GetLargestPossibleRegion().GetSize()[1],
183
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
184
          numberOfIterationsRemaining);
remicres's avatar
remicres committed
185

186
      ShowTime(t);
remicres's avatar
remicres committed
187
188

      }
189
    else // accumulatedMemory > m_Memory
remicres's avatar
remicres committed
190
      {
remicres's avatar
remicres committed
191
      // That means there are no more possible fusions but we can not store the output graph
remicres's avatar
remicres committed
192
193
      // Todo do not clean up temporary directory before copying resulting graph to the output directory
      // In the output directory add an info file to give the number of tiles.
194
      itkExceptionMacro(<< "No more possible fusions, but can not store the output graph");
remicres's avatar
remicres committed
195
196
      }
    }
remicres's avatar
remicres committed
197
  else if (m_TilingMode == LSGRM_TILING_NONE)// tiling_mode is none
remicres's avatar
remicres committed
198
    {
remicres's avatar
remicres committed
199
200
201
202
203
204
205
206
207
#ifdef OTB_USE_MPI
    // Only the master process is doing the next part
    if (otb::MPIConfig::Instance()->GetMyRank() > 0)
      return;
    else
      // Warn that there is some unused MPI processes
      if (otb::MPIConfig::Instance()->GetNbProcs() > 1)
        itkWarningMacro(<< "Only 1 MPI process will be used");
#endif
208
209
210
    // Update input image
    m_InputImage->Update();

remicres's avatar
remicres committed
211
212
213
214
    // Use classic grm
    TSegmenter segmenter;
    segmenter.SetParam(m_SpecificParameters);
    segmenter.SetThreshold(m_Threshold);
215
    segmenter.SetDoFastSegmentation(false);
remicres's avatar
remicres committed
216
217
218
    segmenter.SetNumberOfIterations(m_NumberOfIterations);
    segmenter.SetInput(m_InputImage);
    segmenter.Update();
219

220
    m_OutputGraph = segmenter.m_Graph;
remicres's avatar
remicres committed
221
222
223
224
    }
  else
    {
    itkExceptionMacro(<<"Unknow tiling mode!");
remicres's avatar
remicres committed
225
226
    }

227
228
  // TODO: [MPI] broadcast the graph to other nodes

remicres's avatar
remicres committed
229
230
}

231
232
/*
 * Compute the memory occupied by one node
remicres's avatar
remicres committed
233
234
 * TODO: compute the exact value, e.g. on a given UNIX system,
 * experimental measures shows that
235
 * for one Baatz node (+pixel) memory is about 700-730 bytes...
remicres's avatar
remicres committed
236
 * And our estimation is about 600
237
238
239
240
 */
template<class TSegmenter>
unsigned int Controller<TSegmenter>::GetNodeMemory()
{
241
  // Create a n*n image
242
  const unsigned int n = 100;
243
244
245
246
  typename ImageType::Pointer onePixelImage = ImageType::New();
  typename ImageType::IndexType start;
  start.Fill(0);
  typename ImageType::SizeType size;
247
  size.Fill(n);
248
249
250
251
  typename ImageType::RegionType region(start, size);
  onePixelImage->SetRegions(region);
  onePixelImage->SetNumberOfComponentsPerPixel(m_InputImage->GetNumberOfComponentsPerPixel());
  onePixelImage->Allocate();
252
253

  // Instanciate and initialize a segmenter
254
255
  TSegmenter segmenter;
  segmenter.SetInput(onePixelImage);
256
  grm::GraphOperations<TSegmenter>::InitNodes(onePixelImage,segmenter,FOUR);
257

258
  // Get the memory occupied by the graph, normalize it by n*n
259
  unsigned int memory = segmenter.GetGraphMemory() / (n*n);
260

261
  itkDebugMacro(<<"Size of a node is " << memory);
262
263

  // Get the memory occupied by one pixel of the image
remicres's avatar
remicres committed
264
265
  unsigned int pixelMemory =  sizeof(m_InputImage->GetBufferPointer())
      * m_InputImage->GetNumberOfComponentsPerPixel();
266

267
  itkDebugMacro(<<"Size of an image pixel is " << pixelMemory);
268
269
270

  memory += pixelMemory;

271
  itkDebugMacro(<<"Size of a node+pixel is " << memory);
272

273
  return memory;
274
}
remicres's avatar
remicres committed
275

remicres's avatar
remicres committed
276
template<class TSegmenter>
277
void Controller<TSegmenter>::CheckMemorySize()
remicres's avatar
remicres committed
278
{
279
280
281
282
283
  if (m_Memory == 0)
    {
    m_Memory = getMemorySize();
    assert(m_Memory > 0);
    }
remicres's avatar
remicres committed
284
  m_Memory /= 2; // For safety and can prevent out of memory troubles
285
286
287
288
289
}
/*
 * Compute the maximum number of nodes which can fit in the memory
 */
template<class TSegmenter>
290
std::size_t Controller<TSegmenter>::GetMaximumNumberOfNodesInMemory()
291
292
{
  itkDebugMacro(<< "Computing maximum number of nodes in memory");
remicres's avatar
remicres committed
293

294
  return std::ceil(((float) m_Memory) / ((float) GetNodeMemory()));
remicres's avatar
remicres committed
295
296
297
}

template<class TSegmenter>
298
void Controller<TSegmenter>::ComputeMaximumStabilityMargin(unsigned int width,
remicres's avatar
remicres committed
299
    unsigned int height, unsigned int &niter, unsigned int &margin)
300
    {
remicres's avatar
remicres committed
301
302
303
304
  itkDebugMacro(<< "Computing maximum stability margin");

  // Compute the stability margin. The naive strategy consider a margin value and a stable size equal.
  niter = 1;
305
306
307
  unsigned int maxMargin = std::min(width, height)/2;
  unsigned int currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
  margin = currMargin;
remicres's avatar
remicres committed
308

309
  while(currMargin < maxMargin)
remicres's avatar
remicres committed
310
    {
311
    margin = currMargin;
remicres's avatar
remicres committed
312
    niter++;
313
    currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
remicres's avatar
remicres committed
314
315
316
317
    }
  niter--;

  itkDebugMacro(<< "Number of iterations=" << niter << " margin=" << margin);
remicres's avatar
remicres committed
318

319
    }
remicres's avatar
remicres committed
320

321
/*
322
 * Compute a tiling layout which minimizes a criterion based on tile compactness
remicres's avatar
remicres committed
323
 * and memory usage
324
325
 *
 * TODO: use the lsgrmSplitter to truly compute the largest tile of a given layout
326
 */
remicres's avatar
remicres committed
327
328
329
330
template<class TSegmenter>
void Controller<TSegmenter>::GetAutomaticConfiguration()
{

remicres's avatar
remicres committed
331
332
  itkDebugMacro(<<"Get automatic configuration");

remicres's avatar
remicres committed
333
  // Compute the maximum number of nodes that can fit the memory
remicres's avatar
remicres committed
334
  // TODO: Use the smallest number amongst MPI processes
remicres's avatar
remicres committed
335
  unsigned long int maximumNumberOfNodesInMemory = GetMaximumNumberOfNodesInMemory();
remicres's avatar
remicres committed
336
  itkDebugMacro(<<"Maximum number of nodes in memory is " << maximumNumberOfNodesInMemory);
remicres's avatar
remicres committed
337

remicres's avatar
remicres committed
338
  // Number of nodes in the entire image
339
340
341
  const std::size_t imageWidth = m_InputImage->GetLargestPossibleRegion().GetSize()[0];
  const std::size_t imageHeight = m_InputImage->GetLargestPossibleRegion().GetSize()[1];
  const std::size_t nbOfNodesInImage = imageWidth*imageHeight;
remicres's avatar
remicres committed
342

343
  // Default layout: 1x1
remicres's avatar
remicres committed
344
345
  m_NbTilesX = 1;
  m_NbTilesY = 1;
346

remicres's avatar
remicres committed
347
348
349
350
351
  // Without margins, the number of tiles maximizing memory use
  // is equal to: nbOfNodesInImage / maximumNumberOfNodesInMemory.
  // Actually, there is tile margins. And the best scenario is to have
  // square tiles with margin = width/2, that is tiles 4x larger.
  // Hence the number of tiles maximizing memory use is 4x larger.
352
  unsigned int minimumNumberOfTiles = std::ceil(4.0 * ((float) nbOfNodesInImage) / ((float) maximumNumberOfNodesInMemory));
remicres's avatar
remicres committed
353
354
355
356
  itkDebugMacro(<<"Minimum number of tiles is " << minimumNumberOfTiles);

  // In the following steps, we will optimize tiling layout, starting from a number
  // of tiles equal to "minimumNumberOfTiles", up to a number of tiles equal to
357
  // 4 times the number of tiles (that is double rows/cols)
remicres's avatar
remicres committed
358
359
360
361
362
363
364
365
  unsigned int maximumNumberOfTiles = minimumNumberOfTiles * 4;

  // Search for layout which minimizes the criterion
  // The criterion is the ratio between compactness and memory usage
  // (i.e. tileWidth * tileHeight / maximumNumberOfNodesInMemory)
  itkDebugMacro(<<"Computing layouts properties:");
  float lowestCriterionValue = itk::NumericTraits<float>::max();
  for (unsigned int nbOfTiles = minimumNumberOfTiles ; nbOfTiles <= maximumNumberOfTiles ; nbOfTiles++)
366
    {
remicres's avatar
remicres committed
367
368
    // Get the multiples of k. For each one, compute the criterion of the tiling
    for (unsigned int layoutNCol = 1; layoutNCol<=nbOfTiles; layoutNCol++)
remicres's avatar
remicres committed
369
      {
remicres's avatar
remicres committed
370
371
#ifdef OTB_USE_MPI
      // We want number of tiles which is a multiple of the number of MPI processes
372
373
      if (nbOfTiles % layoutNCol == 0 && // Is it a multiple of the nb of Tiles and nProcs?
          nbOfTiles % otb::MPIConfig::Instance()->GetNbProcs() == 0)
remicres's avatar
remicres committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
#else
        if (nbOfTiles % layoutNCol == 0) // Is it a multiple of the nb of Tiles?
#endif
          {
          // Tiling layout
          unsigned int layoutNRow = nbOfTiles / layoutNCol;
          unsigned int tileWidth = imageWidth / layoutNCol;
          unsigned int tileHeight = imageHeight / layoutNRow;

          // Compute margin for regular tiles of this layout
          unsigned int maxMargin, maxIter;
          ComputeMaximumStabilityMargin(tileWidth, tileHeight, maxIter, maxMargin);
          tileWidth += 2*maxMargin;
          tileHeight += 2*maxMargin;

          // Memory use efficiency
390
          float percentMemory = tileWidth * tileHeight / (float) maximumNumberOfNodesInMemory; // is > 0. Could be greater than 1 in some cases!
remicres's avatar
remicres committed
391
392
393
394

          // Compactness
          float perimeter = tileWidth + tileHeight;
          float surface = tileWidth * tileHeight;
395
          float compactness = perimeter / surface * (float) std::max(tileWidth,tileHeight); // [1,+inf]
remicres's avatar
remicres committed
396
397
398
399

          // Update minimum criterion
          float criterion = compactness / percentMemory; // ]0, +inf]

400
          itkDebugMacro(//<< std::setprecision (2) << std::fixed
remicres's avatar
remicres committed
401
402
403
404
405
406
407
408
409
              << "Nb. tiles=" << nbOfTiles
              << " Layout: " << layoutNRow << "x" << layoutNCol
              << " Mem. use=" << percentMemory
              << " Compactness=" << compactness
              << " Criterion=" << criterion
              << " Size (no margin): " << (tileWidth-2*maxMargin)<< "x"<< (tileHeight-2*maxMargin)
              << " Size (with margin): " << tileWidth << "x" << tileHeight
              << " (margin=" << maxMargin << "/nb. iter=" << maxIter << ")" );

410
          if (criterion < lowestCriterionValue && percentMemory <= 1.0)
remicres's avatar
remicres committed
411
412
413
414
415
416
417
418
            {
            lowestCriterionValue = criterion;
            m_NbTilesX = layoutNCol;
            m_NbTilesY = layoutNRow;
            }
          }
      } // for each multiple of k
    }
remicres's avatar
remicres committed
419
420
421
422

  // Compute the tile size
  m_TileWidth = static_cast<unsigned int>(imageWidth/m_NbTilesX);
  m_TileHeight = static_cast<unsigned int>(imageHeight/m_NbTilesY);
remicres's avatar
remicres committed
423
424
425
426
427
428
  itkDebugMacro(<<"Selected layout: " << m_NbTilesX << "x" << m_NbTilesY
      << " (criterion=" << lowestCriterionValue << ")");

  // Compute the stability margin
  ComputeMaximumStabilityMargin(m_TileWidth, m_TileHeight,m_NumberOfFirstIterations, m_Margin);

429
430
431
  long long unsigned int memoryUsed = GetNodeMemory();
  memoryUsed *= static_cast<long long unsigned int>(m_TileHeight + 2*m_Margin);
  memoryUsed *= static_cast<long long unsigned int>(m_TileWidth + 2*m_Margin);
432
  itkDebugMacro(<< "An amount of " << memoryUsed/(1024.0*1024.0) << " Mbytes of RAM will be used for regular tiles of size "
remicres's avatar
remicres committed
433
      << (m_TileWidth + 2*m_Margin) << "x" << (m_TileHeight + 2*m_Margin) );
remicres's avatar
remicres committed
434

remicres's avatar
remicres committed
435
436
}

remicres's avatar
remicres committed
437
438
439
template <class TSegmenter>
void Controller<TSegmenter>::SetInternalMemoryAvailable(long long unsigned int v) // expecting a value in Mbytes.
{
440
441
442
443
  if (v<=0)
    {
    itkExceptionMacro(<<"Memory value is not valid (value=" << v << ")");
    }
remicres's avatar
remicres committed
444
445
446
447
  m_Memory = v * 1024ul * 1024ul;
}

template<class TSegmenter>
remicres's avatar
remicres committed
448
void Controller<TSegmenter>::SetInputImage(ImageType * inputImage)
remicres's avatar
remicres committed
449
{
remicres's avatar
remicres committed
450
  m_InputImage = inputImage;
remicres's avatar
remicres committed
451
452
453
454
455
456
457
458
}

template<class TSegmenter>
void Controller<TSegmenter>::SetSpecificParameters(const SegmentationParameterType& params)
{
  m_SpecificParameters = params;
}

459
460
461
462
463
464
465
466
467
468
//template<class TSegmenter>
//typename Controller<TSegmenter>::LabelImageType::Pointer
//Controller<TSegmenter>::GetLabeledClusteredOutput()
//{
//#ifdef OTB_USE_MPI
//  // Get the label image from the master process (the one which achieves segmentation)
//  BroadcastImage<typename TSegmenter::LabelImageType>(m_LabelImage);
//#endif
//  return m_LabelImage;
//}
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483

template <class TSegmenter>
std::vector<std::string> Controller<TSegmenter>::GetTemporaryFilesList()
{
  std::vector<std::string> list;
  for (unsigned int i = 0; i < m_Tiles.size(); i++)
    {
    list.push_back(m_Tiles[i].edgeFileName);
    list.push_back(m_Tiles[i].edgeMarginFileName);
    list.push_back(m_Tiles[i].nodeFileName);
    list.push_back(m_Tiles[i].nodeMarginFileName);
    }
  return list;
}

484
template<class TSegmenter>
485
void Controller<TSegmenter>::SetResumingMode()
486
487
488
489
{
  m_Resuming = true;
}

490
491
492
493
494
495
template<class TSegmenter>
void Controller<TSegmenter>::StopResumingMode()
{
  m_Resuming = false;
}

remicres's avatar
remicres committed
496
} // end of namespace lsgrm
497
#endif