otbMultiImageSamplingRate.cxx 16.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbSamplingRateCalculatorList.h"
#include "otbStatisticsXMLFileReader.h"

namespace otb
{
namespace Wrapper
{

class MultiImageSamplingRate : public Application
{
public:
  /** Standard class typedefs. */
  typedef MultiImageSamplingRate        Self;
  typedef Application                   Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

  /** Standard macro */
  itkNewMacro(Self);

  itkTypeMacro(MultiImageSamplingRate, otb::Application);

  /** typedef */
  typedef otb::SamplingRateCalculatorList           RateCalculatorListType;
  
  typedef RateCalculatorListType::ClassCountMapType      ClassCountMapType;
49
  typedef RateCalculatorListType::MapRateType           MapRateType;
50 51 52 53 54 55 56 57 58 59

  typedef itk::VariableLengthVector<float> MeasurementType;
  typedef otb::StatisticsXMLFileReader<MeasurementType> XMLReaderType;

private:
  MultiImageSamplingRate()
    {
    m_CalculatorList = RateCalculatorListType::New();
    }

60
  void DoInit() override
61 62 63 64 65 66 67 68 69 70 71 72
  {
    SetName("MultiImageSamplingRate");
    SetDescription("Compute sampling rate for an input set of images.");

    // Documentation
    SetDocName("Multi-image sampling rate estimation");
    SetDocLongDescription("The application computes sampling rates for a set of"
      " input images. Before calling this application, each pair of image and "
      "training vectors has to be analysed with the application "
      "PolygonClassStatistics. The statistics file is then used to compute the "
      "sampling rates for each class in each image. Several types of sampling "
      " are implemented. Each one is a combination of a mono-image strategy "
73 74
      "and a multi-image mode. The mono-image strategies are :\n\n"
      "  * smallest (default) : select the same number of sample in each "
75
      "class so that the smallest one is fully sampled.\n"
76
      "  * constant : select the same number of samples N in each class "
77
      "(with N below or equal to the size of the smallest class).\n"
78
      "  * byclass : set the required number for each class manually, with an "
79
      "input CSV file (first column is class name, second one is the required "
80 81
      "samples number).\n\n"
      "The multi-image modes (mim) are proportional, equal and custom. The custom "
82 83 84 85
      "mode lets the users choose the distribution of samples among the "
      "images. The different behaviours are described below. Ti(c) and Ni(c) "
      " refers resp. to the total number and needed number of samples in "
      "image i for class c. Let's call L the total number of images.\n"
86 87 88 89 90 91 92 93 94 95 96 97 98 99
      "  * strategy = all\n\n"
      "    - Same behaviour for all modes : take all samples\n\n"
      "  * strategy = constant :"
      " let's call M the global number of samples required per class."
      " For each image i and each class c:\n\n"
      "    - if mim = proportional, then Ni( c ) = M * Ti( c ) / sum_k( Tk(c) )\n\n"
      "    - if mim = equal       , then Ni( c ) = M / L\n\n"
      "    - if mim = custom      , then Ni( c ) = Mi where Mi is the custom requested number of samples for image i\n\n"
      "  * strategy = byClass :"
      " let's call M(c) the global number of samples for class c)."
      " For each image i and each class c:\n\n"
      "    - if mim = proportional, then Ni( c ) = M(c) * Ti( c ) / sum_k( Tk(c) )\n\n"
      "    - if mim = equal       , then Ni( c ) = M(c) / L\n\n"
      "    - if mim = custom      , then Ni( c ) = Mi(c) where Mi(c) is the custom requested number of samples for image i and class c\n\n"
100 101 102 103 104 105 106 107 108 109
      "  * strategy = percent :"
      " For each image i and each class c:\n\n"
      "    - if mim = proportional, then Ni( c ) = p * Ti( c ) where p is the global percentage of samples\n\n"
      "    - if mim = equal       , then Ni( c ) = p * sum_k(Tk(c)]/L where p is the global percentage of samples\n\n"
      "    - if mim = custom      , then Ni( c ) = p(i) * Ti(c) where p(i) is the percentage of samples for image i. c\n\n"
      "  * strategy = total :"
      " For each image i and each class c:\n\n"
      "    - if mim = proportional, then Ni( c ) = total * (sum_k(Ti(k))/sum_kl(Tl(k))) * (Ti(c)/sum_k(Ti(k))) where total is the total number of samples specified.\n\n"
      "    - if mim = equal       , then Ni( c ) = (total / L) * (Ti(c)/sum_k(Ti(k))) where total is the total number of samples specified.\n\n"
      "    - if mim = custom      , then Ni( c ) = total(i) * (Ti(c)/sum_k(Ti(k))) where total(i) is the total number of samples specified for image i. \n\n"                          
110 111 112
      "  * strategy = smallest class\n\n"
      "    - if mim = proportional, then the smallest class size (computed globally) is used for the strategy constant+proportional.\n\n"
      "    - if mim = equal       , then the smallest class size (computed globally) is used for the strategy constant+equal.\n\n"
113
      "    - if mim = custom      , then the smallest class is computed and used for each image separately."
114 115 116 117 118 119 120
      );
    SetDocLimitations("None");
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso(" ");

    AddDocTag(Tags::Learning);

121
    AddParameter(ParameterType_InputFilenameList,  "il",   "Input statistics");
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
    SetParameterDescription("il", "List of statistics files for each input image.");

    AddParameter(ParameterType_OutputFilename, "out", "Output sampling rates");
    SetParameterDescription("out","Output filename storing sampling rates (CSV "
      "format with class name, required samples, total samples, and rate). "
      "The given filename will be used with a suffix to indicate the "
      "corresponding input index (for instance: rates.csv will give rates_1.csv"
      ", rates_2.csv, ...).");

    AddParameter(ParameterType_Choice, "strategy", "Sampling strategy");

    AddChoice("strategy.byclass","Set samples count for each class");
    SetParameterDescription("strategy.byclass","Set samples count for each class");

    AddParameter(ParameterType_InputFilenameList, "strategy.byclass.in", "Number of samples by class");
    SetParameterDescription("strategy.byclass.in", "Number of samples by class "
      "(CSV format with class name in 1st column and required samples in the 2nd)."
      "In the case of the custom multi-image mode, several inputs may be given for each image.");

    AddChoice("strategy.constant","Set the same samples counts for all classes");
    SetParameterDescription("strategy.constant","Set the same samples counts for all classes");

    AddParameter(ParameterType_String, "strategy.constant.nb", "Number of samples for all classes");
    SetParameterDescription("strategy.constant.nb", "Number of samples for all classes."
      "In the case of the custom multi-image mode, several values can be given for each image.");

    AddChoice("strategy.smallest","Set same number of samples for all classes, with the smallest class fully sampled");
    SetParameterDescription("strategy.smallest","Set same number of samples for all classes, with the smallest class fully sampled");

151 152 153 154 155 156 157 158 159 160 161 162
    AddChoice("strategy.percent","Use a percentage of the samples available for each class");
    SetParameterDescription("strategy.percent","Use a percentage of the samples available for each class");

    AddParameter(ParameterType_String,"strategy.percent.p","The percentage(s) to use");
    SetParameterDescription("strategy.percent.p","The percentage(s) to use " "In the case of the custom multi-image mode, several values can be given for each image.");

    AddChoice("strategy.total","Set the total number of samples to generate, and use class proportions.");
    SetParameterDescription("strategy.total","Set the total number of samples to generate, and use class proportions.");

    AddParameter(ParameterType_String,"strategy.total.v","The number of samples to generate");
    SetParameterDescription("strategy.total.v","The number of samples to generate" "In the case of the custom multi-image mode, several values can be given for each image.");
 
163 164 165 166
    AddChoice("strategy.all","Take all samples");
    SetParameterDescription("strategy.all","Take all samples");

    // Default strategy : smallest
167
    SetParameterString("strategy","smallest");
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183

    AddParameter(ParameterType_Choice, "mim", "Multi-Image Mode");

    AddChoice("mim.proportional", "Proportional");
    SetParameterDescription("mim.proportional","Split proportionally the required number of samples");

    AddChoice("mim.equal", "equal");
    SetParameterDescription("mim.equal","Split equally the required number of samples");

    AddChoice("mim.custom", "Custom");
    SetParameterDescription("mim.custom","Split the required number of samples following user choice.");

    // Doc example parameter settings
    SetDocExampleParameterValue("il", "stats_1.xml stats_2.xml");
    SetDocExampleParameterValue("out", "rates.csv");
    SetDocExampleParameterValue("strategy", "smallest");
184
    SetDocExampleParameterValue("mim","proportional");
185

186
    SetOfficialDocLink();
187 188
  }

189
  void DoUpdateParameters() override
190 191 192
  {
  }

193
  void DoExecute() override
194 195
    {
    // Clear state
196
    m_CalculatorList->Clear();
197 198 199 200 201
    std::vector<std::string> inputs = this->GetParameterStringList("il");
    unsigned int nbInputs = inputs.size();
    XMLReaderType::Pointer statReader = XMLReaderType::New();
    for (unsigned int i=0 ; i<nbInputs ; i++ )
      {
202
      m_CalculatorList->PushBack(otb::SamplingRateCalculator::New());
203 204 205 206 207 208
      statReader->SetFileName(inputs[i]);
      ClassCountMapType classCount = statReader->GetStatisticMapByName<ClassCountMapType>("samplesPerClass");
      m_CalculatorList->SetNthClassCount(i,classCount);
      }

    // Cautions : direct mapping between the enum PartitionType and the choice order
209 210
    RateCalculatorListType::PartitionType partitionMode =
      static_cast<RateCalculatorListType::PartitionType>(this->GetParameterInt("mim"));
211 212 213 214 215

    unsigned int minParamSize = 1;
    if (partitionMode == RateCalculatorListType::CUSTOM)
      {
      // Check we have enough inputs for the custom mode
216
      minParamSize = nbInputs;
217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
      }
    
    switch (this->GetParameterInt("strategy"))
      {
      // byclass
      case 0:
        {
        std::vector<std::string> requiredFiles = this->GetParameterStringList("strategy.byclass.in");
        std::vector<ClassCountMapType> requiredCounts;
        if (requiredFiles.size() < minParamSize)
          {
          otbAppLogFATAL("Missing arguments in strategy.byclass.in to process sampling rates");
          }
        otbAppLogINFO("Sampling strategy : set number of samples for each class");
        for (unsigned int i=0 ; i<minParamSize ; i++)
          {
          requiredCounts.push_back(otb::SamplingRateCalculator::ReadRequiredSamples(requiredFiles[i]));
          }
        m_CalculatorList->SetNbOfSamplesByClass(requiredCounts, partitionMode);
        }
      break;
      // constant
      case 1:
        {
        std::vector<itksys::String> parts = itksys::SystemTools::SplitString(this->GetParameterString("strategy.constant.nb"),' ');
        std::vector<unsigned long> countList;
        for (unsigned int i=0 ; i<parts.size() ; i++)
          {
          if (!parts[i].empty())
            {
            std::string::size_type pos1 = parts[i].find_first_not_of(" \t");
            std::string::size_type pos2 = parts[i].find_last_not_of(" \t");
            std::string value(parts[i].substr(pos1, pos2 - pos1 + 1));
            countList.push_back(boost::lexical_cast<unsigned long>(parts[i]));
            }
          }
        if (countList.size() < minParamSize)
          {
          otbAppLogFATAL("Missing arguments in strategy.constant.nb to process sampling rates");
          }
        otbAppLogINFO("Sampling strategy : set a constant number of samples for all classes");
        m_CalculatorList->SetNbOfSamplesAllClasses(countList, partitionMode);
        }
      break;
261
      // smallest class
262
      case 2:
263 264 265 266 267 268 269
        {
        otbAppLogINFO("Sampling strategy : fit the number of samples based on the smallest class");
        m_CalculatorList->SetMinimumNbOfSamplesByClass(partitionMode);
        }
      break;
      // percent
      case 3:
270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297
        {
        std::vector<itksys::String> parts = itksys::SystemTools::SplitString(this->GetParameterString("strategy.percent.p"),' ');
        std::vector<double> percentList;
        for (unsigned int i=0 ; i<parts.size() ; i++)
          {
          if (!parts[i].empty())
            {
            std::string::size_type pos1 = parts[i].find_first_not_of(" \t");
            std::string::size_type pos2 = parts[i].find_last_not_of(" \t");
            std::string value(parts[i].substr(pos1, pos2 - pos1 + 1));
            percentList.push_back(boost::lexical_cast<double>(parts[i]));

            if(percentList.back()<0 || percentList.back()>1)
              {
              otbAppLogFATAL("Percent parameter should be in range [0,1]");
              } 
            }
          }
        if (percentList.size() < minParamSize)
          {
          otbAppLogFATAL("Missing arguments in strategy.percent.p to process sampling rates");
          }
        otbAppLogINFO("Sampling strategy : set a percentage of samples to be used.");
        m_CalculatorList->SetPercentageOfSamples(percentList, partitionMode);        
        }
      break;

      // total
298
      case 4:
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
        {
        std::vector<itksys::String> parts = itksys::SystemTools::SplitString(this->GetParameterString("strategy.total.v"),' ');
        std::vector<unsigned long> totalList;
        for (unsigned int i=0 ; i<parts.size() ; i++)
          {
          if (!parts[i].empty())
            {
            std::string::size_type pos1 = parts[i].find_first_not_of(" \t");
            std::string::size_type pos2 = parts[i].find_last_not_of(" \t");
            std::string value(parts[i].substr(pos1, pos2 - pos1 + 1));
            totalList.push_back(boost::lexical_cast<unsigned long>(parts[i]));
            }
          }
        if (totalList.size() < minParamSize)
          {
          otbAppLogFATAL("Missing arguments in strategy.total.v to process sampling rates");
          }
        otbAppLogINFO("Sampling strategy : set a constant number of samples for all classes");
        m_CalculatorList->SetTotalNumberOfSamples(totalList, partitionMode);
        }
      break;
320
      // all samples
321
      case 5:
322 323 324 325 326 327 328 329 330 331
        {
        otbAppLogINFO("Sampling strategy : take all samples");
        m_CalculatorList->SetAllSamples(partitionMode);
        }
      break;
      default:
        otbAppLogFATAL("Strategy mode unknown :"<<this->GetParameterString("strategy"));
      break;
      }

332
    std::ostringstream oss;    
333
    std::string outputPath(this->GetParameterString("out"));
334 335
    std::string outputBase = outputPath.substr(0, outputPath.find_last_of('.'));
    std::string outputExt = outputPath.substr(outputPath.find_last_of('.'), std::string::npos);
336
    unsigned int overflowCount = 0;
337
    bool noSamples=true;
338 339
    for (unsigned int i=0 ; i<nbInputs ; i++ )
      {
340
      // Print results
341
      oss.str(std::string(""));
342
      oss << " className  requiredSamples  totalSamples  rate\n";
343
      MapRateType rates = m_CalculatorList->GetRatesByClass(i);
344
      if(!rates.empty())
345 346 347
        {
        noSamples = false;
        }
348 349 350 351 352 353 354
      MapRateType::const_iterator itRates = rates.begin();
      for(; itRates != rates.end(); ++itRates)
        {
        otb::SamplingRateCalculator::TripletType tpt = itRates->second;
        oss << itRates->first << "\t" << tpt.Required << "\t" << tpt.Tot << "\t" << tpt.Rate;
        if (tpt.Required > tpt.Tot)
          {
355
          overflowCount++;
356 357 358 359 360 361 362 363 364
          oss << "\t[OVERFLOW]";
          }
        oss << std::endl;
        }
      otbAppLogINFO("Sampling rates for image "<< i+1 <<" : " << oss.str());
      // Output results to disk
      oss.str(std::string(""));
      oss << outputBase << "_" << i+1 << outputExt;
      m_CalculatorList->GetNthElement(i)->Write(oss.str());
365
      }
366 367 368 369
    if (noSamples)
      {
      otbAppLogFATAL("No samples found in the inputs!");
      }
370 371 372 373 374
    if (overflowCount)
      {
      std::string plural(overflowCount>1?"s":"");
      otbAppLogWARNING(<< overflowCount << " case"<<plural<<" of overflow detected! (requested number of samples higher than total available samples)");
      }
375 376 377 378 379 380 381 382 383
  }

  RateCalculatorListType::Pointer m_CalculatorList;
};

} // end namespace Wrapper
} // end namespace otb

OTB_APPLICATION_EXPORT(otb::Wrapper::MultiImageSamplingRate)