UnVignette.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "UnVignette.hh"
#include <iostream>
#include <math.h>
#include <MathAlgo/Lapack.hh>
#include <Base/Common/BIASpragma.hh>
#include <MathAlgo/SVD.hh>
#include <Base/Debug/TimeMeasure.hh>

using namespace BIAS;
using namespace std;
 

UnVignette::UnVignette() 
{
  PrinciplePointX_ = 0.0;
  PrinciplePointY_ = 0.0;
  bIsInitialized_ = false;
  bUseLuT_ = false;
  bLeastSquaresPrepared_ = false;
  dInterPolationMethod_ = LEAST_SQUARES;
}

//if Constructor is called only with CameraParam LuT is not used!
UnVignette::UnVignette(CameraParam &cParam,INTERPOLATION_METHOD method, bool bUseLuT)
{
  dInterPolationMethod_  = method;
  bUseLuT_ = bUseLuT;
  bIsInitialized_ = false;
  bLeastSquaresPrepared_ = false;
  dWidth_ = cParam.GetWidth();
  dHeight_ = cParam.GetHeight();
    
  PrinciplePointX_ = cParam.GetPrincipalX()- dWidth_/2;
  PrinciplePointY_ = cParam.GetPrincipalY() - dHeight_/2;
  
  dWidthHalfMPPX_ = dWidth_/2 - PrinciplePointX_;
  dHeightHalfMPPY_ = dHeight_/2 - PrinciplePointY_; 
  
  ControlPoints_ = cParam.GetIlluCorrX();
  PercentageValues_ = cParam.GetIlluCorrY();
  if(ControlPoints_.size() > 2 && PercentageValues_.size() > 2)
    Init(cParam,method,bUseLuT_);
}

UnVignette::~UnVignette()
{
 
}


void UnVignette::Compute(BIAS::Image<unsigned char> &Image)
{
  if (!bIsInitialized_){
    BIASWARN("Unvignette is not initialised, doing nothing");
    return;
  }
  else{
    if(dWidth_ != Image.GetWidth() || dHeight_ != Image.GetHeight())
    {
      BIASERR("Image size does not fit init size, returning" );
      return;
    }
    float dum=0.0;
    unsigned int newPixelValue =0; 
    dChannelCount_ = Image.GetChannelCount();
   
    ///////////////////////////////////////////
    unsigned char* pId = Image.GetImageData();
    float* pLuId = LuImage_.GetImageData();

    for(unsigned int i=0;i<dWidth_*dHeight_;i++)
    { 
      dum = *pLuId++; // get data at look up image , this is correction factor
      for(unsigned int j=0;j<dChannelCount_;j++)      
      {
        newPixelValue = (unsigned int)( (*pId) * dum );
        if(newPixelValue <= 255)
           *pId++ = newPixelValue;
        else 
         *pId++ = 255;

      
      }
    }  
  }
  return;
}

int UnVignette::Init(const std::vector<double> ControlPoints,
                     std::vector<double> PercentageValues,
                     unsigned int dImageWidth , unsigned int dImageHeight, 
                     int dPPX, int dPPY, INTERPOLATION_METHOD method, bool bUseLuT)
{ 
  bUseLuT_ = bUseLuT; 
  dInterPolationMethod_ = method;
  dWidth_ = dImageWidth;
  dHeight_ = dImageHeight;
  PrinciplePointX_ = dPPX - dWidth_/2;
  PrinciplePointY_ = dPPY - dHeight_/2;
  ControlPoints_ = ControlPoints;
  PercentageValues_ = PercentageValues;
  return Init_();
}

int UnVignette::Init(BIAS::CameraParam cParam,INTERPOLATION_METHOD method, bool bUseLuT)
{
  bUseLuT_ = bUseLuT;
  dInterPolationMethod_ = method;
  dWidth_ = cParam.GetWidth();
  dHeight_ = cParam.GetHeight();
  PrinciplePointX_ = cParam.GetPrincipalX() - dWidth_/2;
  PrinciplePointY_ = cParam.GetPrincipalY() - dHeight_/2;
  ControlPoints_ = cParam.GetIlluCorrX();
  PercentageValues_ = cParam.GetIlluCorrY();
  return Init_();
}

int UnVignette::Init_(){
  int ret =0;
  dWidthHalfMPPX_ = dWidth_/2 - PrinciplePointX_;
  dHeightHalfMPPY_ = dHeight_/2 - PrinciplePointY_; 
  if(ControlPoints_.size() > 2 && PercentageValues_.size() > 2)  
  { 
    switch(dInterPolationMethod_){
      case SPLINE_WITHOUT_LUT:
      case SPLINE_WITH_LUT:
        InterpolatorIllu_.SetKnotPoints(ControlPoints_);  
        InterpolatorIllu_.SetControlPoints(PercentageValues_);
        InterpolatorIllu_.InitSpline(); 
   
        if(bUseLuT_)
          ret = InterpolatorIllu_.PrepareLuTSpline(ControlPoints_[0],
            ControlPoints_[ControlPoints_.size()-1],2000);
        
        break;
      case LEAST_SQUARES:
        ret = PrepareLeastSquaresMin_();
        break;
    }
    if(ret >=0)
      ret = PrepareLuImage_();
    if(ret >=0)
      bIsInitialized_ = true;
    return ret;
  }
  else
    return -1;
}


/////////////// Private functions /////////////////////
double UnVignette::
CalcRadius_(double x, double y)
{
  double radius=0.0;
  radius=sqrt((dWidthHalfMPPX_-x)*(dWidthHalfMPPX_-x) + (dHeightHalfMPPY_-y)*(dHeightHalfMPPY_-y));
  return radius;
}

int UnVignette::
PrepareLuImage_()
{
  // if use the look up image, prepare it!
  double radius=0.0;
  double dum = 0.0;
  //variables for calcradius
  dWidthHalfMPPX_ = dWidth_/2 - PrinciplePointX_;
  dHeightHalfMPPY_ = dHeight_/2 - PrinciplePointY_;
    
  LuImage_.Init(dWidth_,dHeight_,1);
  float ** ida = LuImage_.GetImageDataArray();

  for(unsigned int k=0;k<LuImage_.GetWidth();k++)
  {
    for(unsigned int l=0;l<LuImage_.GetHeight();l++)
    {
      radius = CalcRadius_(k,l);
      if(radius<0.0) radius = 0.0;
      switch(dInterPolationMethod_){
        case SPLINE_WITH_LUT:
        case SPLINE_WITHOUT_LUT :      
          if( bUseLuT_ )
              InterpolatorIllu_.SplineFromLuT(dum, radius/(dWidth_/2));
            else
              InterpolatorIllu_.Spline(dum,radius/(dWidth_/2),3);
        break;
        case LEAST_SQUARES:
           LeastSquares_(dum,radius/(dWidth_/2));
        break;
      } 
     
      // set factor in LuImage_ as correction factor;
      ida[l][k] = (float)dum;
    }
  }
  return 0;
}

///////////////// Prepare the Least Squares polynomial ///////////////
int UnVignette::
PrepareLeastSquaresMin_()
{
  int res =0;
  int rows = PercentageValues_.size();
  A_.newsize(rows,4);
  Y_.newsize(rows);
  for(int row=0;row<rows;row++)
  {
    
    A_[row][0] = 1;
    A_[row][1] = ControlPoints_[row];
    A_[row][2] = pow(ControlPoints_[row],2);
    A_[row][3] = pow(ControlPoints_[row],3);
    Y_[row] = PercentageValues_[row];
  }
  X_.newsize(4);
  
  BIAS::SVD svd;
  res=svd.Solve(A_,Y_,X_);
 // if(res!= 0)
  //  BIASERR("ERROR in SVD: "<<res);
   
  if(res == 0)
    bLeastSquaresPrepared_ = true;
  else
    bLeastSquaresPrepared_ = false;
  return res;
}

//////////////// if least squares is prepared, get value for radius //
int UnVignette::
LeastSquares_(double &res, double radius)
{
  if(!bLeastSquaresPrepared_)
  {  
    BIASERR("Least Squares is not correctly prepared, aborting.");
    return -1;
  }
  //X_ = {a,b,c,d}
  res = X_[3]*pow(radius,3) + X_[2]*pow(radius,2) +X_[1]*radius + X_[0];
  return 0;
}