lsgrmController.txx 15.4 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;
remicres's avatar
remicres committed
21
22
23
24
25
26
27
}

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

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

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

44
45
  CheckMemorySize();

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

    std::cout <<
remicres's avatar
remicres committed
60
61
62
63
64
65
66
        "--- 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
67
68

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

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

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

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

    // Run first partial segmentation
    boost::timer t; t.restart();
    auto accumulatedMemory = RunFirstPartialSegmentation<TSegmenter>(
        m_InputImage,
        m_SpecificParameters,
        m_Threshold,
        m_NumberOfFirstIterations,
95
        numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
96
97
98
99
100
101
102
        m_Tiles,
        m_NbTilesX,
        m_NbTilesY,
        m_TileWidth,
        m_TileHeight,
        isFusion);

103
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
104
    GatherUsefulVariables(accumulatedMemory, isFusion);
105
#endif
remicres's avatar
remicres committed
106
107
108
109
110

    // Time monitoring
    ShowTime(t);

    while(accumulatedMemory > m_Memory && isFusion)
remicres's avatar
remicres committed
111
      {
remicres's avatar
remicres committed
112
113
      isFusion = false;
      accumulatedMemory = RunPartialSegmentation<TSegmenter>(
remicres's avatar
remicres committed
114
115
          m_SpecificParameters,
          m_Threshold,
116
          numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
117
118
119
120
121
122
          m_Tiles,
          m_NbTilesX,
          m_NbTilesY,
          m_InputImage->GetLargestPossibleRegion().GetSize()[0],
          m_InputImage->GetLargestPossibleRegion().GetSize()[1],
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
123
124
          isFusion);

125
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
126
      GatherUsefulVariables(accumulatedMemory, isFusion);
127
#endif
128

remicres's avatar
remicres committed
129
130
      // Time monitoring
      ShowTime(t);
131

132
133
134
      // Update number of remaining iterations
      if (numberOfIterationsRemaining < numberOfIterationsForPartialSegmentations)
        {
135
        break;
136
137
138
139
140
        }
      else
        {
        numberOfIterationsRemaining -= numberOfIterationsForPartialSegmentations;
        }
remicres's avatar
remicres committed
141
142
143
      }

#ifdef OTB_USE_MPI
144
    // Only the master process is doing the next part
remicres's avatar
remicres committed
145
    // TODO: Use the MPI process wich has the largest amount of memory
remicres's avatar
remicres committed
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
    if (otb::MPIConfig::Instance()->GetMyRank() != 0)
      return;
#endif

    if(accumulatedMemory <= m_Memory)
      {
      // Merge all the graphs
      m_LabelImage = MergeAllGraphsAndAchieveSegmentation<TSegmenter>(
          m_SpecificParameters,
          m_Threshold,
          m_Tiles,
          m_NbTilesX,
          m_NbTilesY,
          m_InputImage->GetLargestPossibleRegion().GetSize()[0],
          m_InputImage->GetLargestPossibleRegion().GetSize()[1],
161
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
162
          numberOfIterationsRemaining);
remicres's avatar
remicres committed
163

164
      ShowTime(t);
remicres's avatar
remicres committed
165
166

      }
167
    else // accumulatedMemory > m_Memory
remicres's avatar
remicres committed
168
      {
remicres's avatar
remicres committed
169
      // That means there are no more possible fusions but we can not store the output graph
remicres's avatar
remicres committed
170
171
      // 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.
172
      itkExceptionMacro(<< "No more possible fusions, but can not store the output graph");
remicres's avatar
remicres committed
173
174
      }
    }
remicres's avatar
remicres committed
175
  else if (m_TilingMode == LSGRM_TILING_NONE)// tiling_mode is none
remicres's avatar
remicres committed
176
    {
remicres's avatar
remicres committed
177
178
179
180
181
182
183
184
185
#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
186
187
188
    // Update input image
    m_InputImage->Update();

remicres's avatar
remicres committed
189
190
191
192
    // Use classic grm
    TSegmenter segmenter;
    segmenter.SetParam(m_SpecificParameters);
    segmenter.SetThreshold(m_Threshold);
193
    segmenter.SetDoFastSegmentation(false);
remicres's avatar
remicres committed
194
195
196
    segmenter.SetNumberOfIterations(m_NumberOfIterations);
    segmenter.SetInput(m_InputImage);
    segmenter.Update();
197
198
199

    // Get label image
    m_LabelImage = segmenter.GetLabeledClusteredOutput();
remicres's avatar
remicres committed
200
201
202
203
    }
  else
    {
    itkExceptionMacro(<<"Unknow tiling mode!");
remicres's avatar
remicres committed
204
205
206
207
    }

}

208
209
/*
 * Compute the memory occupied by one node
remicres's avatar
remicres committed
210
211
 * TODO: compute the exact value, e.g. on a given UNIX system,
 * experimental measures shows that
212
 * for one Baatz node (+pixel) memory is about 700-730 bytes...
remicres's avatar
remicres committed
213
 * And our estimation is about 600
214
215
216
217
 */
template<class TSegmenter>
unsigned int Controller<TSegmenter>::GetNodeMemory()
{
218
  // Create a n*n image
219
  const unsigned int n = 100;
220
221
222
223
  typename ImageType::Pointer onePixelImage = ImageType::New();
  typename ImageType::IndexType start;
  start.Fill(0);
  typename ImageType::SizeType size;
224
  size.Fill(n);
225
226
227
228
  typename ImageType::RegionType region(start, size);
  onePixelImage->SetRegions(region);
  onePixelImage->SetNumberOfComponentsPerPixel(m_InputImage->GetNumberOfComponentsPerPixel());
  onePixelImage->Allocate();
229
230

  // Instanciate and initialize a segmenter
231
232
  TSegmenter segmenter;
  segmenter.SetInput(onePixelImage);
233
  grm::GraphOperations<TSegmenter>::InitNodes(onePixelImage,segmenter,FOUR);
234

235
  // Get the memory occupied by the graph, normalize it by n*n
236
  unsigned int memory = segmenter.GetGraphMemory() / (n*n);
237

238
  itkDebugMacro(<<"Size of a node is " << memory);
239
240

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

244
  itkDebugMacro(<<"Size of an image pixel is " << pixelMemory);
245
246
247

  memory += pixelMemory;

248
  itkDebugMacro(<<"Size of a node+pixel is " << memory);
249

250
  return memory;
251
}
remicres's avatar
remicres committed
252

remicres's avatar
remicres committed
253
template<class TSegmenter>
254
void Controller<TSegmenter>::CheckMemorySize()
remicres's avatar
remicres committed
255
{
256
257
258
259
260
  if (m_Memory == 0)
    {
    m_Memory = getMemorySize();
    assert(m_Memory > 0);
    }
remicres's avatar
remicres committed
261
  m_Memory /= 2; // For safety and can prevent out of memory troubles
262
263
264
265
266
267
268
269
}
/*
 * Compute the maximum number of nodes which can fit in the memory
 */
template<class TSegmenter>
long unsigned int Controller<TSegmenter>::GetMaximumNumberOfNodesInMemory()
{
  itkDebugMacro(<< "Computing maximum number of nodes in memory");
remicres's avatar
remicres committed
270

271
  return std::ceil(((float) m_Memory) / ((float) GetNodeMemory()));
remicres's avatar
remicres committed
272
273
274
}

template<class TSegmenter>
275
void Controller<TSegmenter>::ComputeMaximumStabilityMargin(unsigned int width,
remicres's avatar
remicres committed
276
    unsigned int height, unsigned int &niter, unsigned int &margin)
277
    {
remicres's avatar
remicres committed
278
279
280
281
  itkDebugMacro(<< "Computing maximum stability margin");

  // Compute the stability margin. The naive strategy consider a margin value and a stable size equal.
  niter = 1;
282
283
284
  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
285

286
  while(currMargin < maxMargin)
remicres's avatar
remicres committed
287
    {
288
    margin = currMargin;
remicres's avatar
remicres committed
289
    niter++;
290
    currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
remicres's avatar
remicres committed
291
292
293
294
    }
  niter--;

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

296
    }
remicres's avatar
remicres committed
297

298
/*
299
 * Compute a tiling layout which minimizes a criterion based on tile compactness
remicres's avatar
remicres committed
300
 * and memory usage
301
302
 *
 * TODO: use the lsgrmSplitter to truly compute the largest tile of a given layout
303
 */
remicres's avatar
remicres committed
304
305
306
307
template<class TSegmenter>
void Controller<TSegmenter>::GetAutomaticConfiguration()
{

remicres's avatar
remicres committed
308
309
  itkDebugMacro(<<"Get automatic configuration");

remicres's avatar
remicres committed
310
  // Compute the maximum number of nodes that can fit the memory
remicres's avatar
remicres committed
311
  // TODO: Use the smallest number amongst MPI processes
remicres's avatar
remicres committed
312
  unsigned long int maximumNumberOfNodesInMemory = GetMaximumNumberOfNodesInMemory();
remicres's avatar
remicres committed
313
  itkDebugMacro(<<"Maximum number of nodes in memory is " << maximumNumberOfNodesInMemory);
remicres's avatar
remicres committed
314

remicres's avatar
remicres committed
315
  // Number of nodes in the entire image
remicres's avatar
remicres committed
316
317
  const unsigned int imageWidth = m_InputImage->GetLargestPossibleRegion().GetSize()[0];
  const unsigned int imageHeight = m_InputImage->GetLargestPossibleRegion().GetSize()[1];
remicres's avatar
remicres committed
318
319
  const unsigned long int nbOfNodesInImage = imageWidth*imageHeight;

320
  // Default layout: 1x1
remicres's avatar
remicres committed
321
322
  m_NbTilesX = 1;
  m_NbTilesY = 1;
323

remicres's avatar
remicres committed
324
325
326
327
328
  // 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.
329
  unsigned int minimumNumberOfTiles = std::ceil(4.0 * ((float) nbOfNodesInImage) / ((float) maximumNumberOfNodesInMemory));
remicres's avatar
remicres committed
330
331
332
333
  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
334
  // 4 times the number of tiles (that is double rows/cols)
remicres's avatar
remicres committed
335
336
337
338
339
340
341
342
  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++)
343
    {
remicres's avatar
remicres committed
344
345
    // 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
346
      {
remicres's avatar
remicres committed
347
348
#ifdef OTB_USE_MPI
      // We want number of tiles which is a multiple of the number of MPI processes
349
350
      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
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
#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
367
          float percentMemory = tileWidth * tileHeight / (float) maximumNumberOfNodesInMemory; // is > 0. Could be greater than 1 in some cases!
remicres's avatar
remicres committed
368
369
370
371

          // Compactness
          float perimeter = tileWidth + tileHeight;
          float surface = tileWidth * tileHeight;
372
          float compactness = perimeter / surface * (float) std::max(tileWidth,tileHeight); // [1,+inf]
remicres's avatar
remicres committed
373
374
375
376

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

377
          itkDebugMacro(//<< std::setprecision (2) << std::fixed
remicres's avatar
remicres committed
378
379
380
381
382
383
384
385
386
              << "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 << ")" );

387
          if (criterion < lowestCriterionValue && percentMemory <= 1.0)
remicres's avatar
remicres committed
388
389
390
391
392
393
394
395
            {
            lowestCriterionValue = criterion;
            m_NbTilesX = layoutNCol;
            m_NbTilesY = layoutNRow;
            }
          }
      } // for each multiple of k
    }
remicres's avatar
remicres committed
396
397
398
399

  // 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
400
401
402
403
404
405
  itkDebugMacro(<<"Selected layout: " << m_NbTilesX << "x" << m_NbTilesY
      << " (criterion=" << lowestCriterionValue << ")");

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

406
407
408
  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);
409
  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
410
      << (m_TileWidth + 2*m_Margin) << "x" << (m_TileHeight + 2*m_Margin) );
remicres's avatar
remicres committed
411

remicres's avatar
remicres committed
412
413
}

remicres's avatar
remicres committed
414
415
416
template <class TSegmenter>
void Controller<TSegmenter>::SetInternalMemoryAvailable(long long unsigned int v) // expecting a value in Mbytes.
{
417
418
419
420
  if (v<=0)
    {
    itkExceptionMacro(<<"Memory value is not valid (value=" << v << ")");
    }
remicres's avatar
remicres committed
421
422
423
424
  m_Memory = v * 1024ul * 1024ul;
}

template<class TSegmenter>
remicres's avatar
remicres committed
425
void Controller<TSegmenter>::SetInputImage(ImageType * inputImage)
remicres's avatar
remicres committed
426
{
remicres's avatar
remicres committed
427
  m_InputImage = inputImage;
remicres's avatar
remicres committed
428
429
430
431
432
433
434
435
}

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

remicres's avatar
remicres committed
436
437
438
439
template<class TSegmenter>
typename Controller<TSegmenter>::LabelImageType::Pointer
Controller<TSegmenter>::GetLabeledClusteredOutput()
{
remicres's avatar
remicres committed
440
441
442
443
#ifdef OTB_USE_MPI
  // Get the label image from the master process (the one which achieves segmentation)
  BroadcastImage<typename TSegmenter::LabelImageType>(m_LabelImage);
#endif
remicres's avatar
remicres committed
444
445
  return m_LabelImage;
}
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460

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

remicres's avatar
remicres committed
461
} // end of namespace lsgrm
462
#endif