stdalgorithm.hpp

Go to the documentation of this file.

#ifndef STDALGORITHM_HPP
#define STDALGORITHM_HPP

#include "imgclass.hpp"
#include "kube-gustavson-fft.hpp"
#include "matrixutil.hpp"
#include "minmaxheap.hpp"
#include <vector>
#include <cmath>

/* contains algorithms not in imgclass, but very useful.
   I don't think it is adviseable to have functions that 
   merely operate on a const bildfloat as member functions in bildfloat,
   like for example out-place FFT or searching for extrema */

class Parabelzentrum {
  // Funktor, der die Mitte einer Punktmenge=Oberfläche über
  // den Parabelfit bestimmt
  const bildfloat& bild;
  int pktzahl;
  int pktmax;
  double **R;
  double *f;
  
  static const int order; 
  double x_cent, y_cent;
  double x_estimate, y_estimate;
  double s_xcent, s_ycent;
  minheap <double> xmax, ymax;
  maxheap <double> xmin, ymin;
public:
  double *koeff;
  double *errors;
  double *rel_errors;
  Parabelzentrum(const bildfloat& bild_, int pktmax_);
  ~Parabelzentrum();
  void estimate(double x, double y);
  void operator() (int x, int y);
  void clear();
  bool fit();
  double x();
  double y();
  double height();

  // statistische Fehler
  double s_x() {return s_xcent; }
  double s_y() { return s_ycent; }
  double s_height();
};

inline double Parabelzentrum::x() { return x_cent+x_estimate; }

inline double Parabelzentrum::y() { return y_cent+y_estimate; }

inline double Parabelzentrum::height() { 
  return koeff[0]-SQR(koeff[1])/(4*koeff[3])-SQR(koeff[2])/(4*koeff[4]); 
}

inline double Parabelzentrum::s_height() {
  double errsum=SQR(errors[0]);
  errsum+=SQR(SQR(koeff[1])/(4*koeff[3]))*(2*SQR(rel_errors[1])+SQR(rel_errors[3]));
  errsum+=SQR(SQR(koeff[2])/(4*koeff[4]))*(2*SQR(rel_errors[2])+SQR(rel_errors[4]));
  return sqrt(errsum);
}  
  

/* call the functor for every point filled */
template <typename TBild, typename TFunctor, typename TBorder>
void execute_hubbel (const TBild &hubbelbild, region_t &region, int x, int y, const TBorder &border, TFunctor &fun) {
  
  if (region(y,x)) return;
  int xmin=x; int xmax=x;
  while (!border(hubbelbild, region, xmin-1, y)) xmin--;
  while (!border(hubbelbild, region, xmax+1, y)) xmax++;
  for (int x=xmin; x<=xmax; x++) {
    region(y,x)=true;
    fun(x,y);
  }  
  
  
  for (int x=xmin; x<=xmax; x++) {
    if (!border(hubbelbild, region, x,y-1)) execute_hubbel(hubbelbild, region, x,y-1, border, fun);
    if (!border(hubbelbild, region, x,y+1)) execute_hubbel(hubbelbild, region, x,y+1, border, fun);
  }
}  

struct minborder {
  // minimum height for inner region
  float threshold;
  minborder (float thresh) : threshold(thresh) {}
  
  bool operator ()  (const bildfloat& hub, const region_t& region, int x, int y) const {
    if ((x<1) || (y<1) || (x>=region.SIZEX-1) || (y>=region.SIZEY-1)) return true;
    // point outside the field or near border
    if (region(y,x)) return true;
    // already filled
    if (hub(y,x)<threshold) return true;
    // smaller than minimum height
    return false;
  }  
};  

struct maxborder {
  // maximum height for inner region
  float threshold;
  maxborder (float thresh) : threshold(thresh) {}
  
  bool operator ()  (const bildfloat& hub, const region_t& region, int x, int y) const {
    if ((x<1) || (y<1) || (x>=region.SIZEX-1) || (y>=region.SIZEY-1)) return true;
    // point outside the field or near border
    if (region(y,x)) return true;
    // already filled
    if (hub(y,x)>threshold) return true;
    // greater than maximum height
    return false;
  }  
}; 

struct valleyborder {
  bool operator () (const bildfloat& hub, const region_t& region, int x, int y) const {
    if ((x<1) || (y<1) || (x>=region.SIZEX-1) || (y>=region.SIZEY-1)) return true;
    // point outside the field or near border
    if (region(y,x)) return true;
    // already filled
  #define LOCALMIN(y1,x1,y2,x2) if ((hub(y,x)<=hub(y1,x1)) && (hub(y,x) <= hub(y2,x2))) return true 
    LOCALMIN(y,x-1, y,x+1);
    LOCALMIN(y+1,x,y-1,x); 
    LOCALMIN(y+1,x+1, y-1,x-1);
    LOCALMIN(y-1,x+1, y+1, x-1);
    return false;
  }  
};

template <typename TFunctor>
void execute_region (const region_t &region, TFunctor &fun) {
  MAP(region)
    if (region(y,x)) fun(x,y);
}


struct countfunc {
  int sum;
  countfunc() : sum(0) {}
  void operator () (int x, int y) {sum++; }
  int operator() () const {return sum; }
};  

struct hydrosum {
  bildfloat &heightfield;
  float sum;
  hydrosum(bildfloat& rb) : heightfield(rb), sum(0.0) {}
  void operator () (int x, int y) { sum+= SQR(heightfield(y,x)); }
  float operator() () const {return sum; }
};  

struct azimuth_average {
  int anz;
  std::vector<double> rsumme;
  std::vector<int> ranzahl;
  bildfloat &heightfield;
  float x_centr, y_centr;
  
  azimuth_average(bildfloat& rb, float x_centr_, float y_centr_, int anz_) : 
    anz(anz_), rsumme(anz_), ranzahl(anz_), heightfield(rb), x_centr(x_centr_), y_centr(y_centr_) 
  {
         for (int i=0; i<anz; i++) {
           ranzahl[i]=0;
           rsumme[i]=0.0;
         }  
   }
    
  void operator () (int x, int y)
  { 
     int dist=static_cast<int>(sqrt(SQR(x-x_centr)+SQR(y-y_centr)));
     if (dist<anz) {
       ranzahl[dist]++;
       rsumme[dist]+=heightfield(y,x);
     }  
  }
  float operator() (int dist) const {
    if ((dist<anz) && (ranzahl[dist]!=0)) return rsumme[dist]/ranzahl[dist];
    return nanf(""); // falls die Entfernung nicht aufgetaucht ist
  }
}; 


template <typename T>
struct minmaxpos {
  T val;
  int x, y;
  minmaxpos(const T& v, int x_, int y_) : val(v), x(x_), y(y_) {}
  minmaxpos() : val (0), x(0), y(0) {}
  bool operator > (const minmaxpos& v2) const { return val > v2.val; }
  bool operator < (const minmaxpos& v2) const { return val < v2.val; }
};

template <typename T>
struct getmax : public minheap< minmaxpos <T> > {
  const bildfloat& rb;
  getmax(const bildfloat& bfl, int n=1) : minheap< minmaxpos <T> > (n), rb(bfl) {}
  void  operator () (int x, int y) { put(minmaxpos<T>(rb(y,x),x,y)); }
  const minmaxpos<T>  operator () () const { return minheap < minmaxpos <T> >::operator() (); }
};

template <typename T>
struct getmin : public maxheap< minmaxpos <T> > {
  const bildfloat& rb;
  getmin(const bildfloat& bfl, int n=1) : maxheap< minmaxpos <T> > (n), rb(bfl) {}
  void  operator () (int x, int y) { put(minmaxpos<T>(rb(y,x),x,y)); }
  const minmaxpos<T>  operator () () const { return maxheap < minmaxpos <T> >::operator() (); }
};


typedef imgbase<512, 512, COMPLEX> bildcomplex_base;

class bildcomplex : public bildcomplex_base {
public:
  bildcomplex() : bildcomplex_base() { }
  bildcomplex(const bildfloat &realvalue);
  bildcomplex(const bildfloat& real, const bildfloat& imag);
  void decompose(bildfloat& real, bildfloat& imag);
  
  bildcomplex& FFT_forward();
  bildcomplex& FFT_backward();
  bildcomplex& Hamming_window(float width=220, float xm=256, float ym=256);
  bildcomplex& Hamming_unwindow(float width=220, float xm=256, float ym=256);


  bildfloat& norm(bildfloat& ret);
};  

typedef std::pair<float, float> Point;

typedef std::vector<Point> cell;

struct vec2d {
  double x, y;

  vec2d(double x_, double y_): x(x_), y(y_) {}
  vec2d operator - (const vec2d & subtr) const {
    return vec2d(x-subtr.x, y-subtr.y);
  }

  vec2d operator + (const vec2d & summand) {
    return vec2d(x+summand.x, y+summand.y);
  }

  vec2d& operator /= (double div) {
    x/=div;
    y/=div;
    return *this;
  }

  double operator * (const vec2d& v2) {
    // dot product
    return x*v2.x+y*v2.y;
  }

  vec2d operator * (double s) {
    // S-Multiplikation
    return vec2d(x*s, y*s);
  }
  
  vec2d operator / (double s) {
    // S-Multiplikation
    return vec2d(x/s, y/s);
  }  
};

inline double abs(const vec2d& v) {
  return sqrt(v.x*v.x+v.y*v.y);
}  


template <typename TFunctor>
void line_iterate(const vec2d& P1, const vec2d& P2, TFunctor& fun) 
{
  // call a functor for each point on a straight line between P1 and P2
  // could be optimized by Bresenham algorithm
  // however, premature optimization is the root of all evil
  
  vec2d r=P2-P1;
  // r ist Richtungsvektor der Geraden
  
  int stepend, stepinc;
  
  // r ist Richtungsvektor der gesuchten Geraden
  //cerr<<"Richtungsvektor "<<r.x<<", "<<r.y<<endl;
  if (abs(r.x) >abs(r.y)) {
     stepend=abs(r.x);
     r/=abs(r.x);
  } else {
     stepend=abs(r.y);
     r/=abs(r.y);
  }

  // Normierung, so dass immer ein Schritt entlang x oder y gemacht wird 

  for (int step=0; step<=stepend; step++) {
    fun(floor(r.x*step+P1.x+0.5), floor(r.y*step+P1.y+0.5));
  }  
}

struct mean_avg {
  const bildfloat& rb;
  double sum;
  double sqrsum;
  int pktzahl;
  mean_avg(const bildfloat& bild) : rb(bild), sum(0.0), sqrsum(0.0), pktzahl(0) {}
  void operator () (int x, int y) {
    // cerr<<".bild.display create rectangle "<<x<<" "<<y<<" "<<x<<" "<<y<<" -fill green -outline green -tag bound"<<endl;
    sum+=rb(y,x);
    sqrsum+=SQR(rb(y,x));
    pktzahl++;
  }
};  

struct surface {
  double surf;
  const bildfloat& rb;
  surface(const bildfloat& rb_) : surf(0.0), rb(rb_) {}
  void operator () (int x, int y) {
    double dx=(rb(y,   x+1)-rb(y,x-1))/(2*rb.pix_size);
    double dy=(rb(y+1, x) - rb(y-1, x))/(2*rb.pix_size);
    surf+=sqrt(1.0+SQR(dx)+SQR(dy));
  }  
  double operator() () { return surf; }
};


bool pnpoly(int npol, float *xp, float *yp, float x, float y);

class Voronoischnitt {
   const bildfloat& rb;
   float xmitte, ymitte;
   cell zelle;
   float *xp, *yp;
   int npol;
   vec2d r;
 public:
   Voronoischnitt(bildfloat &bfl, float xmitte_, float ymitte_, const char *vorofile);
   ~Voronoischnitt();
   void compute_schnitt(int nr);
   float boundary_average();
   double stddev;
   bool contains(float x, float y);
   vec2d direction();
   float mindist, maxdist;
};


#endif

---