imgclass.cpp

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <cstdlib>
#include <cstring>
using namespace std;

#include "imgclass.hpp"
#include "matrixutil.hpp"


namespace NaN {
  const float fnan = nanf("");
}

using NaN::fnan;


inline double clip(double x, double limit) {
  return (x>limit)?limit:x;
}

//-----------------------------------------------
// alles mit bild16
// 16-Bit Graustufenbilder
// aus der Röntgenanlage
// Ausgabe als pgm-Bitmap möglich
//---------------------------------------------------

const bild16::data_t bild16::intensity_max=0xFFFF;
const bild16::data_t bild16::intensity_min=0x0000;
const streamoff bild16::HIS_START=100;  

bild16::bild16 (const char *fname) THROWSPEC : bild16_base () 
{
  file_header = new header_t;
  if (fname) {
    read_his (fname);
  } else {
    //initialise header
    file_header->FileID = 0x7000;       //1
    file_header->HeaderSize = 68;
    file_header->HeaderVersion = 100;
    file_header->FileSize = 524388;
    file_header->ImageHeaderSize = 32;  //5
    file_header->ULX = 1;
    file_header->ULY = 1;
    file_header->BRX = 512;
    file_header->BRY = 512;
    file_header->NrOfFrames = 1;        //10
    file_header->Correction = 0;
    file_header->IntegrationTime = 0;
    file_header->TypeOfNumbers = 4;     //13
  }


}

bild16::bild16 (unsigned short *data) THROWSPEC : bild16_base ()
{
    file_header = new header_t;
    //initialise header
    file_header->FileID = 0x7000;       //1
    file_header->HeaderSize = 68;
    file_header->HeaderVersion = 100;
    file_header->FileSize = 524388;
    file_header->ImageHeaderSize = 32;  //5
    file_header->ULX = 1;
    file_header->ULY = 1;
    file_header->BRX = 512;
    file_header->BRY = 512;
    file_header->NrOfFrames = 1;        //10
    file_header->Correction = 0;
    file_header->IntegrationTime = 0;
    file_header->TypeOfNumbers = 4;     //13
    
    MAP ((*this)) {
      bild[y][x]=data[SIZEX*y+x];
    }  
         
}

bild16::~bild16() 
{
  delete file_header;  
}

void bild16::write_pgm16 (const char *name) const 
{
  ofstream datei (name, ios::binary);
  datei << "P5" << '\n';        //magic number für PGM
  datei << SIZEX << " " << SIZEY << "\n" << intensity_max << "\n";      //Breite, Höhe und größter Wert für Weiß, 16 Bit 

  for (int y = 0; y < SIZEY; y++)
    for (int x = 0; x < SIZEX; x++) {
      datei << static_cast <unsigned char>(bild[y][x] / 0x100) 
      << static_cast <unsigned char>(bild[y][x] % 0x100);
    }

  // Daten als Block, byte-weise
  // MSB first
 }


void bild16::write_pgm8 (const char *name) const 
{
  ofstream datei (name, ios::binary);
  datei << "P5" << endl;        //magic number für PGM
  datei << SIZEX << " " << SIZEY << "\n" << 0xFF << "\n";       //Breite, Höhe und größter Wert für Weiß, 8 Bit 

  for (int y = 0; y < SIZEY; y++)
    for (int x = 0; x < SIZEX; x++) {
      datei << static_cast < unsigned char >(bild[y][x] / 0x100);
    }

  // Daten als Block, byte-weise
  // Koordinatensystem nicht mathematisch, sondern
  // bildschirmorientiert: y geht nach unten
}

void bild16::write_pgm8_bounds (const char *name, bild16::data_t low, bild16::data_t high) const 
{  if (low==high) STHROW("Error: Zero range (write_pgm8_bounds), low==high=="<<low);
   if (low>high) swap(low,high);
  ofstream datei (name, ios::binary);
  datei << "P5" << endl;        //magic number für PGM
  datei << SIZEX << " " << SIZEY << "\n" << 0xFF << "\n";       //Breite, Höhe und größter Wert für Weiß, 8 Bit 

  for (int y = 0; y < SIZEY; y++)
    for (int x = 0; x < SIZEX; x++) {
      unsigned char byte;
      data_t wert=(*this)(y,x);
      if (wert>high) byte=0xFF;
      else if (wert<low) byte=0x00;
      else byte = static_cast < unsigned char >(0xFF*double((*this)(y,x)-low)/double(high-low));
      datei << byte;
    }

  // Daten als Block, byte-weise
  // Koordinatensystem nicht mathematisch, sondern
  // bildschirmorientiert: y geht nach unten
}

bild16& bild16::read_his (const char *name) THROWSPEC
{
  ifstream his;                 //Unsere Datei-Streams
  his.open (name, ios::binary | ios::in);       // öffnen des hisfiles
  if (!his) STHROW(name<<" not found!\n");
  
  // Einlesen des Fileheaders
  his.read ((char *) file_header, sizeof (*file_header));


  // Einlesen der Daten
  his.seekg (HIS_START);
  for (int y = 0; y < SIZEY; y++) {
    his.read ((char *) bild[y], sizeof (unsigned short) * SIZEX);       
  }

  his.close ();                 //Schliesst das hisfile
  return *this;
}

void bild16::write_his (const char *name) const
{
  ofstream his;         
  his.open (name, ios::binary); // öffnen des hisfiles
  if (!his) {
    STHROW("Cannot write to "<<name); 
  }
  // Schreiben des Fileheaders
  his.write ((char *) file_header, sizeof (*file_header));
  
  // padding 
  for (streamoff i=sizeof(*file_header); i<HIS_START; i++)
    his.put(char(0));

  // Schreiben der Daten
  his.seekp (HIS_START);
  // nur zur Sicherheit ein seek, die Position 
  // sollte eigentlich schon stimmen!
  for (int y = 0; y < SIZEY; y++) {
    his.write ((char *) bild[y], sizeof (unsigned short) * SIZEX);      
  }

  his.close ();                 //Schliesst das hisfile
}

inline int isdecdigit(char dig) {
  return (dig>='0' && dig<='9');
}

inline void eatwhite(istream& is) {
 bool yet=true;
 do {
   switch (is.peek()) {
     case '\n': 
     case '\r':
     case ' ' :
     case '\t':
       is.get();
       break;
     default:
       yet=false;
   }
 } while(yet && !is.eof());
}

static int getpgmnumber(istream &is) THROWSPEC{
  // Eat whitespace, ignore comments, then read decimal number
  
  eatwhite(is);
  while (is.peek()=='#') {
    // Comment, ignore to eol or eof
    int c=is.get();
    while ((c!='\n') && (!is.eof())) c=is.get();
    if (is.eof()) STHROW ("Expected number (PGM)");
    eatwhite(is);
  }
  // next char is no #, should read number

  const int maxlen=10; // maximum tendigits
  char digitbuf[maxlen+1];
  int dnum=0;
  do {
    
    char c=is.get();
    if (is.eof()) STHROW("Premature EOF (PGM)");

    if (isdecdigit(c)) digitbuf[dnum++]=c;
    else {
      is.putback(c);
      break;
    }  
  } while (dnum < maxlen);

  digitbuf[dnum]='\0';

  return atoi(digitbuf);
}
    
bild16& bild16::read_pgm (const char *name) THROWSPEC 
{
  ifstream pgm;                 //Unsere Datei-Streams
  pgm.open (name, ios::binary | ios::in);       // öffnen des hisfiles
  if (!pgm) STHROW(name<<" not found!");
  try {
    // PGM format: P5<whitespace>
    // Dann Breite und Höhe, evtl Kommentare dazwischen
    if (pgm.get()!='P') STHROW("Not a PPM file");
    if (pgm.get()!='5') STHROW("Not a PGM Graymap");
  
    int width=getpgmnumber(pgm);
    if (width!=SIZEX) STHROW("Width doesn't match. Expected"<<SIZEX<<", got"<<width);
    
    int height=getpgmnumber(pgm);
    if (height!=SIZEY) STHROW("Height doesn't match. Expected"<<SIZEY<<", got"<<height);
    
    int maxval=getpgmnumber(pgm);
    pgm.get(); // ein einzelnes Whitespace ist nach Maxval erlaubt
    switch (maxval)
    {
       case 255:
        // 8bit PGM
         for (int y=0; y < SIZEY; y++) {
           for (int x=0; x< SIZEX; x++) {
           unsigned int val8=pgm.get();
           (*this)(y,x)=val8*0x100;
         }
         }
         if (!pgm) STHROW("PGM corrupt");
         break;
       case 65535:
         // 16bit PGM
         for (int y=0; y < SIZEY; y++) {
           for (int x=0; x< SIZEX; x++) {
           unsigned int hi8=pgm.get();
           unsigned int lo8=pgm.get();
           (*this)(y,x)=(hi8 & 0xFF)*0x100+(lo8 & 0xFF);
         }
         }
         if (!pgm) STHROW("PGM corrupt");
         break;
       default:
         STHROW("Can only read 8bit or 16bit PGM (maxval=="<<maxval<<")");
    }     
  
    pgm.close ();                       //Schliesst das hisfile
  } catch (const string& msg) {
      STHROW(name<<": "<<msg);
      // prefix every error by the name
  }    
     
  return *this;
}

void bild16::fromfloat(const bildfloat& bf, float max) {
  MAP (bf) {
      if (bf(y,x)>max) (*this)(y,x)=intensity_max;
      else if (bf(y,x)<0) (*this)(y,x)=intensity_min;
      else (*this)(y,x)=intensity_min+static_cast<data_t>(intensity_max*bf(y,x)/max);
    }  
}

      
bild16& bild16::og_correct(const bild16& offset, const bild16& gain) {
  const double epsilon=1.0; 
  for (int y=0; y<SIZEY; y++) 
    for (int x=0; x<SIZEX; x++) {
       // Offset-Gain-Korrektur
       double denom = gain(y,x)-offset(y,x);
       if (denom <= epsilon) continue;
         //vermeide Division durch Null
       double nomin = bild[y][x] - offset(y,x);
       if (nomin < 0) continue;
         // dunkler als Offset = kaputter Pixel
       bild[y][x] = data_t(clip(double(intensity_max)*nomin/denom, intensity_max));
    }
    return *this;
}    


bild16& bild16::drawx(int x, int y1, int y2) THROWSPEC {
        
  if (y1>y2) swap(y1,y2);
  
  if ((x<0) || (x>=SIZEX)) STHROW("drawx: x out of range (x="<<x<<")");
  if ((y2<0) || (y1>=SIZEY)) STHROW("drawx: line outside image (y1="<<y1<<", y2="<<y2<<")");
  
  /* clip line on border */
  if (y1<0) y1=0;
  if (y2>=SIZEY) y2=SIZEY-1;
  
  for (int y=y1; y<=y2; y++)
    bild[y][x]^=0xFFFF;
  return *this;
}

bild16& bild16::drawy(int y, int x1, int x2) THROWSPEC {
  if (x1>x2) swap(x1,x2);
  
  if ((y<0) || (y>=SIZEY)) STHROW("drawy: y out of range (y="<<y<<")");
  if ((x2<0) || (x1>=SIZEX)) STHROW("drawy: line outside image (x1="<<x1<<", x2="<<x2<<")");
  
  /* clip line on border */

  if (x1<0) x1=0;
  if (x2>=SIZEX) x2=SIZEX-1;
  
  for (int x=x1; x<=x2; x++)
    bild[y][x]^=0xFFFF;
  return *this;  
}

bild16& bild16::kreuz (int x, int y) THROWSPEC {

  if ((x<0) || (x>=SIZEX) || (y<0) || (y>=SIZEY)) STHROW("kreuz: Outside of the image, x="<<x<<", y="<<y<<")");
  
  drawx(x, y-10, y+10);
  drawy(y, x-10, x+10);
  return *this;
} 
//---------------------------------------------------
// Eichkurven
//---------------------------------------------------

eichkurve::eichkurve ()
{
  // reserve memory for the lookup-table
  lookup= new double[0x10000];
}

eichkurve::~eichkurve() {
  // release memory
  delete[] lookup;
}
  
const char eichkurve::magic[]="EICHTABELLE";

void eichkurve::table_init(double hmin, double hmax) {
 //solve the function for all possible input values;
  for (int i=0; i<0x10000; i++) {
    lookup[i]=fsolve(i, hmin, hmax);
  }  

}

double eichkurve::fsolve(bild16::data_t val, double hmin, double hmax) const
{
  const double epsilon=0.1;            // nicht genauer als 1 Bit Höhen-Auflösung
  double vmin=eichfun(hmin);
  double vmax=eichfun(hmax);
  int signum=1;
  if (vmin<vmax) { 
    // Funktion (monoton) steigend
    signum=1; 
    if (val<vmin) return hmin;
    if (val>vmax) return hmax;
  } else {
     //Funktion (monoton) fallend
     signum=-1;
     if (val>vmin) return hmin;
     if (val<vmax) return hmax;
  } 
  //Wert ist zwischen vmin und vmax eingeschlossen
  //Sekantenverfahren (Regula falsi)
  double h;
  double vh;
  do {
    h  = (val-vmin)/(vmax-vmin) * (hmax-hmin) + hmin;
    vh = eichfun(h);
   // cout<<hmax<<"  "<<hmin<<"  "<<h<<"  "<<vh<<"  "<<val<<endl;
    if (signum*(vh-val)>0) {
        vmax = vh;
        hmax = h;
    } else {
        vmin = vh;
        hmin = h;
    }   
    
  } while (abs(vh-val)>epsilon);
  
  return h; 
}  


// Newton-Raphson; konvergiert nicht immer
#if 0 
double eichkurve::fsolve(bild16::data_t val, double hmin, double hmax) const
{
  const double epsilon=0.5;            // nicht genauer als 1 Bit Höhen-Auflösung
  double vmin=eichfun(hmin);
  double vmax=eichfun(hmax);
  if (vmin<vmax) { 
    // Funktion (monoton) steigend
    if (val<vmin) return hmin;
    if (val>vmax) return hmax;
  } else {
     //Funktion (monoton) fallend
     if (val>vmin) return hmin;
     if (val<vmax) return hmax;
  } 
  double x = initguess(val); //Startwert schätzen
  double xalt;
  do {
    xalt = x;
    x = xalt + (val - eichfun(xalt))/deichfun_dx(xalt);
  } while (abs(eichfun(xalt)-eichfun(x))>epsilon);
  //Newton-Raphson Iteration
  
  // Beschränkung, falls Newton außerhalb des erlaubten Bereichs landet
  if (x>hmax) { 
           cout<<"Zu groß: "<<x<<"  "<<val<<endl;
           return hmax;
  }        
  if (x<hmin)  return hmin;
  return x; 
}  
#endif 


void eichkurve::table_write(const char *fname) const {
  ofstream datei(fname, ios::out|ios::binary);
  datei.write(magic, sizeof(magic));
  for (int i=0; i<0x10000; i++)
    datei<<rawfloat(lookup[i]);
}       
//
// Exponentialabschwächung, versch. char. Zeiten
expkurve::expkurve (double  sA, double sg1, double sg2, double sg3, 
                      double sf1, double sf2) :
  eichkurve(), A(sA), g1(sg1), g2(sg2), g3(sg3), f1(sf1), f2(sf2) { 

  table_init(0.0, 40.0);
} 
 
double expkurve::eichfun(double x) const {
    return  A*(f1*exp(-g1*x)+f2*exp(-g2*x)+(1-f1-f2)*exp(-g3*x));
}
  
double expkurve::deichfun_dx(double x) const {
    return -A*(f1*g1*exp(-g1*x)+f2*g2*exp(-g2*x)+(1-f1-f2)*g3*exp(-g3*x));
  
}

double expkurve::initguess(bild16::data_t val) const {
  return - log(double(val)/A)/g1;
}  

// 4 überlagerte Exponentialkurven
exp4kurve::exp4kurve (double sA, double sa1, double sa2, double sa3,
                // rel. Anteile
                double sb1, double sb2, double sb3, double sb4 ) : // charakteristische Tiefen
  eichkurve(), A(sA), a1(sa1), a2(sa2), a3(sa3), b1(sb1), b2(sb2), b3(sb3), b4(sb4) { 

  table_init(0.0, 35.0);
} 
 
double exp4kurve::eichfun(double x) const {
    return A*(a1*exp(-b1*x)+a2*exp(-b2*x)+a3*exp(-b3*x)+(1.0-a1-a2-a3)*exp(-b4*x));
}
  
double exp4kurve::deichfun_dx(double x) const {
    return -A*(a1*b1*exp(-b1*x)+a2*b2*exp(-b2*x)+a3*b3*exp(-b3*x)+(1.0-a1-a2-a3)*b4*exp(-b4*x));
  
}

double exp4kurve::initguess(bild16::data_t val) const {
  return -log(val/(A*a3))/b3;
}  


//
// flatkurve

flatkurve::flatkurve (bool normalize_) : eichkurve(), normalize(normalize_) {

  for (int i=bild16::intensity_min; i<=bild16::intensity_max; i++)
    lookup[i]=initguess(i);
}

double flatkurve::initguess(bild16::data_t  val) const {
  if (normalize) return static_cast<double>(val)/bild16::intensity_max;
  else return static_cast<double>(val);
}

double flatkurve::eichfun(double x) const {
  if (normalize) return x*bild16::intensity_max;
  else return x;
}

double flatkurve::deichfun_dx(double x) const  {
  if (normalize) return static_cast<double>(bild16::intensity_max);
  else return 1.0;
}

// tablekurve schlägt einfach in der Tabelle nach und
// benutzt Dateien, die vorher von anderen Kurven geschrieben wurden
// Das spart die aufwendigen Newton-Iterationen am Anfang
tablekurve::tablekurve (const char *fname) THROWSPEC {
  ifstream datei(fname, ios::binary);
  if (!datei) STHROW(fname<<" not found.");
  char magictest[sizeof(magic)];
  
  for (size_t i=0; i<sizeof(magic); i++) 
    datei.get(magictest[i]);
  
  if (strncmp(magic, magictest, sizeof(magic))!=0)
    STHROW(fname<<" is not a tablekurve file.");
  
  for (int i=0; i< 0x10000; i++) {
    float temp;
    datei>>irawfloat(temp);
    lookup[i]=temp;
  }  
}    


double tablekurve::eichfun (double x) const 
{
  STHROW ("Tablekurve pseudo-function called");
}
double tablekurve::deichfun_dx (double x) const
{
  STHROW ("Tablekurve pseudo-function called");
}

double tablekurve::initguess (bild16::data_t val) const 
{
  STHROW ("Tablekurve pseudo-function called");
}


//---------------------------------------------------
// bildfloat
// 512x512-Bilder in single precision floats
// für alle Rechnungen gut
// Ein/Ausgabe als gnuplot binary file möglich
//---------------------------------------------------




const float bildfloat::pix_size  = static_cast<float>(1606.0 / 1765 * 0.4);

bildfloat::bildfloat(const region_t& regio) : bildfloat_base() {
  MAP (regio)
   if (regio(y,x)) (*this)(y,x)=1.0f;
   else (*this)(y,x) = fnan; 
}      
  
void bildfloat::write_gnuplot (const char *name) const
{
  // write binary gnuplot data file
  ofstream datfile (name, ios::binary);
  datfile << rawfloat (SIZEY);
  for (int y = 0; y < SIZEY; y++)
    datfile << rawfloat (y * pix_size);
    
  for (int y = 0; y < SIZEY; y++) {
    datfile << rawfloat (y * pix_size);
    
    for (int x = 0; x < SIZEX; x++) {
      datfile << rawfloat ((*this)(y, x));
    }  
  }
  datfile.close();
}

bildfloat& bildfloat::read_gnuplot (const char *name) THROWSPEC
{
  // read binary gnuplot data file
  ifstream datfile (name, ios::binary);
  if (!datfile) STHROW (name<<" not readable");
  float dummy;
  datfile >> irawfloat (dummy);
  if (dummy != SIZEY) STHROW ("Not a known gnuplot binary format");
  for (int y = 0; y < SIZEY; y++)
    datfile >> irawfloat (dummy); //ignore these values (y0..yn)
    
  for (int y = 0; y < SIZEY; y++) {
    datfile >> irawfloat (dummy); //ignore x value
    
    for (int x = 0; x < SIZEX; x++) {
      datfile >> irawfloat ((*this)(y, x));
    }  
  }
  datfile.close();
  return *this;
}
     
bildfloat& bildfloat::from_bild16(const bild16& aufnahme) {

  MAP (aufnahme) 
      bild[y][x]=static_cast<float>(aufnahme(y,x))/static_cast<float>(bild16::intensity_max);
  return *this;    
}

bildfloat& bildfloat::from_bild16(const bild16& aufnahme, const eichkurve& ek) {

  MAP (aufnahme) 
      bild[y][x]=ek(aufnahme(y,x));
  return *this;    
}


bildfloat& bildfloat::badpix_correct (const bild16 & bpmap, bild16::data_t threshold)   //badpixel correction
{
#define GOODPIX(Y,X)    if (bpmap((Y),(X)) < threshold) \
        { \
          mean += (*this)((Y),(X)); \
          mean_i++;  \
        }

  for (int y = 1; y < SIZEY - 1; y++)
    for (int x = 1; x < SIZEX - 1; x++) {
      float mean = 0;
      int mean_i = 0;
      if (bpmap(y,x) > threshold)
        //ersetze durch mittelwert der 8 nachbarn, wenn diese ok
      {
          GOODPIX (y - 1, x - 1)
          GOODPIX (y - 1, x)
          GOODPIX (y - 1, x + 1)
          GOODPIX (y, x - 1)
          GOODPIX (y, x + 1)
          GOODPIX (y + 1, x - 1)
          GOODPIX (y + 1, x)
          GOODPIX (y + 1, x + 1)
          if (mean_i > 0) {
            (*this) (y, x) = mean / mean_i;     // mittelwert der nicht defekten Nachbarpix
          } else                        //wenn alle 8 nachbarpix defekt sind?
          //es muss div durch null vermieden werden
          {
            (*this) (y, x) = 0;
          }


      }                         //if-zu
    }                           // for-zu

    return *this;
}

bildfloat& bildfloat::gamma_correct (const bildfloat& gammabild) {
 for (int y = 1; y < SIZEY - 1; y++)
    for (int x = 1; x < SIZEX - 1; x++) {
       if ((*this)(y,x)<=0.0) { 
         (*this)(y,x)=0.0; 
         continue;
       }

       if ((*this)(y,x)>=1.0) {
         (*this)(y,x)=1.0;
         continue;
       }         
       (*this)(y,x)=pow((*this)(y,x), gammabild(y,x));
 }
 return *this;
} 

bildfloat& bildfloat::bquad_correct (const bildfloat& bquadbild) {
 for (int y = 1; y < SIZEY - 1; y++)
    for (int x = 1; x < SIZEX - 1; x++) {
       if ((*this)(y,x)<=0.0) { 
         (*this)(y,x)=0.0; 
         continue;
       }

       if ((*this)(y,x)>=1.0) {
         (*this)(y,x)=1.0;
         continue;
       }         
       (*this)(y,x)=((*this)(y,x)+bquadbild(y,x)*SQR((*this)(y,x)))/(1.0+bquadbild(y,x));
    }
    return *this;
} 


bildfloat& bildfloat::polyfit_correct (const QRfit& korrektor) {
    /* interpoliere durch das Ausgleichspolynom */
     for (int y=0; y<SIZEY; y++) {
       for (int x=0; x<SIZEX; x++) {
        (*this)(y,x)=korrektor(y, x, ((*this)(y,x)) );
       }
     }
   return *this;
}

bildfloat& bildfloat::apply_eichkurve(const eichkurve& ek) {
   MAP ((*this)) {
       float val = (*this)(y,x);
       // ensure we are not out of range
       if (val>1.0) (*this)(y,x)=ek(bild16::intensity_max);
       else if (val<0.0) (*this)(y,x)=ek(bild16::intensity_min);
       else (*this)(y,x)=ek(static_cast<bild16::data_t>(bild16::intensity_max*val));
   }     
  return *this;    
}

void bildfloat::shutter_detect(int &y0, int &y1, float thresh, int radius ) const {
  // looks for the dark shadow 
  const int xmitte=SIZEX/2; const int ymitte=SIZEY/2;
  y0=-1; y1=-1;
  for (int y=SIZEY/2; y < SIZEY-1; y ++) {
     float sum=0; float diffsum=0;
     for (int x=0; x<SIZEX; x++) 
       if ((SQR(x-xmitte)+SQR(y-ymitte))<SQR(radius)) {
       sum+=(*this)(y,x)+(*this)(y+1,x);
       diffsum+=(*this)(y,x)-(*this)(y+1,x);
     }
     if (sum<=0) continue;
     if (fabs(diffsum/sum)>thresh) {
       y1 = y+1;
       break;
     }  
   }

   for (int y=SIZEY/2; y; y--) {
     float sum=0; float diffsum=0;
     for (int x=0; x< SIZEX; x++) 
       if ((SQR(x-xmitte)+SQR(y-ymitte))<SQR(radius)){
       sum+=(*this)(y,x)+(*this)(y+1,x);
       diffsum+=(*this)(y-1,x)-(*this)(y,x);
     }
     
     if (sum<=0) continue;
     if (fabs(diffsum/sum)>thresh) {
       y0 = y-1;
       break;
     } 
   }
 
}

bildfloat& bildfloat::shutter_correct(int y0, int y1, float shadingfactor) {
  for (int y=0; y<=y0; y++) 
      for (int x=0; x< (*this).SIZEX; x++)
        (*this)(y,x)*=shadingfactor;

    for (int y=y1; y<(*this).SIZEY; y++) 
      for (int x=0; x< (*this).SIZEX; x++)
        (*this)(y,x)*=shadingfactor;

   return *this;        
}

bildfloat& bildfloat::convolve (const float * kernel, int sx, int sy, bool normalize) const  THROWSPEC  {
  if ( (sx %2 ==0) || (sy %2 ==0) || (sx<=0) || (sy<=0)) 
      STHROW("Sizes of kernel must be odd (convolve) and positive, sx="<<sx<<", sy="<<sy);

  float *newkernel=NULL;
  if (normalize) {
    // normalize kernel
    float norm=0;
    for (int k=0; k<sx*sy; k++) norm+=kernel[k];
    if (norm==0.0) STHROW("Kernel may not be zero (convolve)");
    newkernel = new float[sx*sy];
    for (int k=0; k<sx*sy; k++) newkernel[k]=kernel[k]/norm;
    kernel=newkernel; // Pointer umbiegen erlaubt
  }
  
  // construct convoluted image
  bildfloat* newbild = new bildfloat;
  
  for (int y=0; y<SIZEY; y++)
    for (int x=0; x<SIZEX; x++) {
      if ((x< sx/2) || (y < sy/2) || (x>= SIZEX - sx/2) || (y>=SIZEY - sy/2)) {
        (*newbild)(y,x)=(*this)(y,x);
        continue;
      } 
      // Randpunkte ausschließen
      float value=0.0;
      for (int ky=0; ky<sy; ky++)
        for (int kx=0; kx<sx; kx++)
          value+=(*this)(y+ky-sy/2,x+kx-sx/2)*kernel[ky*sx+kx];
      (*newbild)(y,x)=value;
    }
  if (normalize) delete[] newkernel;
  //Smart pointers could help, but confuse EP5 people...
  return *newbild;
}  
       
bildfloat& bildfloat::gaussian (float radius) {
   const int size=2*static_cast<int>(radius)+1;
   const int middle=size/2;
   float *kernel= new float[SQR(size)];
   for (int y=0; y<size; y++)   
    for (int x=0; x<size; x++) {
      kernel[y*size+x] = exp(-(SQR(x-middle)+SQR(y-middle))/radius);
     // cerr<<kernel[y*size+x]<<"  ";
    }  
    
    //ugly and SLOW
    //need to copy temporary unnecessarily
    bildfloat *p=&convolve(kernel, size, size, true);
   (*this)=*p;
   delete p;
   delete [] kernel;
   return *this;
 
 } 

bildfloat& bildfloat::mask(const region_t& reg) {
  for (int y=0; y<static_cast<signed>(SIZEY); y++)
    for (int x=0; x<static_cast<signed>(SIZEX); x++)
      if (!reg(y,x)) (*this)(y,x)=fnan;
  return *this;
}

const char interpol::magic[]="INTERREF";

interpol::interpol(const char *fname) : interpol_base() {
  read(fname);
}  

interpol::interpol(bild16 * refbilder, float* refw, int nrefimg): 
  interpol_base(), nref(nrefimg), refwerte(refw, nrefimg) {

  /* Kopiere die Referenzbilder in eine "transponierte" Struktur */
  /* eine Matrix 512x512 von valarrays, die die Kurve enthält */
  for (int y=0; y<SIZEY; y++) {
    for (int x=0; x<SIZEX; x++) {
      (*this)(y,x).resize(nref);

      for (int r=0; r<nref; r++)
          (*this)(y,x)[r] = static_cast<float>(refbilder[r] (y,x))/bild16::intensity_max;
     }    
   }

}  

void interpol::write (const char* fname) {
  ofstream datei(fname, ios::out|ios::binary);
  datei.write(magic, sizeof(magic));
  datei.put(static_cast<char>(nref));
  for (int r=0; r<nref; r++) {
    datei<<rawfloat(refwerte[r]);
  }

  for (int y=0; y<SIZEY; y++)
    for (int x=0; x<SIZEX; x++)
      for (int r=0; r<nref; r++)
        datei<<rawfloat((*this)(y,x)[r]);
}       
    
  

void interpol::read(const char* fname) throw(const char*){
  ifstream datei(fname, ios::binary);
  if (!datei) throw("Nicht gefunden.");
  char magictest[sizeof(magic)];
  
  for (size_t i=0; i<sizeof(magic); i++) 
    datei.get(magictest[i]);
  
  if (strncmp(magic, magictest, sizeof(magic))!=0)
    throw("Kein Interpol-Reffile");
  
  char c;
  datei.get(c);
  nref=c;
  
  refwerte.resize(nref);
  for (int r=0; r<nref; r++)
    datei>>irawfloat(refwerte[r]);
  
  for (int y=0; y<SIZEY; y++)
    for (int x=0; x<SIZEX; x++) {
      (*this)(y,x).resize(nref);
      for (int r=0; r<nref; r++) 
        datei>>irawfloat((*this)(y,x)[r]);
    }
    
}    
  
float interpol::operator () (int y, int x, float val) {
   if (val>=(*this)(y,x)[nref-1]) {
     return 1.0;
   }
   if (val<=(*this)(y,x)[0]) {
     return 0.0;
   }
   int refindex=0;
   while (val > (*this)(y,x)[refindex+1]) refindex++;
   
   /* zu interpolierender Wert liegt zwischen refindex und refindex+1 */

   /* vermeide Division durch Null */
   if ((*this)(y,x)[refindex+1]==(*this)(y,x)[refindex]) {
      return val;
   }     

   /* interpoliere linear */
   return  refwerte[refindex] + (refwerte[refindex+1] - refwerte[refindex]) * 
                (val-(*this)(y,x)[refindex])/((*this)(y,x)[refindex+1]-(*this)(y,x)[refindex]);
}               


const char QRfit::magic[]="QRfit";

QRfit::QRfit(bild16 *refbilder, float *refw, int nrefimg, int order_) :
   interpol_base(), nref(nrefimg), order(order_), refwerte(refw, nrefimg) {
  
  reserve_2D(R, order, nref);
  d=new double[order];
  b=new double[nref];
  koeff=new double[order];
  double *bw = new double[nref];
  for (int y=0; y<SIZEY; y++) {
    for (int x=0; x<SIZEX; x++) {
      (*this)(y,x).resize(order);

      for (int k=0; k<nref; k++)
        bw[k]=double((refbilder[k])(y,x))/double(bild16::intensity_max);
        
      polyfit(bw);

      for (int k=0; k<order; k++)
        (*this)(y,x)[k]=koeff[k];
      
      if ((x==100) && (y==100)) {
        for (int r=0; r<order; r++)
           cout<<koeff[r]<<"  ";
         cout<<endl; 
      }
      // Fitte ein Polynom an die gewünschte Übertragungskurve
      // und speichere in dem valarray ab
    }  
  }
  delete[] bw;
} 

QRfit::QRfit(const char *fname) : interpol_base() {
/*  reserve_2D(R, order, nref);
  d=new double[order];
  b=new double[nref];
  koeff=new double[order];
  */
  R=NULL;
  d=NULL;
  b=NULL;
  koeff=NULL;
  read(fname);
}

QRfit::~QRfit() {
  if (R) delete_2D(R);
  if (d) delete[] d;
  if (b) delete[] b;
  if (koeff) delete[] koeff;
}

void QRfit::write (const char* fname) const {
  ofstream datei(fname, ios::out|ios::binary);
  datei.write(magic, sizeof(magic));

  datei.put(static_cast<char>(nref));
  datei.put(static_cast<char>(order));
  
  for (int r=0; r<nref; r++) {
    datei<<rawfloat(refwerte[r]);
  }

  for (int y=0; y<SIZEY; y++)
    for (int x=0; x<SIZEX; x++)
      for (int k=0; k<order; k++)
        datei<<rawfloat((*this)(y,x)[k]);
}       
    
  

void QRfit::read(const char* fname) THROWSPEC {
  ifstream datei(fname, ios::binary);
  if (!datei) STHROW(fname<<" not found.");
  char magictest[sizeof(magic)];
 
  /* Testen, ob die Datei mit dem Magicstring beginnt */
  
  for (size_t i=0; i<sizeof(magic); i++) 
    datei.get(magictest[i]);
  
  if (strncmp(magic, magictest, sizeof(magic))!=0)
    STHROW(fname<<" is not a QRfit file");
  
  /* Anzahl der Referenzwerte und Ordnung des Polynoms lesen */
  char c;
  datei.get(c);
  nref=c;
  
  datei.get(c);
  order=c;
  
  /* Referenzwerte lesen */
  refwerte.resize(nref);
  for (int r=0; r<nref; r++)
    datei>>irawfloat(refwerte[r]);
  
  /* Koeffizienten lesen */
  for (int y=0; y<SIZEY; y++)
    for (int x=0; x<SIZEX; x++) {
      (*this)(y,x).resize(order);
      for (int k=0; k<order; k++) 
        datei>>irawfloat((*this)(y,x)[k]);
    }
    
}    
 
float QRfit::operator () (int y, int x, float val) const {
  float sum=(*this)(y,x)[0];
  float var=1.0;
  /* werte Polynom aus */
  for (int k=1; k<order; k++) {
    var*=val;
    sum+=var*(*this)(y,x)[k];
  }
   if (sum<0.0) return 0.0;
   if (sum>1.0) return 1.0;
  return sum;
}

float QRfit::get_coeff(int x, int y, int n) const {
  if (n<0 || n>order) STHROW("Index must be positive and < "<<order);
  return (*this)(y,x)[n];
}


void QRfit::polyfit(double *bildwerte) {
  
  // Fülle die Matrix mit den Bildwerten^ordnung
  // Für Polynomfit
  for (int nr=0; nr<nref; nr++) {
    double r=bildwerte[nr];
    R[nr][0]=1;
    
    for (int o=1; o<order; o++) {
       R[nr][o]=R[nr][o-1]*r;
    }

  }

  // Fülle die rechte Seite mit den Referenzwerten
  for (int nr=0; nr<nref; nr++)
    b[nr]=refwerte[nr];
  
  // Nun löse mit dem QR-Algorithmus
  // || R x - b || = minimal nach x

  if (QR(nref, order, R, d, b, 1e-10))  {
    //  STHROW("Singuläre Matrix");
    // Nicht gleich Exception werfen
    // "Notlösung" Offset-Gain für diesen Pixel
    cerr<<".";
    double offset=refwerte[0]; 
    double gain=refwerte[nref-1];
    if (offset==gain) {
      // jetzt ist wirklich alles singulär
      cerr<<":";
      koeff[0]=0.0;
      koeff[1]=1.0;
    } else {  
      koeff[0]=- offset/(gain-offset);
      koeff[1]=1.0/(gain-offset);
      // Standard Offset-Gain
    }
    for (int o=2; o<order; o++) koeff[o]=0.0;

  }  else  { 
   
    // Fit war nichtsingulär -- löse das QR-Ergebnis auf
    rdsolve(order, R, d, b, koeff);

  }  
}    

---