//---------------------------- s_index ----------------------------

// Compute the simplified sharpness index (S) of a numerical image
//
//               author: Lionel Moisan
//
//  This program is freely available on the web page
//
//     http://www.mi.parisdescartes.fr/~moisan/sharpness/
//
//  It implements the S Index described in the paper
//
//  A. Leclaire, L. Moisan, No-reference image quality assessment 
//  and blind deblurring with sharpness metrics exploiting Fourier 
//  phase information, preprint MAP5 n° 2014-03, 2014.
//
//  If you use it for a publication, please mention this paper
//  If you modify this program, please indicate in the code that you 
//  did so and leave this message.
//  You can report bugs or suggestions to lionel.moisan [AT] parisdescartes.fr
//
// usage:    s = s_index(u)   or   s = s_index(u,pmode) 
//
// u is an image (a 2D array = a Scilab matrix)
//
// available preprocessing modes are:
//
// pmode = 0             raw S index of u
// pmode = 1             S index of the periodic component of u 
// pmode = 2             S index of the 1/2,1/2-translated of u 
// pmode = 3 (default)   S index of the 1/2,1/2-translated 
//                       of the periodic component of u
//
// Default mode (pmode = 3) should be used, unless you want to work on very
// specific images that are naturally periodic or not quantized (see paper)
//
// note: this function also works for 1D signals (line or column vectors)
//
// v1.0 (01/2014): initial version from sharpness_index.sci v1.0
// v1.1 (01/2014): improved erfc estimate 
// v1.2 (02/2014): exact (up to machine precision) logerfc function

function s = s_index(u,pmode)

// log of erfc function 
// (special code that avoids underflow for large values of x)
function y = logerfc(x)
  if x<20 
    y = log(erfc(x))
  else
    z = x.^(-2)
    s = 1
    for k=8:-1:1
      s = 1-(k-0.5)*z.*s
    end
    y = -0.5*log(%pi)-x.^2+log(s./x)
  end
endfunction

// Dequantization with (1/2,1/2) translation in Fourier domain 
function v = dequant(u)
  [ny,nx] = size(u)
  mx = int(nx/2); my = int(ny/2)
  Tx = exp(-%i*%pi/nx*( modulo(mx:mx+nx-1,nx)-mx ))
  Ty = exp(-%i*%pi/ny*( modulo(my:my+ny-1,ny)-my ))
  v = real(ifft(fft(u).*(Ty.'*Tx))) 
endfunction
 
// Periodic plus Smooth Image Decomposition
// see http://www.mi.parisdescartes.fr/~moisan/p+s
function [p,s] = perdecomp(u)
  [ny,nx] = size(u); X = 1:nx; Y = 1:ny
  v = zeros(u)
  v(1,X)  = u(1,X)-u(ny,X)
  v(ny,X) = -v(1,X)
  v(Y,1 ) = v(Y,1 )+u(Y,1)-u(Y,nx)
  v(Y,nx) = v(Y,nx)-u(Y,1)+u(Y,nx)
  fx = ones(ny,1)*cos(2.*%pi*(X-1)/nx)
  fy = cos(2.*%pi*(Y'-1)/ny)*ones(1,nx)
  fx(1,1)=0.   // avoid division by 0 in the line below
  s = real(ifft(fft(v)*0.5./(2.-fx-fy)))
  p = u-s
endfunction


// MAIN FUNCTION

if argn(2)<2, pmode = 3; end
if pmode==1 | pmode==3, u = perdecomp(u); end
if pmode==2 | pmode==3, u = dequant(u); end
[ny,nx] = size(u)
gx = u(:,[2:nx,1])-u
gy = u([2:ny,1],:)-u
tv = sum(abs(gx)+abs(gy))
fu = fft(u)
p = 0:nx-1; q = 0:ny-1
P = 2*%pi*ones(q')*p/nx; Q = 2*%pi*q'*ones(p)/ny
fgx2 = fu.*sin(P/2); fgx2 = real(4*fgx2.*conj(fgx2))
fgy2 = fu.*sin(Q/2); fgy2 = real(4*fgy2.*conj(fgy2))
fgxx2 = sum(fgx2.^2)
fgyy2 = sum(fgy2.^2)
fgxy2 = sum(fgx2.*fgy2)
vara = 0 // the simplified variance
axx = sum(gx.^2);    if axx>0, vara = vara+fgxx2/axx; end
ayy = sum(gy.^2);    if ayy>0, vara = vara+fgyy2/ayy; end
axy = sqrt(axx*ayy); if axy>0, vara = vara+2*fgxy2/axy; end
vara = vara/(%pi*nx*ny)
if vara>0 
  // t = ( E(TV)-tv )/sqrt(vara))
  t = ( (sqrt(axx)+sqrt(ayy))*sqrt(2*nx*ny/%pi) - tv )/sqrt(vara)
  // s = -log10(P(N(0,1)>t))
  s = -logerfc(t/sqrt(2))/log(10)+log10(2)
else 
  s = 0
end
endfunction
