#ifndef _OTB_TRISKELE_SD_ATTRIBUTES_TPP
#define _OTB_TRISKELE_SD_ATTRIBUTES_TPP

using namespace boost::chrono;

inline
SDAttributes::SDAttributes (const Tree &tree, const AreaAttributes &areaAttributes)
  : CompAttribute<double> (tree) {
  compute (areaAttributes);
}

inline
SDAttributes::~SDAttributes () {
}

template<typename PixelT>
inline void
SDAttributes::cut (vector<vector<PixelT> > &allBands, const AttributeProfiles<PixelT> &attributeProfiles,
		   const vector<double> &thresholds) const {
  DEF_LOG ("SDAttributes::cut", "thresholds:" << thresholds.size ());
  if (!thresholds.size ())
    return;
  auto start = high_resolution_clock::now ();
  double maxValue = CompAttribute<double>::getMaxValue ();
  cerr << "sd max value:" << maxValue << endl;
  CompAttribute<double>::cut (allBands, attributeProfiles, 0,
			       CompAttribute<double>::getScaledThresholds (thresholds, maxValue));
  globalTreeStats.addTime (filteringStats, duration_cast<duration<double> > (high_resolution_clock::now ()-start).count ());
}

inline void
SDAttributes::compute (const AreaAttributes &areaAttributes) {
  auto start = high_resolution_clock::now ();
  const DimImg *areas = areaAttributes.getValues ();
  computeSameCompLevel ([this, &areas] (const DimImg &parentId) {
      double tmpSD = 0.;
      tree.forEachChildTI (parentId, [this, &tmpSD, &areas, &parentId] (const bool &isLeaf, const DimImg &childId) {
	  if (isLeaf)
	    return;
	  double diff = areas [childId]-areas [parentId];
	  tmpSD += (diff * diff)*areas[childId];
	});
      CompAttribute<double>::values[parentId] = tmpSD / areas[parentId];

      // XXX on peut virer pour gagner du temps et calculer SD*SD
      CompAttribute<double>::values[parentId] = sqrt (CompAttribute<double>::values[parentId]);
    });
  globalTreeStats.addTime (sdStats, duration_cast<duration<double> > (high_resolution_clock::now ()-start).count ());
}

#endif // _OTB_TRISKELE_SD_ATTRIBUTES_TPP