imgclass.hpp

Go to the documentation of this file.

#ifndef imgclass_hpp
#define imgclass_hpp

#include <iostream>
#include <valarray>
#include <string>
#include <sstream>
#include "matrixutil.hpp"
#include "vc6compat.h"

#define TRACE(x) 
#define TRACEX(x,y)


#define STHROW(msg) { \
         std::ostringstream err;\
         err<<msg; \
         throw err.str(); }

#define THROWSPEC throw(std::string)



#define MAP(img) for (int y=0; y<img.SIZEY; y++) \
                  for (int x=0; x<img.SIZEX; x++)


#define MAP_CIRCLE(img, radius) MAP(img) \
                     if (SQR(x-img.SIZEX/2)+ SQR(y-img.SIZEY/2) <= SQR(radius))


#define MAP_REGION(region) MAP(region) if (region(y,x)) 


#define MAP_CIRCLE_XY(img, radius, xm, ym) \
               for (int y=ym-radius; y<=ym+radius; y++) \
                  for (int x=xm-radius; x<=xm+radius; x++) \
                     if ((x>=0) && (x<img.SIZEX) && (y>=0) && (y<img.SIZEY)) \
                      if (SQR(x-xm)+ SQR(y-ym) <= SQR(radius))



#define MAPF(img, Functor) MAP(img) { Functor(img(y,x)); }


#define MAPFXY(img, Functor) MAP(img) { Functor(x, y, img(y,x)); }

/* bequemes Iterieren über ein Bild */


 template <typename T> inline T SQR (T x) {
    return x * x;
  } 
 

  template < int s_x, int s_y, typename T > 
  class imgbase {
  public:

    typedef T data_t;

    const int SIZEX;

    const int SIZEY;
    
    imgbase ();  
    imgbase (const imgbase& im);
    ~imgbase (); 

    T& operator  () (int y, int x) { return bild[y][x]; }
    const T & operator  () (int y, int x) const { return bild[y][x]; }

    imgbase<s_x, s_y, T>& operator =  (const imgbase<s_x, s_y, T>& what); 
    imgbase<s_x, s_y, T>& operator /= (const imgbase<s_x, s_y, T>& what); 
    imgbase<s_x, s_y, T>& operator /= (const T& what);
    imgbase<s_x, s_y, T>& operator *= (const imgbase<s_x, s_y, T>& what); 
    imgbase<s_x, s_y, T>& operator *= (const T& what);
    imgbase<s_x, s_y, T>& operator -= (const imgbase<s_x, s_y, T>& what);
    imgbase<s_x, s_y, T>& operator -= (const T& what);
    imgbase<s_x, s_y, T>& operator += (const imgbase<s_x, s_y, T>& what);
    imgbase<s_x, s_y, T>& operator += (const T& what);
  protected:
    data_t ** bild; 
  };



  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>::imgbase () : 
     SIZEX (s_x), SIZEY (s_y) {
    // reserve memory for bild
     reserve_2D(bild, s_x, s_y);  
  }



  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>::imgbase (const imgbase &im) : 
     SIZEX (s_x), SIZEY (s_y) {
    // reserve memory for bild
    reserve_2D(bild, s_x, s_y);  
    operator = (im);  
   }
  

  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>::~imgbase () {
    // release memory
    delete_2D(bild);
  }
  
  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator = (const imgbase<s_x, s_y, T>& what)
  {
     MAP(what)
        bild[y][x]=what(y,x);
    return *this;       
  }     
  

  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator /= (const imgbase<s_x, s_y, T>& what)
  {
    MAP(what) {
        T val=what(y,x);
        if (val!=0) bild[y][x]/=val;
      } 
    return *this;       
  }     
  
  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator *= (const imgbase<s_x, s_y, T>& what)
  {
    MAP(what) {
        T val=what(y,x);
        bild[y][x]*=val;
      } 
    return *this;       
  }     


  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator -= (const imgbase<s_x, s_y, T>& what)
  {
    MAP((*this)) {
        T val=what(y,x);
        bild[y][x]-=val;
      } 
    return *this;       
  }     

  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator /= (const T& what)
  {
    MAP((*this)) {
        if (what!=0) bild[y][x]/=what;
      } 
    return *this;       
  }     
   
  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator *= (const T& what)
  {
    MAP((*this)) {
        bild[y][x]*=what;
      } 
    return *this;       
  }

  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator -= (const T& what)
  {
    MAP((*this)) {
        bild[y][x]-=what;
      } 
    return *this;       
  }     


  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator += (const imgbase<s_x, s_y, T>& what)
  {
     MAP(what) {
        T val=what(y,x);
        bild[y][x]+=val;
      } 
    return *this;       
  }     
 
  template < int s_x, int s_y, typename T >
  imgbase<s_x, s_y, T>& imgbase<s_x, s_y, T>::operator += (const T& what)
  {
     MAP((*this)) {
        bild[y][x]+=what;
      } 
    return *this;       
  }     

  
  template < int s_x, int s_y, typename T >
  T average(const imgbase<s_x, s_y, T> &img) {
     T temp=0;
     for (int y=0; y< s_y; y++)
       for (int x=0; x< s_x; x++)
         temp += img(y,x);
     return temp/(s_x*s_y);      
  }   
    

  typedef imgbase<512,512,unsigned short> bild16_base;
  
  class bildfloat; // forward declaration
                   // because of dependencies 
  

  class bild16 : public bild16_base {
    static const std::streamoff HIS_START;  // 0-100 header des His-files
  public:
    static const data_t intensity_max;
    static const data_t intensity_min;

  #pragma pack(push,1)          // stellt sicher, dass structur groesse auf 1 byte genau bestimmt wird
    typedef struct {
      unsigned short FileID;    //1
      unsigned short HeaderSize;
      unsigned short HeaderVersion;
      unsigned long FileSize;
      unsigned short ImageHeaderSize;   //5
      unsigned short ULX;
      unsigned short ULY;
      unsigned short BRX;
      unsigned short BRY;
      unsigned short NrOfFrames;        //10
      unsigned short Correction;
      double IntegrationTime;
      unsigned short TypeOfNumbers;     //13
    }  header_t;
  
  #pragma pack(pop) // stellt alten Wert wieder her
      
      bild16(const char *fname=NULL) THROWSPEC;
      bild16(unsigned short *data) THROWSPEC;
      ~bild16();
      header_t* file_header;
      void write_pgm16 (const char *name) const;
      void write_pgm8 (const char *name) const;
      void write_pgm8_bounds (const char *name, bild16::data_t low, bild16::data_t high) const;
      bild16& read_his (const char *name) THROWSPEC;
      void write_his (const char *name) const;
      bild16& read_pgm(const char *name) THROWSPEC;
      bild16& og_correct(const bild16& offset, const bild16& gain);
      bild16& drawx(int x, int y1, int y2) THROWSPEC;
      bild16& drawy(int y, int x1, int x2) THROWSPEC;
      bild16& kreuz(int x, int y) THROWSPEC;
      void fromfloat(const bildfloat &bf, float max=1.0);
      bild16& operator = (const bild16& b) {
        bild16_base::operator = (b);
        *file_header=*(b.file_header);
        return *this;
      } 
  };


   class rawfloat {
   public: 
     rawfloat (const float& f): f_ (f) {}
   private:
     static float dummy;
     const float& f_;
     friend std::ostream& operator << (std::ostream& os, const rawfloat& r) {
       
       return os.write (reinterpret_cast <const char *>(&(r.f_)), sizeof (r.f_));
     }
  };
  

  class irawfloat {
   public: 
     irawfloat (float& f) : f_(f) {}
   private:
     float& f_;
     
     friend std::istream& operator >> (std::istream& is, const irawfloat& r) { 
       return is.read(reinterpret_cast<char *>(&(r.f_)), sizeof(r.f_)); 
     }
  };
  

  class eichkurve {
    public:
      eichkurve();
      virtual ~eichkurve();
      double operator() (bild16::data_t val) const;
      virtual double eichfun(double x) const=0;
      virtual double deichfun_dx(double x) const=0;

      virtual double initguess(bild16::data_t val) const=0;
      void table_write(const char *name) const;
    protected:
      static const char magic[];
      double *lookup;
      double fsolve(bild16::data_t val, double hmin, double hmax) const;
      void   table_init(double hmin, double hmax);
      // Newton-Raphson-Verfahren
  };
  
  inline double eichkurve::operator() (bild16::data_t val) const 
  {
    return lookup[val];
  }
  

  class expkurve : public eichkurve {
    public:
      expkurve (double sA, double sg1, double sg2, double sg3,
                // Abschwächungskoeffizienten
                double sf1, double sf2
                // Mischparameter
              );
      double eichfun(double x) const;
      double deichfun_dx(double x) const;
      double initguess(bild16::data_t val) const;
              
    private:
      double A, g1, g2, g3, f1, f2;
  };



  class exp4kurve : public eichkurve {
    public:
      exp4kurve (double sA, double sa1, double sa2, double sa3,
                // rel. Anteile
                double sb1, double sb2, double sb3, double sb4
                // charakteristische Tiefen
              );
      double eichfun(double x) const;
      double deichfun_dx(double x) const;
      double initguess(bild16::data_t val) const;
              
    private:
      double A, a1, a2, a3, b1, b2, b3, b4;
  };

  

  class tablekurve : public eichkurve {
    public:
      tablekurve (const char *fname) THROWSPEC;
      double eichfun(double x) const;
      double deichfun_dx(double x) const;
      double initguess(bild16::data_t val) const;
  };


  class flatkurve : public eichkurve {
    public:
      flatkurve (bool normalize_=true);
      virtual double eichfun(double x) const;
      virtual double deichfun_dx(double x) const;
      virtual double initguess(bild16::data_t val) const;
    private:
      bool normalize;
  };


  typedef imgbase<512,512,bool> region_base;

  typedef region_base region_t; //später vielleicht abgeleitet

  class QRfit;
  
  typedef imgbase<512, 512, float> bildfloat_base;


  class bildfloat : public bildfloat_base {
    public:
      bildfloat();
      bildfloat(const bild16& aufnahme, const eichkurve& ek);
      bildfloat(const bild16& aufnahme);
      bildfloat(const region_t& regio);
      bildfloat(const char* name);
      ~bildfloat();
      static const float pix_size;
      void write_gnuplot (const char *name) const;
      bildfloat& read_gnuplot (const char *name) THROWSPEC;
      
      bildfloat& from_bild16(const bild16& aufnahme, const eichkurve& ek);
      bildfloat& from_bild16(const bild16& aufnahme);
      bildfloat& apply_eichkurve(const eichkurve& ek);
      bildfloat& convolve (const float * kernel, int sx, int sy, bool normalize=true) const THROWSPEC;
      bildfloat& gaussian(float radius); //Glätten durch Faltung mit einer Gaußglocke
      void shutter_detect(int  &y0, int &y1, float thresh, int radius) const;
      bildfloat& shutter_correct(int y0, int y1, float shadefactor);
      bildfloat& gamma_correct(const bildfloat& gammabild);
      bildfloat& bquad_correct(const bildfloat& bquadbild);
      bildfloat& polyfit_correct(const QRfit& korrektor);
      bildfloat& badpix_correct (const bild16& bpmap, bild16::data_t threshold=bild16::intensity_max/2);
      bildfloat& mask(const region_t& reg); // set all points not in reg to NaN
  };    
  
  inline bildfloat::bildfloat() : bildfloat_base() { }
  
  inline bildfloat::bildfloat(const char *name) : bildfloat_base() 
  { 
    read_gnuplot(name);
  }
  
  inline bildfloat::bildfloat(const bild16& aufnahme, const eichkurve& ek) : bildfloat_base() 
  { 
    from_bild16(aufnahme, ek); 
  }
  
  inline bildfloat::bildfloat(const bild16& aufnahme) : bildfloat_base() 
  { 
    from_bild16(aufnahme); 
  }
  

  inline bildfloat::~bildfloat()  { }

  typedef std::valarray<float> lut;
  typedef imgbase<512, 512, lut> interpol_base;
  

  class interpol : public interpol_base {
   public: 
    interpol(bild16 * refbilder, float* refw, int nrefimg);
    interpol(const char *fname);
    void write(const char *fname);
    void read (const char *fname) throw (const char*);
    float operator() (int y, int x, float val);
    const lut& operator() (int y, int x) const { return interpol_base::operator()(y,x); }
    
   private:
    static const char magic[];
    int nref;
    lut refwerte;
    lut& operator() (int y, int x) { return interpol_base::operator() (y,x); }
  };
  

class QRfit : public interpol_base {
 public:
  QRfit(bild16 *refbilder, float *refw, int nrefimg, int order_=4);
  QRfit(const char *fname);
  ~QRfit();
  void write(const char *fname) const;
  void read (const char *fname) THROWSPEC;

  float operator() (int y, int x, float val) const;
  const lut& operator() (int y, int x) const { return interpol_base::operator()(y,x); }
  float get_coeff(int x, int y, int n) const;
 private:
  int nref, order;
  lut refwerte;
  lut& operator() (int y, int x) { return interpol_base::operator() (y,x); }

  void polyfit(double *bildwerte);
  double **R;
  double *d, *b, *koeff;
  static const char magic[];
};

 
#endif

---