//------------------------- sharpness_index -------------------------

//    Compute the sharpness index (SI) 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 Sharpness Index (SI) described in the paper
//
//  G. Blanchet, L. Moisan, An explicit Sharpness Index related to
//  Global Phase Coherence, proceedings of the IEEE International Conference 
//  on Acoustics, Speech, and Signal Processing (ICASSP), pp. 1065-1068, 2012. 
//
//  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:    si = sharpness_index(u)   or   si = sharpness_index(u,pmode) 
//
// u is an image (a 2D array = a Scilab matrix)
//
// available preprocessing modes are:
//
// pmode = 0             raw sharpness index of u
// pmode = 1             sharpness index of the periodic component of u 
// pmode = 2             sharpness index of the 1/2,1/2-translated of u 
// pmode = 3 (default)   sharpness 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 (09/2013): initial Scilab version
// v1.1 (01/2014): added preprocessing modes
// v1.2 (01/2014): avoid division by zero for degenerate images 
// v1.3 (01/2014): improved erfc estimate 
// v1.4 (02/2014): exact (up to machine precision) logerfc function

function si = sharpness_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; fgx = fft(gx)
gy = u([2:ny,1],:)-u; fgy = fft(gy)
tv = sum(abs(gx)+abs(gy))
Gxx = real(ifft(fgx.*conj(fgx)))
Gyy = real(ifft(fgy.*conj(fgy)))
Gxy = real(ifft(fgx.*conj(fgy)))
deff("y=oomega(t)","y=real(t.*asin(t)+sqrt(1-t.^2)-1)")
var = 0
axx = Gxx(1,1);      if axx>0, var = var+  axx*sum(oomega(Gxx/axx)); end
ayy = Gyy(1,1);      if ayy>0, var = var+  ayy*sum(oomega(Gyy/ayy)); end
axy = sqrt(axx*ayy); if axy>0, var = var+2*axy*sum(oomega(Gxy/axy)); end
var = var*2/%pi
if var>0 
  // t = ( E(TV)-tv )/sqrt(var(TV))
  t = ( (sqrt(axx)+sqrt(ayy))*sqrt(2*nx*ny/%pi) - tv )/sqrt(var)
  // si = -log10(P(N(0,1)>t))
  si = -logerfc(t/sqrt(2))/log(10)+log10(2)
else 
  si = 0
end
endfunction
