lsgrmController.txx 13.2 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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
      // Time monitoring
      ShowTime(t);
      }


#ifdef OTB_USE_MPI
    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],
137
          m_InputImage->GetNumberOfComponentsPerPixel(),
remicres's avatar
remicres committed
138
          numberOfIterationsRemaining);
remicres's avatar
remicres committed
139

140
      //      ShowTime(t);
remicres's avatar
remicres committed
141
142

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

}

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

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

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

195
  itkDebugMacro(<<"Size of a node is " << memory);
196

197
  return memory;
198
}
remicres's avatar
remicres committed
199

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

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

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

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

}

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

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

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

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

240
    }
remicres's avatar
remicres committed
241

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

remicres's avatar
remicres committed
252
253
  itkDebugMacro(<<"Get automatic configuration");

remicres's avatar
remicres committed
254
255
  // Compute the maximum number of nodes that can fit the memory
  unsigned long int maximumNumberOfNodesInMemory = GetMaximumNumberOfNodesInMemory();
remicres's avatar
remicres committed
256
  itkDebugMacro(<<"Maximum number of nodes in memory is " << maximumNumberOfNodesInMemory);
remicres's avatar
remicres committed
257

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

263
  // Default layout: 1x1
remicres's avatar
remicres committed
264
265
  m_NbTilesX = 1;
  m_NbTilesY = 1;
266

remicres's avatar
remicres committed
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
  // 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++)
286
    {
remicres's avatar
remicres committed
287
288
    // 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
289
      {
remicres's avatar
remicres committed
290
291
#ifdef OTB_USE_MPI
      // We want number of tiles which is a multiple of the number of MPI processes
292
293
      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
294
295
296
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
#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
339
340
341
342

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

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

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

remicres's avatar
remicres committed
355
356
}

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

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

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

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