lsgrmController.txx 13.4 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
  m_Margin = 0;
  m_NumberOfIterations = 0;
  m_NumberOfFirstIterations = 0;
  m_TileHeight = 0;
  m_TileWidth = 0;
  m_NbTilesX = 0;
  m_NbTilesY = 0;
17
  m_CleanTemporaryFiles = true;
remicres's avatar
remicres committed
18
19
  m_Threshold = 75;
  m_Memory = 0;
remicres's avatar
remicres committed
20
21
22
23
24
25
26
27
28
29
30
}

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

template<class TSegmenter>
void Controller<TSegmenter>::RunSegmentation()
{

31
  if (m_TilingMode == LSGRM_TILING_AUTO || m_TilingMode == LSGRM_TILING_USER)
remicres's avatar
remicres committed
32
    {
remicres's avatar
remicres committed
33
34
35

    const unsigned int numberOfIterationsForPartialSegmentations = 3; // TODO: find a smart value
    unsigned int numberOfIterationsRemaining = m_NumberOfIterations;
36
37
38
39
40
41
42
43
44
45
46

    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
47
48
49
50
51
52
53
        "--- 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
54
55

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

remicres's avatar
remicres committed
59
60
61
62
63
64
65
66
67
68
    // 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,
69
        numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
70
71
72
73
74
75
76
        m_Tiles,
        m_NbTilesX,
        m_NbTilesY,
        m_TileWidth,
        m_TileHeight,
        isFusion);

77
    // Update the given number of iterations
remicres's avatar
remicres committed
78
    numberOfIterationsRemaining -= m_NumberOfFirstIterations;
79

80
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
81
82
    // Gathering useful variables
    GatherUsefulVariables(accumulatedMemory, isFusion);
83
#endif
remicres's avatar
remicres committed
84
85
86
87

    // Time monitoring
    ShowTime(t);

remicres's avatar
remicres committed
88
89
    std::cout << "accumulatedMemory=" << accumulatedMemory << std::endl;

remicres's avatar
remicres committed
90
    while(accumulatedMemory > m_Memory && isFusion)
remicres's avatar
remicres committed
91
      {
remicres's avatar
remicres committed
92
93
94

      isFusion = false;
      accumulatedMemory = RunPartialSegmentation<TSegmenter>(
remicres's avatar
remicres committed
95
96
          m_SpecificParameters,
          m_Threshold,
97
          numberOfIterationsForPartialSegmentations,
remicres's avatar
remicres committed
98
99
100
101
102
103
          m_Tiles,
          m_NbTilesX,
          m_NbTilesY,
          m_InputImage->GetLargestPossibleRegion().GetSize()[0],
          m_InputImage->GetLargestPossibleRegion().GetSize()[1],
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
104
105
          isFusion);

106
107
108
109
110
      // Update the number of iterations remaining
      if (numberOfIterationsRemaining > numberOfIterationsForPartialSegmentations)
        {
        numberOfIterationsRemaining -= numberOfIterationsForPartialSegmentations;
        }
111

112
#ifdef OTB_USE_MPI
remicres's avatar
remicres committed
113
114
      // Gathering useful variables
      GatherUsefulVariables(accumulatedMemory, isFusion);
115
#endif
remicres's avatar
remicres committed
116
117
118
119
120
121
      // Time monitoring
      ShowTime(t);
      }


#ifdef OTB_USE_MPI
122
    // Only the master process is doing the next part
remicres's avatar
remicres committed
123
    // TODO: Use the MPI process wich has the largest amount of memory
remicres's avatar
remicres committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
    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],
139
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
140
          numberOfIterationsRemaining);
remicres's avatar
remicres committed
141

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

      }
145
    else // accumulatedMemory > m_Memory
remicres's avatar
remicres committed
146
      {
remicres's avatar
remicres committed
147
      // That means there are no more possible fusions but we can not store the output graph
remicres's avatar
remicres committed
148
149
      // 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.
150
      itkExceptionMacro(<< "No more possible fusions, but can not store the output graph");
remicres's avatar
remicres committed
151
152
      }
    }
remicres's avatar
remicres committed
153
  else if (m_TilingMode == LSGRM_TILING_NONE)// tiling_mode is none
remicres's avatar
remicres committed
154
    {
remicres's avatar
remicres committed
155
156
157
158
159
160
161
162
163
164
165
166
    // 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();
    }
  else
    {
    itkExceptionMacro(<<"Unknow tiling mode!");
remicres's avatar
remicres committed
167
168
169
170
    }

}

171
172
173
174
175
176
/*
 * Compute the memory occupied by one node
 */
template<class TSegmenter>
unsigned int Controller<TSegmenter>::GetNodeMemory()
{
177
  // Create a n*n image
178
  const unsigned int n = 100;
179
180
181
182
  typename ImageType::Pointer onePixelImage = ImageType::New();
  typename ImageType::IndexType start;
  start.Fill(0);
  typename ImageType::SizeType size;
183
  size.Fill(n);
184
185
186
187
  typename ImageType::RegionType region(start, size);
  onePixelImage->SetRegions(region);
  onePixelImage->SetNumberOfComponentsPerPixel(m_InputImage->GetNumberOfComponentsPerPixel());
  onePixelImage->Allocate();
188
189

  // Instanciate and initialize a segmenter
190
191
  TSegmenter segmenter;
  segmenter.SetInput(onePixelImage);
192
  grm::GraphOperations<TSegmenter>::InitNodes(onePixelImage,segmenter,FOUR);
193

194
  // Get the memory occupied by the graph, normalize it by n*n
195
  unsigned int memory = segmenter.GetGraphMemory() / (n*n);
196

197
  itkDebugMacro(<<"Size of a node is " << memory);
198

199
  return memory;
200
}
remicres's avatar
remicres committed
201

202
203
204
/*
 * Compute the maximum number of nodes which can fit in the system memory
 */
remicres's avatar
remicres committed
205
template<class TSegmenter>
remicres's avatar
remicres committed
206
long unsigned int Controller<TSegmenter>::GetMaximumNumberOfNodesInMemory()
remicres's avatar
remicres committed
207
{
remicres's avatar
remicres committed
208
209
  itkDebugMacro(<< "Computing maximum number of nodes in memory");

remicres's avatar
remicres committed
210
211
212
213
214
  m_Memory = getMemorySize();
  assert(m_Memory > 0);

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

215
  return std::ceil(((float) m_Memory) / ((float) GetNodeMemory()));
remicres's avatar
remicres committed
216
217
218
219

}

template<class TSegmenter>
220
221
void
Controller<TSegmenter>::ComputeMaximumStabilityMargin(unsigned int width,
remicres's avatar
remicres committed
222
    unsigned int height, unsigned int &niter, unsigned int &margin)
223
    {
remicres's avatar
remicres committed
224
225
226
227
  itkDebugMacro(<< "Computing maximum stability margin");

  // Compute the stability margin. The naive strategy consider a margin value and a stable size equal.
  niter = 1;
228
229
230
  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
231

232
  while(currMargin < maxMargin)
remicres's avatar
remicres committed
233
    {
234
    margin = currMargin;
remicres's avatar
remicres committed
235
    niter++;
236
    currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
remicres's avatar
remicres committed
237
238
239
240
    }
  niter--;

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

242
    }
remicres's avatar
remicres committed
243

244
/*
245
 * Compute a tiling layout which minimizes a criterion based on tile compactness
remicres's avatar
remicres committed
246
 * and memory usage
247
248
 *
 * TODO: use the lsgrmSplitter to truly compute the largest tile of a given layout
249
 */
remicres's avatar
remicres committed
250
251
252
253
template<class TSegmenter>
void Controller<TSegmenter>::GetAutomaticConfiguration()
{

remicres's avatar
remicres committed
254
255
  itkDebugMacro(<<"Get automatic configuration");

remicres's avatar
remicres committed
256
  // Compute the maximum number of nodes that can fit the memory
remicres's avatar
remicres committed
257
  // TODO: Use the smallest number amongst MPI processes
remicres's avatar
remicres committed
258
  unsigned long int maximumNumberOfNodesInMemory = GetMaximumNumberOfNodesInMemory();
remicres's avatar
remicres committed
259
  itkDebugMacro(<<"Maximum number of nodes in memory is " << maximumNumberOfNodesInMemory);
remicres's avatar
remicres committed
260

remicres's avatar
remicres committed
261
  // Number of nodes in the entire image
remicres's avatar
remicres committed
262
263
  const unsigned int imageWidth = m_InputImage->GetLargestPossibleRegion().GetSize()[0];
  const unsigned int imageHeight = m_InputImage->GetLargestPossibleRegion().GetSize()[1];
remicres's avatar
remicres committed
264
265
  const unsigned long int nbOfNodesInImage = imageWidth*imageHeight;

266
  // Default layout: 1x1
remicres's avatar
remicres committed
267
268
  m_NbTilesX = 1;
  m_NbTilesY = 1;
269

remicres's avatar
remicres committed
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
  // 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++)
289
    {
remicres's avatar
remicres committed
290
291
    // 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
292
      {
remicres's avatar
remicres committed
293
294
#ifdef OTB_USE_MPI
      // We want number of tiles which is a multiple of the number of MPI processes
295
296
      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
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
#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;
          float compactness = perimeter / surface * (float) vcl_max(tileWidth,tileHeight); // [1,+inf]

          // 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
342
343
344
345

  // 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
346
347
348
349
350
351
  itkDebugMacro(<<"Selected layout: " << m_NbTilesX << "x" << m_NbTilesY
      << " (criterion=" << lowestCriterionValue << ")");

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

352
353
354
  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);
355
  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
356
      << (m_TileWidth + 2*m_Margin) << "x" << (m_TileHeight + 2*m_Margin) );
remicres's avatar
remicres committed
357

remicres's avatar
remicres committed
358
359
}

remicres's avatar
remicres committed
360
361
362
363
364
365
366
367
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
368
void Controller<TSegmenter>::SetInputImage(ImageType * inputImage)
remicres's avatar
remicres committed
369
{
remicres's avatar
remicres committed
370
  m_InputImage = inputImage;
remicres's avatar
remicres committed
371
372
373
374
375
376
377
378
}

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

remicres's avatar
remicres committed
379
380
381
382
383
384
template<class TSegmenter>
typename Controller<TSegmenter>::LabelImageType::Pointer
Controller<TSegmenter>::GetLabeledClusteredOutput()
{
  return m_LabelImage;
}
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399

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