ExampleHdr.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
*/


/**
 * @example ExampleHdr.cpp
   @deprecated Schedule for removal!!
   @relates SVD
   @brief Example for Shit
   @ingroup g_examples
   @author  MIP
*/

#include <Base/Image/Image.hh>
#include <Base/Math/Matrix.hh>
#include <Base/Math/Vector.hh>
#include <MathAlgo/SVD.hh>

using namespace std;
using namespace BIAS;


//////////////////////////////////
//////////////////////////////////
//                              //
//  THIS IS UNDER CONSTRUCTION  //
//                              //
//////////////////////////////////
//////////////////////////////////




int main(int argc, char *argv[]) {

  return 0;
}

//
// gsolve.m Solve for imaging system response function
//
// Given a set of pixel values observed for several pixels in several
// images with different exposure times, this function returns the
// imaging systems response function g as well as the log film irradiance
// values for the observed pixels.
//
// Assumes:
//
//   Zmin = 0
//   Zmax = 255
//
// Arguments:
//
//   Z(i,j)   is the pixel values of pixel location number i in image j
//   B(j)     is the log delta t, or log shutter speed, for image j
//   l        is lamdba, the constant that determines the amount of smoothness
//   w(z)     is the weighting function value for pixel value z
//
// Returns:
//
//   g(z)     is the log exposure corresponding to pixel value z
//   lE(i)    is the log film irradiance at pixel location i
//


void gsolve(BIAS::Matrix<unsigned char> Z, BIAS::Vector<double> B, double l,
            BIAS::Vector<int> w) {
  const int n = sizeof(unsigned char);

  // prepare matrix for linear equation system
  BIAS::Matrix<double> A(Z.num_rows() * Z.num_cols() + n + 1,
                        n + Z.num_rows(),
                        0.0);

  // prepare vector for linear equation system
  BIAS::Vector<double> b(A.num_rows(), 1.0);

  // include the data-fitting equations
  int k = 0;
  double wij; // temp value
  for (int i = 0; i < Z.num_rows(); i++) {
    for (int j = 0; j < Z.num_cols(); j++) {
      wij = w[Z[i][j] + 1]; // +1 ???
      A[k][Z[i][j] + 1] = wij; // +1 ???
      A[k][n + i] = -wij;
      b[k] = wij * B[j];
      k++;
    }
  }

  // Fix the curve by setting its middle value to 0
  A[k][128] = 1;
  k++;

  //Include the smoothness equations
  for (int i = 0; i < n - 2; i++) {
    A[k][i] = l * w[i+1]; // +1 ???
    A[k][i+1] = -2 * l * w[i+1]; // +1 ???
    A[k][i+2] = l * w[i+1]; // +1 ???
    k++;
  }

  // Solve the system using SVD
  SVD svd;
  BIAS::Vector<double> x;
  svd.Solve(A, b, x);

  // get cam response function and irradiance from solution vector
  BIAS::Vector<double> g(n);
  BIAS::Vector<double> lE(x.size() - n);

  for (int i = 0; i < n; i++) {
    g[i] = x[i];
  }

  for (int i = n; i < x.size(); i++) {
    lE[i - n] = x[i];
  }

  // TODO: return g and lE
}



/*
function [g,lE]=gsolve(Z,B,l,w)
//n = 256;
//A = zeros(size(Z,1)*size(Z,2)+n+1,n+size(Z,1));
//b = zeros(size(A,1),1);
%% Include the data fitting equations
k = 1;
for i=1:size(Z,1)
  for j=1:size(Z,2)
      wij = w(Z(i,j)+1);
      A(k,Z(i,j)+1) = wij;
      A(k,n+i) =  wij;
      b(k,1) = wij * B(i,j); %%%%%%%%%%%% <-------- may be wrong! B is a vector
      k=k+1;
  end
end
%% Fix the curve by setting its middle value to 0
A(k,129) = 1;
k=k+1;
%% Include the smoothness equations
for i=1:n2
  A(k,i)=l*w(i+1);
  A(k,i+1)=2*l*w(i+1);
  A(k,i+2)=l*w(i+1);
  k=k+1;
end
%% Solve the system using SVD
x = A\b;
g = x(1:n);
lE = x(n+1:size(x,1));
*/