lsgrmController.txx 14 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();
}

remicres's avatar
remicres committed
33
34
35
36
template<class TSegmenter>
void Controller<TSegmenter>::RunSegmentation()
{

37
  if (m_TilingMode == LSGRM_TILING_AUTO || m_TilingMode == LSGRM_TILING_USER)
remicres's avatar
remicres committed
38
    {
remicres's avatar
remicres committed
39

40
41
42
43
44
45
46
47
48
49
    if(m_TilingMode == LSGRM_TILING_AUTO)
      {
      this->GetAutomaticConfiguration();
      }
    else if (m_TilingMode == LSGRM_TILING_USER)
      {
      m_Margin = static_cast<unsigned int>(pow(2, m_NumberOfFirstIterations + 1) - 2);
      }

    std::cout <<
remicres's avatar
remicres committed
50
51
52
53
54
55
56
        "--- 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
57
58

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

62
    // If there is only one tile, then fallback to LSGRM_TILING_NONE case
63
64
65
66
67
68
    if (m_Tiles.size() == 1)
      {
      std::cout << "Only one tile is needed. Fallback to tiling=none." << std::endl;
      SetTilingModeNone();
      }

69
70
71
72
73
74
75
76
    }

  if (m_TilingMode == LSGRM_TILING_AUTO || m_TilingMode == LSGRM_TILING_USER)
    {

    const unsigned int numberOfIterationsForPartialSegmentations = 3; // TODO: find a smart value
    unsigned int numberOfIterationsRemaining = m_NumberOfIterations;

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

95
    // Update the given number of iterations
remicres's avatar
remicres committed
96
    numberOfIterationsRemaining -= m_NumberOfFirstIterations;
97

98
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
99
100
    // Gathering useful variables
    GatherUsefulVariables(accumulatedMemory, isFusion);
101
#endif
remicres's avatar
remicres committed
102
103
104
105

    // Time monitoring
    ShowTime(t);

remicres's avatar
remicres committed
106
107
    std::cout << "accumulatedMemory=" << accumulatedMemory << std::endl;

remicres's avatar
remicres committed
108
    while(accumulatedMemory > m_Memory && isFusion)
remicres's avatar
remicres committed
109
      {
remicres's avatar
remicres committed
110
111
112

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

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

remicres's avatar
remicres committed
129
130
      // Time monitoring
      ShowTime(t);
131
132
133
134
135
136
137
138
139
140

      // Update the number of iterations remaining
      if (numberOfIterationsRemaining > numberOfIterationsForPartialSegmentations)
        {
        numberOfIterationsRemaining -= numberOfIterationsForPartialSegmentations;
        }
      else
        {
        break;
        }
remicres's avatar
remicres committed
141
142
143
144
      }


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

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

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

remicres's avatar
remicres committed
181
182
183
184
185
186
187
188
    // 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();
189
190
191

    // Get label image
    m_LabelImage = segmenter.GetLabeledClusteredOutput();
remicres's avatar
remicres committed
192
193
194
195
    }
  else
    {
    itkExceptionMacro(<<"Unknow tiling mode!");
remicres's avatar
remicres committed
196
197
198
199
    }

}

200
201
/*
 * Compute the memory occupied by one node
202
 * TODO: compute the exact value
203
204
205
206
 */
template<class TSegmenter>
unsigned int Controller<TSegmenter>::GetNodeMemory()
{
207
  // Create a n*n image
208
  const unsigned int n = 100;
209
210
211
212
  typename ImageType::Pointer onePixelImage = ImageType::New();
  typename ImageType::IndexType start;
  start.Fill(0);
  typename ImageType::SizeType size;
213
  size.Fill(n);
214
215
216
217
  typename ImageType::RegionType region(start, size);
  onePixelImage->SetRegions(region);
  onePixelImage->SetNumberOfComponentsPerPixel(m_InputImage->GetNumberOfComponentsPerPixel());
  onePixelImage->Allocate();
218
219

  // Instanciate and initialize a segmenter
220
221
  TSegmenter segmenter;
  segmenter.SetInput(onePixelImage);
222
  grm::GraphOperations<TSegmenter>::InitNodes(onePixelImage,segmenter,FOUR);
223

224
  // Get the memory occupied by the graph, normalize it by n*n
225
  unsigned int memory = segmenter.GetGraphMemory() / (n*n);
226

227
  itkDebugMacro(<<"Size of a node is " << memory);
228

229
  return memory;
230
}
remicres's avatar
remicres committed
231

232
233
234
/*
 * Compute the maximum number of nodes which can fit in the system memory
 */
remicres's avatar
remicres committed
235
template<class TSegmenter>
remicres's avatar
remicres committed
236
long unsigned int Controller<TSegmenter>::GetMaximumNumberOfNodesInMemory()
remicres's avatar
remicres committed
237
{
remicres's avatar
remicres committed
238
239
  itkDebugMacro(<< "Computing maximum number of nodes in memory");

remicres's avatar
remicres committed
240
241
242
243
244
  m_Memory = getMemorySize();
  assert(m_Memory > 0);

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

245
  return std::ceil(((float) m_Memory) / ((float) GetNodeMemory()));
remicres's avatar
remicres committed
246
247
248
249

}

template<class TSegmenter>
250
void Controller<TSegmenter>::ComputeMaximumStabilityMargin(unsigned int width,
remicres's avatar
remicres committed
251
    unsigned int height, unsigned int &niter, unsigned int &margin)
252
    {
remicres's avatar
remicres committed
253
254
255
256
  itkDebugMacro(<< "Computing maximum stability margin");

  // Compute the stability margin. The naive strategy consider a margin value and a stable size equal.
  niter = 1;
257
258
259
  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
260

261
  while(currMargin < maxMargin)
remicres's avatar
remicres committed
262
    {
263
    margin = currMargin;
remicres's avatar
remicres committed
264
    niter++;
265
    currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
remicres's avatar
remicres committed
266
267
268
269
    }
  niter--;

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

271
    }
remicres's avatar
remicres committed
272

273
/*
274
 * Compute a tiling layout which minimizes a criterion based on tile compactness
remicres's avatar
remicres committed
275
 * and memory usage
276
277
 *
 * TODO: use the lsgrmSplitter to truly compute the largest tile of a given layout
278
 */
remicres's avatar
remicres committed
279
280
281
282
template<class TSegmenter>
void Controller<TSegmenter>::GetAutomaticConfiguration()
{

remicres's avatar
remicres committed
283
284
  itkDebugMacro(<<"Get automatic configuration");

remicres's avatar
remicres committed
285
  // Compute the maximum number of nodes that can fit the memory
remicres's avatar
remicres committed
286
  // TODO: Use the smallest number amongst MPI processes
remicres's avatar
remicres committed
287
  unsigned long int maximumNumberOfNodesInMemory = GetMaximumNumberOfNodesInMemory();
remicres's avatar
remicres committed
288
  itkDebugMacro(<<"Maximum number of nodes in memory is " << maximumNumberOfNodesInMemory);
remicres's avatar
remicres committed
289

remicres's avatar
remicres committed
290
  // Number of nodes in the entire image
remicres's avatar
remicres committed
291
292
  const unsigned int imageWidth = m_InputImage->GetLargestPossibleRegion().GetSize()[0];
  const unsigned int imageHeight = m_InputImage->GetLargestPossibleRegion().GetSize()[1];
remicres's avatar
remicres committed
293
294
  const unsigned long int nbOfNodesInImage = imageWidth*imageHeight;

295
  // Default layout: 1x1
remicres's avatar
remicres committed
296
297
  m_NbTilesX = 1;
  m_NbTilesY = 1;
298

remicres's avatar
remicres committed
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
  // 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++)
318
    {
remicres's avatar
remicres committed
319
320
    // 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
321
      {
remicres's avatar
remicres committed
322
323
#ifdef OTB_USE_MPI
      // We want number of tiles which is a multiple of the number of MPI processes
324
325
      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
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
#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;
347
          float compactness = perimeter / surface * (float) std::max(tileWidth,tileHeight); // [1,+inf]
remicres's avatar
remicres committed
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370

          // 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
371
372
373
374

  // 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
375
376
377
378
379
380
  itkDebugMacro(<<"Selected layout: " << m_NbTilesX << "x" << m_NbTilesY
      << " (criterion=" << lowestCriterionValue << ")");

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

381
382
383
  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);
384
  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
385
      << (m_TileWidth + 2*m_Margin) << "x" << (m_TileHeight + 2*m_Margin) );
remicres's avatar
remicres committed
386

remicres's avatar
remicres committed
387
388
}

remicres's avatar
remicres committed
389
390
391
392
393
394
395
396
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
397
void Controller<TSegmenter>::SetInputImage(ImageType * inputImage)
remicres's avatar
remicres committed
398
{
remicres's avatar
remicres committed
399
  m_InputImage = inputImage;
remicres's avatar
remicres committed
400
401
402
403
404
405
406
407
}

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

remicres's avatar
remicres committed
408
409
410
411
412
413
template<class TSegmenter>
typename Controller<TSegmenter>::LabelImageType::Pointer
Controller<TSegmenter>::GetLabeledClusteredOutput()
{
  return m_LabelImage;
}
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428

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