hitfinder.cpp

Go to the documentation of this file.

// Copyright (C) 2013,2018 Lutz Foucar

/** @file hitfinder.cpp contains processors that will extract pixels of
 *                      interrest from 2d histograms.
 * @author Lutz Foucar
 */

#include <QtCore/QString>
#include <iterator>
#include <algorithm>

#include "cass.h"
#include "hitfinder.h"
#include "convenience_functions.h"
#include "cass_settings.h"
#include "log.h"
#include "geom_parser.h"

using namespace std;
using namespace cass;
using tr1::bind;
using tr1::placeholders::_1;
using tr1::placeholders::_2;


// ********** processor 203: subtract local background ************

pp203::pp203(const name_t &name)
  : Processor(name)
{
  loadSettings(0);
}

void pp203::loadSettings(size_t)
{
  CASSSettings s;

  s.beginGroup("Processor");
  s.beginGroup(QString::fromStdString(name()));

  // size of box for median
  _boxSize = make_pair(s.value("BoxSizeX", 10).toUInt(),
                       s.value("BoxSizeY",10).toUInt());
  _sectionSize = make_pair(s.value("SectionSizeX", 1024).toUInt(),
                           s.value("SectionSizeY",512).toUInt());

  setupGeneral();

  // Get the input
  _hist = setupDependency("ImageName");
  bool ret (setupCondition());
  if (!(_hist && ret)) return;

  // Create the output
  createHistList(_hist->result().clone());

  Log::add(Log::INFO,"Processor '" + name() +
           "' removes median of a box with size x'" + toString(_boxSize.first) +
           "' and y '" + toString(_boxSize.second) + "' from individual pixels" +
           " of hist in pp '"+ _hist->name() + "'. Condition is '" +
           _condition->name() + "'");
}

void pp203::process(const CASSEvent& evt, result_t &result)
{
  // Get the input
  const result_t &image(_hist->result(evt.id()));
  QReadLocker lock(&image.lock);

  const size_t xsize(image.axis(result_t::xAxis).nBins);
  const size_t ysize(image.axis(result_t::yAxis).nBins);
  const size_t xsectionsize(_sectionSize.first);
  const size_t ysectionsize(_sectionSize.second);
  const size_t xnbrsections(xsize/xsectionsize);
  const size_t ynbrsections(ysize/ysectionsize);
  const size_t xboxsize(_boxSize.first);
  const size_t yboxsize(_boxSize.second);

  vector<float> box;
  for (size_t sy=0; sy<ynbrsections; ++sy)
  {
    for (size_t sx=0; sx<xnbrsections; ++sx)
    {
      const size_t xsectionbegin(sx*xsectionsize);
      const size_t xsectionend(sx*xsectionsize + xsectionsize);
      const size_t ysectionbegin(sy*ysectionsize);
      const size_t ysectionend(sy*ysectionsize + ysectionsize);

//      cout << xsectionbegin << " "<<xsectionend << " "<<ysectionbegin << " "<<ysectionend<<endl;
      for (size_t y=ysectionbegin; y<ysectionend; ++y)
      {
        for (size_t x=xsectionbegin; x<xsectionend; ++x)
        {
          const size_t pixAddr(y*xsize+x);
          const float pixel(image[pixAddr]);

          if (fuzzycompare(pixel,0.f) )
          {
            result[pixAddr] = pixel;
            continue;
          }

          const size_t xboxbegin(max(static_cast<int>(xsectionbegin),static_cast<int>(x)-static_cast<int>(xboxsize)));
          const size_t xboxend(min(xsectionend,x+xboxsize));
          const size_t yboxbegin(max(static_cast<int>(ysectionbegin),static_cast<int>(y)-static_cast<int>(yboxsize)));
          const size_t yboxend(min(ysectionend,y+yboxsize));

//          cout <<x << " " << xboxbegin << " "<<xboxend <<" : " <<y<< " "<<yboxbegin << " "<<yboxend<<endl;
          box.clear();
          for (size_t yb=yboxbegin; yb<yboxend;++yb)
          {
            for (size_t xb=xboxbegin; xb<xboxend;++xb)
            {
              const size_t pixAddrBox(yb*xsize+xb);
              const float pixel_box(image[pixAddrBox]);
              if ( ! fuzzyIsNull(pixel_box) )
                box.push_back(pixel_box);
            }
          }

          if (box.empty())
          {
            result[pixAddr] = pixel;
          }
          else
          {
            const size_t mid(0.5*box.size());
            nth_element(box.begin(), box.begin() + mid, box.end());
            const float median = box[mid];
            result[pixAddr] = median;
//            const float backgroundsubstractedPixel(pixel-median);
//            image_out[pixAddr] = backgroundsubstractedPixel;
          }
        }
      }
    }
  }
}



// ********** processor 204: find bragg peaks ************

pp204::pp204(const name_t &name)
  : Processor(name)
{
  loadSettings(0);
}

void pp204::loadSettings(size_t)
{
  CASSSettings s;

  s.beginGroup("Processor");
  s.beginGroup(QString::fromStdString(name()));

  // size of box for median
  const int peakRadius(s.value("BraggPeakRadius",2).toInt());
  _peakRadiusSq = peakRadius*peakRadius;
  const int goodBoxSize(sqrt(3.1415) * peakRadius);
  _box = make_pair(s.value("BoxSizeX", goodBoxSize).toUInt(),
                   s.value("BoxSizeY",goodBoxSize).toUInt());
  _section = make_pair(s.value("SectionSizeX", 1024).toUInt(),
                       s.value("SectionSizeY",512).toUInt());
  _threshold = s.value("Threshold",300).toFloat();
  _minSnr = s.value("MinSignalToNoiseRatio",20).toFloat();
  _minNeighbourSNR = s.value("MinNeighbourSNR",3).toFloat();
  _minBckgndPixels = s.value("MinNbrBackgrndPixels",10).toInt();

  setupGeneral();

  _hist = setupDependency("ImageName");

  bool ret (setupCondition());
  if (!(_hist && ret))
    return;

  /** Create the resulting table */
  createHistList(result_t::shared_pointer(new result_t(nbrOf,0)));

  /** log what the user was requesting */
  string output("Processor '" + name() + "' finds bragg peaks." +
                "'. Boxsize '" + toString(_box.first)+"x"+ toString(_box.second)+
                "'. SectionSize '" + toString(_section.first)+"x"+ toString(_section.second)+
                "'. Threshold '" + toString(_threshold) +
                "'. MinSignalToNoiseRatio '" + toString(_minSnr) +
                "'. MinNbrBackgrndPixels '" + toString(_minBckgndPixels) +
                "'. Square BraggPeakRadius '" + toString(_peakRadiusSq) +
                "'. Using input histogram :" + _hist->name() +
                "'. Condition is '" + _condition->name() + "'");
  Log::add(Log::INFO,output);
}


int pp204::getBoxStatistics(result_t::const_iterator pixel,
                            const imagepos_t ncols,
                            pixelval_t &mean, pixelval_t &stdv, int &count)
{
  enum{good,skip};
  CummulativeStatisticsCalculator<pixelval_t> stat;
  for (imagepos_t bRow=-_box.second; bRow <= _box.second; ++bRow)
  {
    for (imagepos_t bCol=-_box.first; bCol <= _box.first; ++bCol)
    {
      /** only consider pixels that are not in the origin of the box */
      if (bRow == 0 && bCol == 0)
        continue;

      /** check if the current box pixel value is bigger than the center pixel
       *  value. If so skip this pixel
       */
      const imagepos_t bLocIdx(bRow*ncols+bCol);
      const pixelval_t bPixel(pixel[bLocIdx]);
      if(*pixel < bPixel )
        return skip;

      /** if box pixel is outside the radius, add it to the statistics, if it
       *  is inside the radius bad then skip, as no pixel that could potentially
       *  be part of the peak should be a bad pixel.
       */
      const bool pixIsBad(qFuzzyIsNull(bPixel));
      const int radiussq(bRow*bRow + bCol*bCol);
      if (_peakRadiusSq < radiussq)
        stat.addDatum(bPixel);
      else if (pixIsBad)
        return skip;
    }
  }
  count = stat.count();
  mean = stat.mean();
  stdv = stat.stdv();
  return good;
}

void pp204::process(const CASSEvent& evt, result_t &result)
{
  const result_t &image(_hist->result(evt.id()));
  QReadLocker lock(&image.lock);

  /** reset the nbr of y entries */
  result.resetTable();

  const imagepos_t ncols(image.shape().first);
  const imagepos_t nrows(image.shape().second);

  imagepos_t neigbourOffsetArray[] =
  {
    ncols-1,
    ncols,
    ncols+1,
    -1,
    1,
    -ncols-1,
    -ncols,
    -ncols+1,
  };
  vector<imagepos_t> neighboursOffsets(neigbourOffsetArray,
                                neigbourOffsetArray + sizeof(neigbourOffsetArray)/sizeof(int));
  table_t peak(nbrOf,0);
  vector<bool> checkedPixels(image.size(),false);

  vector<bool>::iterator checkedPixel(checkedPixels.begin());
  result_t::const_iterator pixel(image.begin());
  imagepos_t idx(0);
  for (;pixel != image.end(); ++pixel,++idx, ++checkedPixel)
  {
    /** check if pixel has been touched before */
    if (*checkedPixel)
      continue;

    /** check above a set threshold */
    if (*pixel < _threshold)
      continue;

    /** get coordinates of pixel and tell that it is touched */
    const uint16_t col(idx % ncols);
    const uint16_t row(idx / ncols);

    /** make sure that pixel is located such that the box will not conflict with
     *  the image and section boundaries
     */
    if (col < _box.first  || ncols - _box.first  < col || //within box in x
        row < _box.second || nrows - _box.second < row || //within box in y
        (col - _box.first)  / _section.first  != (col + _box.first)  / _section.first || //x within same section
        (row - _box.second) / _section.second != (row + _box.second) / _section.second)  //y within same section
      continue;

    /** check wether current pixel value is highest and generate background values */
    pixelval_t mean,stdv;
    int count;
    if (getBoxStatistics(pixel,ncols,mean,stdv,count))
      continue;

    /** check that we took more than the minimum pixels for the background */
    if (count < _minBckgndPixels)
      continue;

    /** check if signal to noise ration is good for the current pixel */
    const pixelval_t snr((*pixel - mean) / stdv);
    if (snr < _minSnr)
      continue;

    /** from the central pixel look around and see which pixels are direct
     *  part of the peak. If a neighbour is found, mask it such that one does
     *  not use it twice.
     *
     *  Create a list that should contain the indizes of the pixels that are
     *  part of the peak. Go through that list and check for each found
     *  neighbour whether it also has a neighbour. If so add it to the list, but
     *  only if it is part of the box.
     */
    vector<imagepos_t> peakIdxs;
    peakIdxs.push_back(idx);
    *checkedPixel = true;
    for (size_t pix=0; pix < peakIdxs.size(); ++pix)
    {
      const imagepos_t pixpos(peakIdxs[pix]);
      vector<imagepos_t>::const_iterator nOffset(neighboursOffsets.begin());
      vector<imagepos_t>::const_iterator neighboursEnd(neighboursOffsets.end());
      while(nOffset != neighboursEnd)
      {
        const size_t nIdx(pixpos + *nOffset++);
        const imagepos_t nCol(nIdx % ncols);
        const imagepos_t nRow(nIdx / ncols);
        if (checkedPixels[nIdx] ||
            _box.first < abs(col - nCol) ||
            _box.second < abs(row - nRow))
          continue;
        const pixelval_t nPixel(image[nIdx]);
        const pixelval_t nPixelWOBckgnd(nPixel - mean);
        const pixelval_t nSNR(nPixelWOBckgnd / stdv);
        if (_minNeighbourSNR < nSNR)
        {
          peakIdxs.push_back(nIdx);
          checkedPixels[nIdx] = true;
        }
      }
    }

    /** go through all pixels in the box, find out which pixels are part of the
     *  peak and centroid them. Mask all pixels in the box as checked so one
     *  does not check them again
     */
    pixelval_t integral = 0;
    pixelval_t weightCol = 0;
    pixelval_t weightRow = 0;
    imagepos_t nPix = 0;
    imagepos_t max_radiussq=0;
    imagepos_t min_radiussq=max(_box.first,_box.second)*max(_box.first,_box.second);
    for (int bRow=-_box.second; bRow <= _box.second; ++bRow)
    {
      for (int bCol=-_box.first; bCol <= _box.first; ++bCol)
      {
        const imagepos_t bLocIdx(bRow*ncols+bCol);
        const pixelval_t bPixel(pixel[bLocIdx]);
        const pixelval_t bPixelWOBckgnd(bPixel - mean);
        if (checkedPixel[bLocIdx])
        {
          const imagepos_t radiussq(bRow*bRow + bCol*bCol);
          if (radiussq > max_radiussq)
            max_radiussq = radiussq;
          if (radiussq < min_radiussq)
            min_radiussq = radiussq;
          integral += bPixelWOBckgnd;
          weightCol += (bPixelWOBckgnd * (bCol+col));
          weightRow += (bPixelWOBckgnd * (bRow+row));
          ++nPix;
        }
        checkedPixel[bLocIdx] = true;
      }
    }
    peak[centroidColumn] = weightCol / integral;
    peak[centroidRow] = weightRow / integral;
    peak[Intensity] = integral;
    peak[nbrOfPixels] = nPix;
    peak[SignalToNoise] = snr;
    peak[MaxRadius] = sqrt(max_radiussq);
    peak[MinRadius] = sqrt(min_radiussq);
    peak[Index] = idx;
    peak[Column] = col;
    peak[Row] = row;
    peak[MaxADU] = *pixel;
    peak[LocalBackground] = mean;
    peak[LocalBackgroundDeviation] = stdv;
    peak[nbrOfBackgroundPixels] = count;

    result.appendRows(peak);

  }

}







// ********** processor 205: display peaks ************

pp205::pp205(const name_t &name)
  : Processor(name)
{
  loadSettings(0);
}

void pp205::loadSettings(size_t)
{
  CASSSettings s;

  s.beginGroup("Processor");
  s.beginGroup(QString::fromStdString(name()));

  setupGeneral();

  // Get the input
  _hist = setupDependency("ImageName");
  _table = setupDependency("TableName");
  bool ret (setupCondition());
  if (!(_hist && ret && _table))
    return;

  _boxsize = make_pair(s.value("BoxSizeX", 10).toUInt(),
                       s.value("BoxSizeY",10).toUInt());
  _drawVal = s.value("DrawPixelValue",16000.f).toFloat();
  _drawInner = s.value("DrawInnerPixel",false).toBool();
  _drawInnerValue = s.value("InnerPixelValue",16000.f).toFloat();
  _radius = s.value("Radius",2.f).toFloat();
  _idxCol = s.value("IndexColumn").toUInt();
  _drawCircle = s.value("DrawCircle",true).toBool();
  _drawBox = s.value("DrawBox",true).toBool();

  size_t maxIdx(_table->result().axis(result_t::xAxis).nBins);
  if (_idxCol >= maxIdx)
    throw runtime_error("pp205::loadSettings(): '" + name() + "' The requested " +
                        "column '" + toString(_idxCol) + " 'exeeds the " +
                        "maximum possible value of '" + toString(maxIdx) + "'");

  // Create the output
  createHistList(_hist->result().clone());

  Log::add(Log::INFO,"Processor '" + name() +
           "' displays the peaks that are listed in '" + _table->name() +
           "' and that were found in '" + _hist->name() + "'. Condition is '" +
           _condition->name() + "'");
}

void pp205::process(const CASSEvent& evt, result_t &result)
{
  // Get the input
  const result_t &image(_hist->result(evt.id()));
  QReadLocker lock1(&image.lock);
  const result_t &table(_table->result(evt.id()));
  QReadLocker lock2(&table.lock);

  /** copy the input image to the resulting image */
  copy(image.begin(),image.end(),result.begin());

  /** extract the nbr of rows and columns of this table and the image */
  const size_t nTableCols(table.shape().first);
  const size_t nTableRows(table.shape().second);
  const size_t nImageCols(image.shape().first);

  /** go through all rows in table */
  for (size_t row=0; row < nTableRows; ++row)
  {
    /** extract the column with the global index of the peak center */
    const int idx(table[row*nTableCols + _idxCol]);
    result_t::iterator centerpixel(result.begin()+idx);

    /** draw user requested info (extract other stuff from table) */
    //box
    if (_drawBox)
    {
//      //lower row
//      for (int bCol=-_boxsize.first; bCol <= _boxsize.first; ++bCol)
//      {
//        const int bRow = -_boxsize.second;
//        const int bLocIdx(bRow*nImageCols+bCol);
//        centerpixel[bLocIdx] = _drawVal;
//      }
//      //upper row
//      for (int bCol=-_boxsize.first; bCol <= _boxsize.first; ++bCol)
//      {
//        const int bRow = _boxsize.second;
//        const int bLocIdx(bRow*nImageCols+bCol);
//        centerpixel[bLocIdx] = _drawVal;
//      }
//      //left col
//      for (int bRow=-_boxsize.second; bRow <= _boxsize.second; ++bRow)
//      {
//        const int bCol = -_boxsize.first;
//        const int bLocIdx(bRow*nImageCols+bCol);
//        centerpixel[bLocIdx] = _drawVal;
//      }
//      //right col
//      for (int bRow=-_boxsize.second; bRow <= _boxsize.second; ++bRow)
//      {
//        const int bCol = _boxsize.first;
//        const int bLocIdx(bRow*nImageCols+bCol);
//        centerpixel[bLocIdx] = _drawVal;
//      }

      for (int bRow = -_boxsize.second; bRow <= _boxsize.second; ++bRow)
      {
        for (int bCol = -_boxsize.first; bCol <= _boxsize.first; ++bCol)
        {
          const int bLocIdx(idx + bRow*nImageCols+bCol);
          if (bLocIdx >= static_cast<int>(result.size()) || bLocIdx < 0)
          {
            string tableRow;
            for (size_t i(0); i<nTableCols; ++i)
                tableRow += toString(i) + "'" +
                            toString(table[row*nTableCols + i])
                            +"', ";
            throw logic_error("pp205 '" + name() +
                              "': Calculated index out of bounds: " +
                              "nTableRows'" + toString(nTableRows) + "', " +
                              "tableRow'" + toString(row) + "', " +
                              "tableRowContent'" + tableRow + "', " +
                              "bRow'" + toString(bRow) + "', " +
                              "bCol'" + toString(bCol) + "', " +
                              "size'" + toString(result.size()) + "', " +
                              "idx'" + toString(idx) + "', " +
                              "nImageCols'" + toString(nImageCols) + "', " +
                              "localIdx'" + toString(bLocIdx) + "'");
          }
          const bool border = (bCol == _boxsize.first ||
                               bCol == -_boxsize.first ||
                               bRow == _boxsize.second ||
                               bRow == -_boxsize.second);
          if (border)
            result[bLocIdx] = _drawVal;
          else if (_drawInner)
            result[bLocIdx] = _drawInnerValue;
        }
      }
    }

    //circle
    if (_drawCircle)
    {
      for(size_t angle_deg = 0; angle_deg <360; ++angle_deg)
      {
        const float angle_rad(3.14 * static_cast<float>(angle_deg)/180.);
        const int cCol (static_cast<size_t>(round(_radius*sin(angle_rad))));
        const int cRow (static_cast<size_t>(round(_radius*cos(angle_rad))));
        const int cLocIdx(cRow*nImageCols+cCol);
        centerpixel[cLocIdx] = _drawVal;
      }
    }

  }
}






// ********** processor 206: find pixels of bragg peaks ************

pp206::pp206(const name_t &name)
  : Processor(name)
{
  loadSettings(0);
}

void pp206::loadSettings(size_t)
{
  CASSSettings s;

  s.beginGroup("Processor");
  s.beginGroup(QString::fromStdString(name()));

  // size of box for median
  _box = make_pair(s.value("BoxSizeX", 10).toUInt(),
                   s.value("BoxSizeY",10).toUInt());
  _section = make_pair(s.value("SectionSizeX", 1024).toUInt(),
                       s.value("SectionSizeY",512).toUInt());
  _threshold = s.value("Threshold",300).toFloat();
  _multiplier = s.value("Multiplier",4).toFloat();

  setupGeneral();

  _imagePP = setupDependency("ImageName");
  _noisePP = setupDependency("NoiseName");

  bool ret (setupCondition());
  if (!(_imagePP && ret && _noisePP))
    return;

  /** Create the result output */
  createHistList(result_t::shared_pointer(new result_t(pp204::nbrOf,0)));

  /** log what the user was requesting */
  string output("Processor '" + name() + "' finds bragg peaks." +
                "'. Boxsize '" + toString(_box.first)+"x"+ toString(_box.second)+
                "'. SectionSize '" + toString(_section.first)+"x"+ toString(_section.second)+
                "'. Threshold '" + toString(_threshold) +
                "'. Image Histogram :" + _imagePP->name() +
                "'. Noise Histogram :" + _noisePP->name() +
                "'. Condition is '" + _condition->name() + "'");
  Log::add(Log::INFO,output);
}

void pp206::process(const CASSEvent& evt, result_t &result)
{
  const result_t &image(_imagePP->result(evt.id()));
  QReadLocker lock1(&(image.lock));
  const result_t &noisemap(_noisePP->result(evt.id()));
  QReadLocker lock2(&(noisemap.lock));

  /** reset the table */
  result.resetTable();


  table_t peak(pp204::nbrOf,0);

  result_t::const_iterator pixel(image.begin());
  result_t::const_iterator imageEnd(image.end());
  result_t::const_iterator noise(noisemap.begin());
  size_t idx(0);
  vector<float> box;
  const uint16_t ncols(image.shape().first);
  const uint16_t nrows(image.shape().second);
  for (; pixel != imageEnd; ++pixel, ++noise, ++idx)
  {
//    if(*noise * _multiplier < *pixel)
    if (_threshold < *pixel)
    {
      const uint16_t x(idx % ncols);
      const uint16_t y(idx / ncols);

      const uint16_t xboxbegin(max(static_cast<int>(0),static_cast<int>(x)-static_cast<int>(_box.first)));
      const uint16_t xboxend(min(ncols,static_cast<uint16_t>(x+_box.first)));
      const uint16_t yboxbegin(max(static_cast<int>(0),static_cast<int>(y)-static_cast<int>(_box.second)));
      const uint16_t yboxend(min(nrows,static_cast<uint16_t>(y+_box.second)));

      box.clear();
      for (size_t yb=yboxbegin; yb<yboxend;++yb)
      {
        for (size_t xb=xboxbegin; xb<xboxend;++xb)
        {
          const size_t pixAddrBox(yb*ncols+xb);
          const float pixel_box(image[pixAddrBox]);
          /** check if current sourrounding pixel is a bad pixel (0.),
           *  if so we should disregard the pixel as a candiate and check the
           *  next pixel.
           */
          if (qFuzzyCompare(pixel_box,0.f) )
            goto NEXTPIXEL;
          else
            box.push_back(pixel_box);
        }
      }

      if (box.size() > 1)
      {
        const size_t mid(0.5*box.size());
        nth_element(box.begin(), box.begin() + mid, box.end());
        const float bckgnd = box[mid];
        const float clrdpixel(*pixel - bckgnd);
//        if(*noise * _multiplier < clrdpixel)
        if (_threshold < clrdpixel)
        {
          peak[pp204::Column] = x;
          peak[pp204::Row] = y;
          peak[pp204::MaxADU] = *pixel;
          peak[pp204::Intensity] = clrdpixel;

          result.appendRows(peak);
        }
      }
NEXTPIXEL:;
    }
  }
}









// ********** processor 208: find bragg peaks ************

pp208::pp208(const name_t &name)
  : Processor(name)
{
  loadSettings(0);

}

void pp208::loadSettings(size_t)
{
  CASSSettings s;

  s.beginGroup("Processor");
  s.beginGroup(QString::fromStdString(name()));

  _section = make_pair(s.value("SectionSizeX", 1024).toUInt(),
                       s.value("SectionSizeY",512).toUInt());
  _minSnr = s.value("MinSignalToNoiseRatio",4).toFloat();
  _minRatio = s.value("MinRatio",3).toFloat();

  const int peakDiameter(s.value("BraggPeakDiameter",2).toInt());
  /** area of box should at least be ratio times area of peak.
   * \f{eqnarray*}{
   * A_{box}&=& ratio \times A_{circle} \\
   * size^2 &=& ratio \times \pi r^2 \\
   * size   &=& \sqrt{\pi r^2 \times ratio} \\
   *        &=& r \sqrt{\pi \times ratio} \\
   *        &=& 0.5 d \sqrt{\pi \times ratio}
   * \f}
   *  If \f$ size \f$ goes from  \f$-s ... s\f$ then \f$ size = 2s+1 \f$
   * \f{eqnarray*}{
   * 2s + 1 &=& 0.5 d \sqrt{\pi \times ratio} \\
   *      s &=& 0.5 \times 0.5 d \sqrt{\pi \times ratio} - 0.5 \\
   *        &=& 0.25 d \sqrt{\pi \times ratio} - 0.5
   * \f}
   * convert to integer
   * \f{eqnarray*}{
   * (int)s &=& 0.25d\sqrt{\pi \times ratio} - 0.5 + 0.5 \\
   *        &=& 0.25d\sqrt{\pi \times ratio}
   * \f}
   */
  const int bsize(2 * peakDiameter );
  _box = make_pair(s.value("BoxSizeX", bsize).toInt(),
                   s.value("BoxSizeY",bsize).toInt());

  /** min nbr of pixels should reflect the area under the bragg peak
   *  \f$ nbrPeaks = 0.25 \pi d \f$
   */
  _minNbrPixels = s.value("MinNbrPixels",0.25*3.14159*square(peakDiameter)).toInt();

  setupGeneral();
  bool ret (setupCondition());

  /** use fixed value for wavelength if value can be converted to double,
   *  otherwise use the wavelength from the processor
   */
  bool wlIsDouble(false);
  QString wlkey("Wavelength_A");
  QString wlparam(s.value(wlkey,"1").toString());
  double wlval(wlparam.toDouble(&wlIsDouble));
  if (wlIsDouble)
  {
    _wavelength = wlval;
    _getLambda = std::tr1::bind(&pp208::lambdaFromConstant,this,_1);
  }
  else
  {
    _wavelengthPP = setupDependency(wlkey.toStdString());
    ret = _wavelengthPP && ret;
    _getLambda = std::tr1::bind(&pp208::lambdaFromProcessor,this,_1);
  }
  /** use fixed value for detector distance if value can be converted to double,
   *  otherwise use the detector distance from the processor
   */
  bool ddIsDouble(false);
  QString ddkey("DetectorDistance_m");
  QString ddparam(s.value(ddkey,"60e-2").toString());
  double ddval(ddparam.toDouble(&ddIsDouble));
  if (ddIsDouble)
  {
    _detdist = ddval;
    _getDistance = std::tr1::bind(&pp208::distanceFromConstant,this,_1);
  }
  else
  {
    _detdistPP = setupDependency(ddkey.toStdString());
    ret = _detdistPP && ret;
    _getDistance = std::tr1::bind(&pp208::distanceFromProcessor,this,_1);
  }
  /** use fixed value for threshold if value can be converted to double,
   *  otherwise use the detector distance from the processor
   */
  bool threshIsDouble(false);
  QString threshkey("Threshold");
  QString threshparam(s.value(threshkey,"0").toString());
  const double threshval(threshparam.toDouble(&threshIsDouble));
  if (threshIsDouble)
  {
    _threshold = threshval;
    _thresh = std::tr1::bind(&pp208::thresholdFromConstant,this,_1);
  }
  else
  {
    _threshPP = setupDependency(threshkey.toStdString());
    ret = _threshPP && ret;
    _thresh = std::tr1::bind(&pp208::thresholdFromProcessor,this,_1);
  }

  _imagePP = setupDependency("ImageName");

  if (!(_imagePP && ret))
    return;

  /** check if the input processors have the correct type */
  if (_imagePP->result().dim() != 2)
    throw invalid_argument("pp208:loadSettings '" +name() +
                           "': input processor '" + _imagePP->name() +
                           "' is not a 2d result");
  if (!wlIsDouble && _wavelengthPP->result().dim() != 0)
    throw invalid_argument("pp208:loadSettings '" +name() +
                           "': wavelength processor '" + wlparam.toStdString() +
                           "' is not a 0d result");
  if (!ddIsDouble && _detdistPP->result().dim() != 0)
    throw invalid_argument("pp208:loadSettings '" +name() +
                           "': detector distance processor '" + ddparam.toStdString() +
                           "' is not a 0d result");
  if (!threshIsDouble && _threshPP->result().dim() != 0)
    throw invalid_argument("pp208:loadSettings '" +name() +
                           "': threshold processor '" + threshparam.toStdString()+
                           "' is not a 0d result");

  /** generate the lookuptable for the radia, if no geom file is provided,
   *  set the radia to 0
   */
  const string filename = s.value("GeometryFilename","wrong_file").toString().toStdString();
  const bool convertCheetahToCASSLayout = s.value("ConvertCheetahToCASSLayout",true).toBool();
  const double pixsizeX_m = s.value("PixelSizeX_m",109.92e-6).toDouble();
  const double pixsizeY_m = s.value("PixelSizeY_m",109.92e-6).toDouble();
  const result_t &srcImageHist(_imagePP->result());
  if (filename != "wrong_file")
  {
    GeometryInfo::conversion_t src2lab =
        GeometryInfo::generateConversionMap(filename,
                                            srcImageHist.size(),
                                            srcImageHist.shape().first,
                                            convertCheetahToCASSLayout);
    _src2labradius.resize(src2lab.size());
    for (size_t i=0; i < src2lab.size(); ++i)
    {
      const double x_m(src2lab[i].x *pixsizeX_m);
      const double y_m(src2lab[i].y *pixsizeY_m);
      _src2labradius[i] = sqrt(x_m*x_m + y_m*y_m);
    }
  }
  else
  {
    _src2labradius.resize(srcImageHist.size(),0);
  }

  /** set up the neighbouroffset list */
  _imageShape = srcImageHist.shape();
  _neighbourOffsets.clear();
//  _neighbourOffsets.push_back(+_imageShape.first-1);     //up left
  _neighbourOffsets.push_back(+_imageShape.first+0);     //up
//  _neighbourOffsets.push_back(+_imageShape.first+1);     //up right
  _neighbourOffsets.push_back(-1);                       //left
  _neighbourOffsets.push_back(+1);                       //right
//  _neighbourOffsets.push_back(-_imageShape.first-1);     //low left
  _neighbourOffsets.push_back(-_imageShape.first+0);     //low
//  _neighbourOffsets.push_back(-_imageShape.first+1);     //low right

  /** Create the result output */
  createHistList(result_t::shared_pointer(new result_t(nbrOf,0)));

  /** log what the user was requesting */
  string output("Processor '" + name() + "' finds bragg peaks." +
                "'. Boxsize '" + toString(_box.first)+"x"+ toString(_box.second)+
                "'. SectionSize '" + toString(_section.first)+"x"+ toString(_section.second)+
                "'. Threshold '" + threshparam.toStdString() +
                "'. MinSignalToNoiseRatio '" + toString(_minSnr) +
                "'. MinPixels '" + toString(_minNbrPixels) +
                "'. MinFraction '" + toString(_minRatio) +
                "'. Using input histogram :" + _imagePP->name() +
                "'. Wavelength :" + wlparam.toStdString() +
                "'. Detector Distance :" + ddparam.toStdString() +
                "'. Geomfile :" + filename +
                "'. Convert Cheetah to CASS :" + (convertCheetahToCASSLayout? "true":"false") +
                "'. Pixelsize [m] :" + toString(pixsizeX_m) + "x" + toString(pixsizeY_m) +
                "'. Condition is '" + _condition->name() + "'");
  Log::add(Log::INFO,output);
}

double pp208::lambdaFromProcessor(const CASSEvent::id_t& id)
{
  const result_t &wavelength(_wavelengthPP->result(id));
  QReadLocker lock(&wavelength.lock);
  return wavelength.getValue();
}

double pp208::distanceFromProcessor(const CASSEvent::id_t& id)
{
  const result_t &detdist(_detdistPP->result(id));
  QReadLocker lock(&detdist.lock);
  return detdist.getValue();
}

pp208::pixelval_t pp208::thresholdFromProcessor(const CASSEvent::id_t& id)
{
  const result_t &thresh(_threshPP->result(id));
  QReadLocker lock(&thresh.lock);
  return thresh.getValue();
}

pp208::pixelval_t pp208::thresholdFromConstant(const CASSEvent::id_t&)
{
  return _threshold;
}

int pp208::getBoxStatistics(result_t::const_iterator pixel,
                            const index_t linIdx, const shape_t &box, stat_t &stat)
{
  enum{use,skip};

  /** get coordinates of pixel from the linearized index */
  const index_t col(linIdx % _imageShape.first);
  const index_t row(linIdx / _imageShape.first);

  /** make sure that pixel is located such that the box will not conflict with
   *  the image and section boundaries. If it does continue with next pixel
   */
  if (col < box.first  || _imageShape.first - box.first  < col || //within box in x
      row < box.second || _imageShape.second - box.second < row || //within box in y
      (col - box.first)  / _section.first  != (col + box.first)  / _section.first || //x within same section
      (row - box.second) / _section.second != (row + box.second) / _section.second)  //y within same section
    return skip;

  /** go through all pixels defined by the box from -rows ... rows, -cols ... cols */
  for (shape_t::second_type bRow = -box.second; bRow <= box.second; ++bRow)
  {
    for (shape_t::first_type bCol = -box.first; bCol <= box.first; ++bCol)
    {
      /** check if the current box pixel value is bigger than the center pixel
       *  value. If so skip this pixel
       */
      const index_t bLocIdx(bRow*_imageShape.first+bCol);
      const pixelval_t bPixel(pixel[bLocIdx]);
      if(*pixel < bPixel )
        return skip;

      /** if its not a bad pixel add pixel to statistics */
      if (!qFuzzyIsNull(bPixel))
        stat.addDatum(bPixel);
    }
  }
  return use;
}

int pp208::isNotHighest(result_t::const_iterator pixel,
                        const index_t linIdx, shape_t box, stat_t &stat)
{
  enum{use,skip};

  bool boxsizeincreased(false);
  do
  {
    /** check whether current pixel value is highest and generate background values */
    stat.reset();
    if (getBoxStatistics(pixel,linIdx,box,stat))
      return skip;

    /** increase the box size and start over if the fraction of outliers to
     *  points used in the statistics is smaller than requested.
     */
    if (stat.nbrPointsUsed() < _minRatio * stat.nbrUpperOutliers())
    {
      ++(box.first);
      ++(box.second);
      boxsizeincreased = true;
    }
    else
    {
      boxsizeincreased = false;
    }
  }
  while(boxsizeincreased);

  /** skip this pixel if there are not enough pixels that could potentially be
   *  part of the bragg peak.
   */
  if (stat.nbrUpperOutliers() < _minNbrPixels)
    return skip;

  return use;
}


void pp208::process(const CASSEvent & evt, result_t &result)
{
  const result_t &image(_imagePP->result(evt.id()));
  QReadLocker hLock(&image.lock);

  /** clear the resulting table to fill it with the values of this image */
  result.resetTable();

  /** get a table row that we can later add to the table */
  table_t peak(nbrOf,0);

  /** set up a mask that we can see which pixels have been treated */
  vector<bool> checkedPixels(image.size(),false);

  /** get iterators for the mask and the image with which we can iterate through
   *  the image, also rember which linearized index we are working on right
   *  now, to be able to retrieve the column and row that the current pixel
   *  corresponsed to.
   */
  const pixelval_t thresh(_thresh(evt.id()));
  vector<bool>::iterator checkedPixel(checkedPixels.begin());
  result_t::const_iterator pixel(image.begin());
  result_t::const_iterator ImageEnd(image.end());
  index_t idx(0);
  for (;pixel != ImageEnd; ++pixel,++idx, ++checkedPixel)
  {
    /** check if pixel should be treated, when it has been treated before
     *  continue with the next pixel
     */
    if (*checkedPixel)
      continue;

    /** check if pixel is above the threshold, otherwise continue with the next
     *  pixel
     */
    if (*pixel < thresh)
      continue;

    /** check if pixel is highest within box. If so, check if there are enough
     *  pixels that are not outliers. If there are not enough pixels in the box
     *  increase the box size and do the checks all over again. If the pixel is
     *  not highest within the bigger box continue with the next pixel.
     */
    shape_t box(_box);
    stat_t stat(_minSnr);
    if (isNotHighest(pixel, idx, box, stat))
      continue;

    /** retrive statistical values from the calculator */
    const stat_t::value_type mean(stat.mean());
    const stat_t::value_type stdv(stat.stdv());

    /** get coordinates of pixel in image from the linearized index */
    const index_t col(idx % _imageShape.first);
    const index_t row(idx / _imageShape.first);

    /** from the central pixel look around and see which pixels are direct
     *  part of the peak. If a neighbour is found, mask it such that one does
     *  not use it twice.
     *
     *  Create a list that should contain the indizes of the pixels that are
     *  part of the peak. Go through that list and check for each found
     *  neighbour whether it also has a neighbour. If so add it to the list, but
     *  only if it is part of the box.
     */
    vector<index_t> peakIdxs;
    peakIdxs.push_back(idx);
    *checkedPixel = true;
    for (size_t pix=0; pix < peakIdxs.size(); ++pix)
    {
      const index_t pixpos(peakIdxs[pix]);
      neighbourList_t::const_iterator ngbrOffset(_neighbourOffsets.begin());
      neighbourList_t::const_iterator neighboursEnd(_neighbourOffsets.end());
      while(ngbrOffset != neighboursEnd)
      {
        const index_t ngbrIdx(pixpos + *ngbrOffset++);
        const index_t ngbrCol(ngbrIdx % _imageShape.first);
        const index_t ngbrRow(ngbrIdx / _imageShape.first);
        /** @note need to check if the pixel has been used after verifiying that
         *  the pixel is within the box. Otherwise in rare occasions it will
         *  happen that the pixelindex is not within the boolean array which
         *  will cause a segfault.
         */
        if (box.first < abs(col - ngbrCol) ||   //pixel not inside box col
            box.second < abs(row - ngbrRow))    //pixel not inside box row
          continue;
        if (checkedPixels[ngbrIdx])  //pixel has been used
          continue;
        const pixelval_t ngbrPixel(image[ngbrIdx]);
        const pixelval_t ngbrPixelWOBckgnd(ngbrPixel - mean);
        if ((_minSnr * stdv) < ngbrPixelWOBckgnd)
        {
          peakIdxs.push_back(ngbrIdx);
          checkedPixels[ngbrIdx] = true;
        }
      }
    }

    /** if the peak doesn't have enough pixels continue with next pixel */
    if (static_cast<int>(peakIdxs.size()) < _minNbrPixels)
      continue;

    /** go through all pixels in the box, find out which pixels are part of the
     *  peak and centroid them. Mask all pixels in the box as checked so one
     *  does not check them again
     */
    pixelval_t integral = 0;
    pixelval_t weightCol = 0;
    pixelval_t weightRow = 0;
    index_t max_radiussq=0;
    index_t min_radiussq=max(box.first,box.second)*max(box.first,box.second);
    for (int bRow = -box.second; bRow <= box.second; ++bRow)
    {
      for (int bCol = -box.first; bCol <= box.first; ++bCol)
      {
        const index_t bLocIdx(bRow*_imageShape.first+bCol);
        if (checkedPixel[bLocIdx])
        {
          const pixelval_t bPixel(pixel[bLocIdx]);
          const pixelval_t bPixelWOBckgnd(bPixel - mean);
          const index_t radiussq(bRow*bRow + bCol*bCol);
          if (radiussq > max_radiussq)
            max_radiussq = radiussq;
          if (radiussq < min_radiussq)
            min_radiussq = radiussq;
          integral += bPixelWOBckgnd;
          weightCol += (bPixelWOBckgnd * (bCol+col));
          weightRow += (bPixelWOBckgnd * (bRow+row));
        }
        checkedPixel[bLocIdx] = true;
      }
    }

    /** set the peak's properties and add peak to the list of found peaks */
    peak[centroidColumn] = weightCol / integral;
    peak[centroidRow] = weightRow / integral;
    peak[Intensity] = integral;
    peak[nbrOfPixels] = peakIdxs.size();
    peak[SignalToNoiseHighestPixel] = (*pixel-mean)/stdv;
    peak[SignalToNoiseSpot] = integral/stdv;
    peak[MaxRadius] = sqrt(max_radiussq);
    peak[MinRadius] = sqrt(min_radiussq);
    peak[Index] = idx;
    peak[Column] = col;
    peak[Row] = row;
    peak[MaxADU] = *pixel;
    peak[LocalBackground] = mean;
    peak[LocalBackgroundDeviation] = stdv;
    peak[nUpOutliers] = stat.nbrUpperOutliers();
    peak[Resolution] = 1.0 / ((2.0/_getLambda(evt.id())) *
                              sin(0.5*atan(_src2labradius[idx]/_getDistance(evt.id()))));

    result.appendRows(peak);
  }

  /** tell that only the result of one event (image) is present in the table */
}





//******************** processor 209: cluster pixels ************


pp209::pp209(const name_t &name)
  : Processor(name)
{
  loadSettings(0);
}

void pp209::loadSettings(size_t)
{
  CASSSettings s;

  s.beginGroup("Processor");
  s.beginGroup(QString::fromStdString(name()));

  _factor = s.value("Factor",1).toFloat();
  setupGeneral();
  bool ret (setupCondition());
  _imagePP = setupDependency("ImageName");
  ret = _imagePP && ret;
  _threshPP = setupDependency("Threshold");
  ret = _threshPP && ret;

  if (!ret)
    return;

  /** get references to the input results */
  const result_t &srcImageHist(_imagePP->result());
  const result_t &srcThreshHist(_threshPP->result());

  /** get the shape of the input results */
  _imageShape = srcImageHist.shape();
  shape_t threshShape = srcThreshHist.shape();


  /** check if the input processors have the correct type */
  if (srcImageHist.dim() != 2)
    throw invalid_argument("pp209:loadSettings '" +name() +
                           "': image processor '" + _imagePP->name() +
                           "' doesn't contain a 2d result");
  if (srcThreshHist.dim() != 2)
    throw invalid_argument("pp209:loadSettings '" +name() +
                           "': threshold processor '" + _threshPP->name() +
                           "' doesn't contain a 2d result");
  if (_imageShape.first != threshShape.first)
    throw invalid_argument("pp209:loadSettings '" +name() +
                           "': threshold processor '" + _threshPP->name() +
                           "' x-size '" + toString(threshShape.first) +
                           "' while image processor '" + _imagePP->name() +
                           "' has x-size '" + toString(_imageShape.first));
  if (_imageShape.second != threshShape.second)
    throw invalid_argument("pp209:loadSettings '" +name() +
                           "': threshold processor '" + _threshPP->name() +
                           "' y-size '" + toString(threshShape.second) +
                           "' while image processor '" + _imagePP->name() +
                           "' has y-size '" + toString(_imageShape.second));

  /** set up the neighbouroffset list */
  _neighbourOffsets.clear();
  _neighbourOffsets.push_back(+_imageShape.first+0);     //up
  _neighbourOffsets.push_back(-1);                       //left
  _neighbourOffsets.push_back(+1);                       //right

  /** Create the result output */
  createHistList(result_t::shared_pointer(new result_t(nbrOf,0)));

  /** log what the user was requesting */
  string output("Processor '" + name() + "' clusters pixels above threshold in '" +
                _imagePP->name() +
                "'. Threshold '" + _threshPP->name() +
                "'. Factor '" + toString(_factor) +
                "'. Condition is '" + _condition->name() + "'");
  Log::add(Log::INFO,output);
}

void pp209::process(const CASSEvent & evt, result_t &result)
{
  const result_t &image(_imagePP->result(evt.id()));
  QReadLocker imageLock(&image.lock);

  const result_t &thresholds(_threshPP->result(evt.id()));
  QReadLocker threshLock(&image.lock);

  /** clear the resulting table to fill it with the values of this image */
  result.resetTable();

  /** get a table row that we can later add to the table */
  table_t peak(nbrOf,0);

  /** set up a mask that we can see which pixels have been treated */
  vector<bool> checkedPixels(image.size(),false);

  /** get iterators for the mask and the image with which we can iterate through
   *  the image, also rember which linearized index we are working on right
   *  now, to be able to retrieve the column and row that the current pixel
   *  corresponsed to.
   */
  vector<bool>::iterator checkedPixel(checkedPixels.begin());
  result_t::const_iterator pixel(image.begin());
  result_t::const_iterator ImageEnd(image.begin() + image.datasize());
  result_t::const_iterator thresh(thresholds.begin());
  index_t idx(0);
  for (;pixel != ImageEnd; ++pixel, ++idx, ++checkedPixel, ++thresh)
  {
    /** check if pixel should be treated, when it has been treated before
     *  continue with the next pixel
     */
    if (*checkedPixel)
      continue;

    /** if it wasn't checked, then its checked now */
    *checkedPixel = true;

    /** check if pixel is above the pixelwise threshold, which is increased
     *  by the user provided factor, this allows to check the pixels using
     *  a noise map.
     */
    if (*pixel < (_factor * (*thresh)))
      continue;

    /** from the pixel look around and see which pixels are also above the
     *  threshold. If a neighbour is found, mask it such that one does
     *  not use it twice.
     *
     *  Create a list that should contain the indizes of the pixels that are
     *  part of the peak. Go through that list and check for each found
     *  neighbour whether it also has a neighbour. If so add it to the list, but
     *  only if its above the threshold
     */
    vector<index_t> pixIdxs;
    pixIdxs.push_back(idx);
    for (size_t pix(0); pix < pixIdxs.size(); ++pix)
    {
      const index_t pixpos(pixIdxs[pix]);
      neighbourList_t::const_iterator ngbrOffset(_neighbourOffsets.begin());
      neighbourList_t::const_iterator neighboursEnd(_neighbourOffsets.end());
      while(ngbrOffset != neighboursEnd)
      {
        const index_t ngbrIdx(pixpos + *ngbrOffset++);
        const index_t ngbrCol(ngbrIdx % _imageShape.first);
        const index_t ngbrRow(ngbrIdx / _imageShape.first);
        /** check if the neighbour is within the image shape */
        if ( ngbrCol < 0 ||                  // neighbour is to the left of image
             ngbrCol >= _imageShape.first || // neighbour is to the right of image
             ngbrRow < 0 ||                  // neighbour is below the image
             ngbrRow >= _imageShape.second)  // neighbour is above the image
          continue;
        /** check if neighbour was already treated */
        if (checkedPixels[ngbrIdx])
          continue;
        /** when the checks have passed it will be treated */
        checkedPixels[ngbrIdx] = true;
        /** check whether its above the threshold, if so add it to the pixellist */
        if ((_factor * thresholds[ngbrIdx])< image[ngbrIdx])
          pixIdxs.push_back(ngbrIdx);
      }
    }

    /** now that we found all the pixels that are connected to the starting
     *  pixel, go through all pixels in the list and centroid them.
     */
    pixelval_t integral(0);
    pixelval_t weightCol(0);
    pixelval_t weightRow(0);
    pixelval_t maxPixVal(0);
    index_t maxPixIdx(0);
    index_t maxPixCol(0);
    index_t maxPixRow(0);
    index_t minCol(0);
    index_t maxCol(_imageShape.first);
    index_t minRow(0);
    index_t maxRow(_imageShape.second);
    vector<index_t>::const_iterator pixIdx(pixIdxs.begin());
    vector<index_t>::const_iterator pixIdxEnd(pixIdxs.end());
    for (; pixIdx != pixIdxEnd; ++pixIdx)
    {
      /** get the pixel parameters */
      const index_t pixCol(*pixIdx % _imageShape.first);
      const index_t pixRow(*pixIdx / _imageShape.first);
      const pixelval_t pixVal(image[*pixIdx]);

      /** calculate the peak properties */
      integral += pixVal;
      weightCol += (pixVal * static_cast<pixelval_t>(pixCol));
      weightRow += (pixVal * static_cast<pixelval_t>(pixRow));

      /** find the maximum pixel */
      if (maxPixVal < pixVal)
      {
        maxPixVal = pixVal;
        maxPixIdx = *pixIdx;
        maxPixCol = pixCol;
        maxPixRow = pixRow;
      }

      /** find min and max col and row */
      if (maxCol < pixCol)
        maxCol = pixCol;
      if (pixCol < minCol)
        minCol = pixCol;
      if (maxRow < pixRow)
        maxRow = pixRow;
      if (pixRow < minRow)
        minRow = pixRow;
    }

    /** set the peak's properties and add peak to the list of found peaks */
    peak[Integral] = integral;
    peak[CentroidColumn] = weightCol / integral;
    peak[CentroidRow] = weightRow / integral;
    peak[MaxADU] = maxPixVal;
    peak[Index] = maxPixIdx;
    peak[Column] = maxPixCol;
    peak[Row] = maxPixRow;
    peak[MaxColumn] = maxCol;
    peak[MinColumn] = minCol;
    peak[ColumnSize] = maxCol - minCol;
    peak[MaxRow] = maxRow;
    peak[MinRow] = minRow;
    peak[RowSize] = maxRow - minRow;
    peak[NbrOfPixels] = pixIdxs.size();

    result.appendRows(peak);
  }

  /** tell that only the result of one event (image) is present in the table */
}





//******************** processor 210: Find Bragg peaks using MAD **************

pp210::pp210(const name_t &name)
  : Processor(name)
{
  loadSettings(0);
}

void pp210::loadSettings(size_t)
{
  CASSSettings s;

  s.beginGroup("Processor");
  s.beginGroup(QString::fromStdString(name()));

  _threshold = s.value("Threshold",10).toFloat();
  _minSnr = s.value("MinSignalToNoiseRatio",8).toFloat();
  _minNbrSnr = s.value("MinSignalToNoiseRatioForNeighbours",4).toFloat();
  _maxAdditionalPixels = s.value("MaxAdditionalPixels",2).toFloat();
  _minNbrPixels = s.value("MinNbrPixels",2).toInt();
  _badPixVal = s.value("BadPixelValue",0).toFloat();
  _box = make_pair(s.value("BoxSizeX", 5).toInt(),
                   s.value("BoxSizeY",5).toInt());
  _section = make_pair(s.value("SectionSizeX", 1024).toUInt(),
                       s.value("SectionSizeY",512).toUInt());

  setupGeneral();
  bool ret (setupCondition());

  _imagePP = setupDependency("ImageName");

  /** return if any of the dependecies has not been loaded */
  if (!(_imagePP && ret))
    return;

  const result_t &srcImageHist(_imagePP->result());
  _imageShape = srcImageHist.shape();

  _neighbourOffsets.clear();
//  _neighbourOffsets.push_back(+_imageShape.first-1);     //up left
  _neighbourOffsets.push_back(+_imageShape.first+0);     //up
//  _neighbourOffsets.push_back(+_imageShape.first+1);     //up right
  _neighbourOffsets.push_back(-1);                       //left
  _neighbourOffsets.push_back(+1);                       //right
//  _neighbourOffsets.push_back(-_imageShape.first-1);     //low left
  _neighbourOffsets.push_back(-_imageShape.first+0);     //low
//  _neighbourOffsets.push_back(-_imageShape.first+1);     //low right

  /** Create the result output */
  createHistList(result_t::shared_pointer(new result_t(nbrOf,0)));

  /** log what the user was requesting */
  string output("Processor '" + name() + "' finds bragg peaks using mad." +
                ". Boxsize '" + to_string(_box.first)+"x"+ to_string(_box.second)+
                "'. SectionSize '" + to_string(_section.first)+"x"+ to_string(_section.second)+
                "'. Threshold '" + to_string(_threshold) +
                "'. MinSignalToNoiseRatio '" + to_string(_minSnr) +
                "'. MinPixels '" + to_string(_minNbrPixels) +
                "'. Condition is '" + _condition->name() + "'");
  Log::add(Log::INFO,output);
}

void pp210::process(const CASSEvent & evt, result_t &result)
{
  const result_t &image(_imagePP->result(evt.id()));
  QReadLocker hLock(&image.lock);

  /** clear the resulting table to fill it with the values of this image */
  result.resetTable();

  /** get a table row that we can later add to the table */
  table_t peak(nbrOf,0);

  /** set up a mask that we can see which pixels have been treated */
  vector<bool> checkedPixels(image.size(),false);
  vector<bool>::iterator checkedPixel(checkedPixels.begin());

  /** get iterators to the beginning and end of the image */
  result_t::const_iterator pixel(image.begin());
  result_t::const_iterator ImageEnd(image.end());
  index_t idx(0);
  for (;pixel != ImageEnd; ++pixel,++idx, ++checkedPixel)
  {
    /** check if pixel should be treated, when it has been treated before
     *  continue with the next pixel
     */
    if (*checkedPixel)
      continue;

    /** check if pixel is above the threshold, otherwise continue with the next
     *  pixel
     */
    if (*pixel < _threshold)
      continue;

    /** caluclate where we are within the image */
    const index_t col(idx % _imageShape.first);
    const index_t row(idx / _imageShape.first);

    /** get the right box boundaries, so that we're within the image and the
     *  sections
     */
    const shape_t::first_type  colLow (col - _box.first);
    const shape_t::first_type  colHigh(col + _box.first);
    const shape_t::second_type rowLow (row - _box.second);
    const shape_t::second_type rowHigh(row + _box.second);
    if (colLow < 0 || _imageShape.first  < colHigh ||
        rowLow < 0 || _imageShape.second < rowHigh)
      continue;

    /** calc the section that the box edges are in */
    const int colLowSection (colLow  / _section.first);
    const int colHighSection(colHigh / _section.first);
    const int rowLowSection (rowLow  / _section.second);
    const int rowHighSection(rowHigh / _section.second);
    if (colLowSection != colHighSection ||
        rowLowSection != rowHighSection)
      continue;

    /** check if pixel is highest within box. If so, check if there are enough
     *  pixels that are not outliers. If there are not enough pixels in the box
     *  increase the box size and do the checks all over again. If the pixel is
     *  not highest within the bigger box continue with the next pixel.
     */
    /** now check if its the highest pixel */
    vector<float> medCont;
    bool isHighest(true);
    for (shape_t::second_type bRow(rowLow);(isHighest) && (bRow < rowHigh); ++bRow)
    {
      for (shape_t::first_type bCol(colLow); (isHighest) && (bCol < colHigh); ++bCol)
      {
        /** check if we're still inside the image */
        const float bPixel(image[bCol + bRow*_imageShape.first]);
        if ( *pixel < bPixel)
        {
          isHighest = false;
          break;
        }
        medCont.push_back(bPixel);
      }
    }
    if (!isHighest)
      continue;

    /** calculate the MAD */
    const size_t medianpos(0.5*medCont.size());
    const auto target(medCont.begin() + medianpos);
    std::nth_element(medCont.begin(), target, medCont.end());
    const float  median(*target);
    for (auto &val : medCont) val = fabs(val-median);
    std::nth_element(medCont.begin(), target, medCont.end());
    const float  mad((*target)*1.4862);

    /** check if the median is postive */
    if ( median < 0 )
      continue;

    /** check if pixel satifies the snr requirement */
    const float snr((*pixel - median) / mad);
    if (snr < _minSnr)
      continue;

//    /** retrive statistical values from the calculator */
//    const stat_t::value_type mean(stat.mean());
//    const stat_t::value_type stdv(stat.stdv());

    /** from the central pixel look around and see which pixels are direct
     *  part of the peak. If a neighbour is found, mask it such that one does
     *  not use it twice.
     *
     *  Create a list that should contain the indizes of the pixels that are
     *  part of the peak. Go through that list and check for each found
     *  neighbour whether it also has a neighbour. If so add it to the list, but
     *  only if it is part of the box.
     */
    vector<index_t> peakIdxs;
    peakIdxs.push_back(idx);
    *checkedPixel = true;
    bool isGoodSpot(true);
    for (size_t pix=0; isGoodSpot && (pix < peakIdxs.size()); ++pix)
    {
      const index_t pixpos(peakIdxs[pix]);
      neighbourList_t::const_iterator ngbrOffset(_neighbourOffsets.begin());
      neighbourList_t::const_iterator neighboursEnd(_neighbourOffsets.end());
      while(ngbrOffset != neighboursEnd)
      {
        const index_t ngbrIdx(pixpos + *ngbrOffset++);
        const index_t ngbrCol(ngbrIdx % _imageShape.first);
        const index_t ngbrRow(ngbrIdx / _imageShape.first);
        /** @note need to check if the pixel has been used after verifiying that
         *  the pixel is within the box. Otherwise in rare occasions it will
         *  happen that the pixelindex is not within the boolean array which
         *  will cause a segfault.
         */
        if (_box.first < abs(col - ngbrCol) ||   //pixel not inside box col
            _box.second < abs(row - ngbrRow))    //pixel not inside box row
          continue;
        if (checkedPixels[ngbrIdx])  //pixel has been used
          continue;
        checkedPixels[ngbrIdx] = true;
        const pixelval_t ngbrPixel(image[ngbrIdx]);
        /** check if the potential neighbour is a bad pixel, if so ignore this
         *   bragg spot
         */
        if (abs(ngbrPixel - _badPixVal) < std::numeric_limits<pixelval_t>::epsilon())
        {
          isGoodSpot = false;
          break;
        }
        const pixelval_t ngbrPixelWOBckgnd(ngbrPixel - median);
        if ((_minNbrSnr * mad) < ngbrPixelWOBckgnd)
        {
          peakIdxs.push_back(ngbrIdx);
        }
      }
    }

    if (!isGoodSpot)
      continue;

    /** if the peak doesn't have enough pixels continue with next pixel */
    if (static_cast<int>(peakIdxs.size()) < _minNbrPixels)
      continue;

    /** go through all pixels in the box, find out which pixels are part of the
     *  peak and centroid them. Mask all pixels in the box as checked so one
     *  does not check them again
     */
    pixelval_t integral = 0;
    pixelval_t weightCol = 0;
    pixelval_t weightRow = 0;
    index_t max_radiussq=0;
    index_t min_radiussq=max(_box.first,_box.second)*max(_box.first,_box.second);
    size_t badcount(0);
    for (int bRow(-_box.second); bRow <= _box.second; ++bRow)
    {
      for (int bCol(-_box.first); bCol <= _box.first; ++bCol)
      {
        const index_t bLocIdx(bRow*_imageShape.first+bCol);
        const pixelval_t bPixel(pixel[bLocIdx]);
        const pixelval_t bPixelWOBckgnd(bPixel - median);
        if (checkedPixel[bLocIdx])
        {
          const index_t radiussq(bRow*bRow + bCol*bCol);
          if (radiussq > max_radiussq)
            max_radiussq = radiussq;
          if (radiussq < min_radiussq)
            min_radiussq = radiussq;
          integral += bPixelWOBckgnd;
          weightCol += (bPixelWOBckgnd * (bCol+col));
          weightRow += (bPixelWOBckgnd * (bRow+row));
        }
        /** count how many pixels are not connected to the center one but still
         * above the snr
         */
        else if ((_minNbrSnr * mad) < bPixelWOBckgnd)
        {
          ++badcount;
        }
        checkedPixel[bLocIdx] = true;
      }
    }

    /** check how many pixels that are not part of the main peak there are */
    if (_maxAdditionalPixels < badcount)
      continue;


    /** set the peak's properties and add peak to the list of found peaks */
    peak[centroidColumn] = weightCol / integral;
    peak[centroidRow] = weightRow / integral;
    peak[Intensity] = integral;
    peak[nbrOfPixels] = peakIdxs.size();
    peak[SignalToNoiseHighestPixel] = snr;
    peak[SignalToNoiseSpot] = integral/mad;
    peak[MaxRadius] = sqrt(max_radiussq);
    peak[MinRadius] = sqrt(min_radiussq);
    peak[Index] = idx;
    peak[Column] = col;
    peak[Row] = row;
    peak[MaxADU] = *pixel;
    peak[LocalBackground] = median;
    peak[LocalBackgroundDeviation] = mad;
    peak[AdditionalPixels] = badcount;

    result.appendRows(peak);
  }
}