lsgrmController.txx 14.9 KB
Newer Older
remicres's avatar
remicres committed
1
#include "lsgrmController.h"
2
3
4
5

namespace lsgrm
{

remicres's avatar
remicres committed
6
7
8
template<class TSegmenter>
Controller<TSegmenter>::Controller()
{
9
  m_TilingMode = LSGRM_TILING_AUTO;
remicres's avatar
remicres committed
10
11
12
13
14
15
16
17
18
  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
19
20
21
22
23
24
25
}

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

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

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

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

    std::cout <<
remicres's avatar
remicres committed
56
57
58
59
60
61
62
        "--- 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
63
64

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

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

  if (m_TilingMode == LSGRM_TILING_AUTO || m_TilingMode == LSGRM_TILING_USER)
    {
78
    unsigned int numberOfIterationsForPartialSegmentations = 3; // TODO: find a smart value
79
80
    unsigned int numberOfIterationsRemaining = m_NumberOfIterations;

remicres's avatar
remicres committed
81
82
83
84
85
86
87
88
89
90
    // 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,
91
        numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
92
93
94
95
96
97
98
        m_Tiles,
        m_NbTilesX,
        m_NbTilesY,
        m_TileWidth,
        m_TileHeight,
        isFusion);

99
    // Update the given number of iterations
remicres's avatar
remicres committed
100
    numberOfIterationsRemaining -= m_NumberOfFirstIterations;
101

102
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
103
104
    // Gathering useful variables
    GatherUsefulVariables(accumulatedMemory, isFusion);
105
#endif
remicres's avatar
remicres committed
106
107
108
109

    // Time monitoring
    ShowTime(t);

remicres's avatar
remicres committed
110
111
    std::cout << "accumulatedMemory=" << accumulatedMemory << std::endl;

remicres's avatar
remicres committed
112
    while(accumulatedMemory > m_Memory && isFusion)
remicres's avatar
remicres committed
113
      {
114
115
116
117
118
119
120
121
122
123
      if (numberOfIterationsRemaining < numberOfIterationsForPartialSegmentations)
        {
        numberOfIterationsForPartialSegmentations = numberOfIterationsRemaining;
        numberOfIterationsRemaining = 0;
        }
      else
        {
        numberOfIterationsRemaining -= numberOfIterationsForPartialSegmentations;
        }

remicres's avatar
remicres committed
124
125
      isFusion = false;
      accumulatedMemory = RunPartialSegmentation<TSegmenter>(
remicres's avatar
remicres committed
126
127
          m_SpecificParameters,
          m_Threshold,
128
          numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
129
130
131
132
133
134
          m_Tiles,
          m_NbTilesX,
          m_NbTilesY,
          m_InputImage->GetLargestPossibleRegion().GetSize()[0],
          m_InputImage->GetLargestPossibleRegion().GetSize()[1],
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
135
136
          isFusion);

137
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
138
139
      // Gathering useful variables
      GatherUsefulVariables(accumulatedMemory, isFusion);
140
#endif
141

remicres's avatar
remicres committed
142
143
      // Time monitoring
      ShowTime(t);
144

145
146
      if (numberOfIterationsRemaining == 0)
        break;
remicres's avatar
remicres committed
147
148
149
      }

#ifdef OTB_USE_MPI
150
    // Only the master process is doing the next part
remicres's avatar
remicres committed
151
    // TODO: Use the MPI process wich has the largest amount of memory
remicres's avatar
remicres committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
    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],
167
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
168
          numberOfIterationsRemaining);
remicres's avatar
remicres committed
169

170
      ShowTime(t);
remicres's avatar
remicres committed
171
172

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

remicres's avatar
remicres committed
195
196
197
198
199
200
201
202
    // Use classic grm
    TSegmenter segmenter;
    segmenter.SetParam(m_SpecificParameters);
    segmenter.SetThreshold(m_Threshold);
    segmenter.SetDoFastSegmentation(true);
    segmenter.SetNumberOfIterations(m_NumberOfIterations);
    segmenter.SetInput(m_InputImage);
    segmenter.Update();
203
204
205

    // Get label image
    m_LabelImage = segmenter.GetLabeledClusteredOutput();
remicres's avatar
remicres committed
206
207
208
209
    }
  else
    {
    itkExceptionMacro(<<"Unknow tiling mode!");
remicres's avatar
remicres committed
210
211
212
213
    }

}

214
215
/*
 * Compute the memory occupied by one node
216
 * TODO: compute the exact value
217
218
219
220
 */
template<class TSegmenter>
unsigned int Controller<TSegmenter>::GetNodeMemory()
{
221
  // Create a n*n image
222
  const unsigned int n = 100;
223
224
225
226
  typename ImageType::Pointer onePixelImage = ImageType::New();
  typename ImageType::IndexType start;
  start.Fill(0);
  typename ImageType::SizeType size;
227
  size.Fill(n);
228
229
230
231
  typename ImageType::RegionType region(start, size);
  onePixelImage->SetRegions(region);
  onePixelImage->SetNumberOfComponentsPerPixel(m_InputImage->GetNumberOfComponentsPerPixel());
  onePixelImage->Allocate();
232
233

  // Instanciate and initialize a segmenter
234
235
  TSegmenter segmenter;
  segmenter.SetInput(onePixelImage);
236
  grm::GraphOperations<TSegmenter>::InitNodes(onePixelImage,segmenter,FOUR);
237

238
  // Get the memory occupied by the graph, normalize it by n*n
239
  unsigned int memory = segmenter.GetGraphMemory() / (n*n);
240

241
  itkDebugMacro(<<"Size of a node is " << memory);
242

243
  return memory;
244
}
remicres's avatar
remicres committed
245

246
247
248
/*
 * Compute the maximum number of nodes which can fit in the system memory
 */
remicres's avatar
remicres committed
249
template<class TSegmenter>
remicres's avatar
remicres committed
250
long unsigned int Controller<TSegmenter>::GetMaximumNumberOfNodesInMemory()
remicres's avatar
remicres committed
251
{
remicres's avatar
remicres committed
252
253
  itkDebugMacro(<< "Computing maximum number of nodes in memory");

remicres's avatar
remicres committed
254
255
256
257
258
  m_Memory = getMemorySize();
  assert(m_Memory > 0);

  m_Memory /= 2; // For safety and can prevent out of memory troubles

259
  return std::ceil(((float) m_Memory) / ((float) GetNodeMemory()));
remicres's avatar
remicres committed
260
261
262
263

}

template<class TSegmenter>
264
void Controller<TSegmenter>::ComputeMaximumStabilityMargin(unsigned int width,
remicres's avatar
remicres committed
265
    unsigned int height, unsigned int &niter, unsigned int &margin)
266
    {
remicres's avatar
remicres committed
267
268
269
270
  itkDebugMacro(<< "Computing maximum stability margin");

  // Compute the stability margin. The naive strategy consider a margin value and a stable size equal.
  niter = 1;
271
272
273
  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
274

275
  while(currMargin < maxMargin)
remicres's avatar
remicres committed
276
    {
277
    margin = currMargin;
remicres's avatar
remicres committed
278
    niter++;
279
    currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
remicres's avatar
remicres committed
280
281
282
283
    }
  niter--;

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

285
    }
remicres's avatar
remicres committed
286

287
/*
288
 * Compute a tiling layout which minimizes a criterion based on tile compactness
remicres's avatar
remicres committed
289
 * and memory usage
290
291
 *
 * TODO: use the lsgrmSplitter to truly compute the largest tile of a given layout
292
 */
remicres's avatar
remicres committed
293
294
295
296
template<class TSegmenter>
void Controller<TSegmenter>::GetAutomaticConfiguration()
{

remicres's avatar
remicres committed
297
298
  itkDebugMacro(<<"Get automatic configuration");

remicres's avatar
remicres committed
299
  // Compute the maximum number of nodes that can fit the memory
remicres's avatar
remicres committed
300
  // TODO: Use the smallest number amongst MPI processes
remicres's avatar
remicres committed
301
  unsigned long int maximumNumberOfNodesInMemory = GetMaximumNumberOfNodesInMemory();
remicres's avatar
remicres committed
302
  itkDebugMacro(<<"Maximum number of nodes in memory is " << maximumNumberOfNodesInMemory);
remicres's avatar
remicres committed
303

remicres's avatar
remicres committed
304
  // Number of nodes in the entire image
remicres's avatar
remicres committed
305
306
  const unsigned int imageWidth = m_InputImage->GetLargestPossibleRegion().GetSize()[0];
  const unsigned int imageHeight = m_InputImage->GetLargestPossibleRegion().GetSize()[1];
remicres's avatar
remicres committed
307
308
  const unsigned long int nbOfNodesInImage = imageWidth*imageHeight;

309
  // Default layout: 1x1
remicres's avatar
remicres committed
310
311
  m_NbTilesX = 1;
  m_NbTilesY = 1;
312

remicres's avatar
remicres committed
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
  // 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.
  unsigned int minimumNumberOfTiles = std::ceil(4 * nbOfNodesInImage / ((float) maximumNumberOfNodesInMemory));
  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
  // twice the number of tiles (that is memory usage about 50%)
  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++)
332
    {
remicres's avatar
remicres committed
333
334
    // 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
335
      {
remicres's avatar
remicres committed
336
337
#ifdef OTB_USE_MPI
      // We want number of tiles which is a multiple of the number of MPI processes
338
339
      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
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
#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
          float percentMemory = tileWidth * tileHeight / (float) maximumNumberOfNodesInMemory; // ]0, 1]

          // Compactness
          float perimeter = tileWidth + tileHeight;
          float surface = tileWidth * tileHeight;
361
          float compactness = perimeter / surface * (float) std::max(tileWidth,tileHeight); // [1,+inf]
remicres's avatar
remicres committed
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384

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

          itkDebugMacro(<< std::setprecision (2) << std::fixed
              << "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 << ")" );

          if (criterion < lowestCriterionValue)
            {
            lowestCriterionValue = criterion;
            m_NbTilesX = layoutNCol;
            m_NbTilesY = layoutNRow;
            }
          }
      } // for each multiple of k
    }
remicres's avatar
remicres committed
385
386
387
388

  // 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
389
390
391
392
393
394
  itkDebugMacro(<<"Selected layout: " << m_NbTilesX << "x" << m_NbTilesY
      << " (criterion=" << lowestCriterionValue << ")");

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

395
396
397
  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);
398
  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
399
      << (m_TileWidth + 2*m_Margin) << "x" << (m_TileHeight + 2*m_Margin) );
remicres's avatar
remicres committed
400

remicres's avatar
remicres committed
401
402
}

remicres's avatar
remicres committed
403
404
405
406
407
408
409
410
template <class TSegmenter>
void Controller<TSegmenter>::SetInternalMemoryAvailable(long long unsigned int v) // expecting a value in Mbytes.
{
  assert(v > 0);
  m_Memory = v * 1024ul * 1024ul;
}

template<class TSegmenter>
remicres's avatar
remicres committed
411
void Controller<TSegmenter>::SetInputImage(ImageType * inputImage)
remicres's avatar
remicres committed
412
{
remicres's avatar
remicres committed
413
  m_InputImage = inputImage;
remicres's avatar
remicres committed
414
415
416
417
418
419
420
421
}

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

remicres's avatar
remicres committed
422
423
424
425
template<class TSegmenter>
typename Controller<TSegmenter>::LabelImageType::Pointer
Controller<TSegmenter>::GetLabeledClusteredOutput()
{
remicres's avatar
remicres committed
426
427
428
429
#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
430
431
  return m_LabelImage;
}
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446

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
447
} // end of namespace lsgrm