//----------------------- global_phase_coherence -----------------------
//
//    Compute the Global Phase Coherence (GPC) 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, B. Rougé. Measuring the global phase
//  coherence of an image, proceedings of the IEEE International Conference 
//  on Image Processing (ICIP), pp. 1176-1179, 2008.
//
//  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:       gpc = global_phase_coherence(u)   
//          or  gpc = global_phase_coherence(u,N)
//          or  gpc = global_phase_coherence(u,N,pmode)
//
// u is an image (a 2D array = a Scilab matrix)
//
// N is the number of Monte-Carlo iterations used to estimate gpc
// (default value: 100)
//
// available preprocessing modes are:
//
// pmode = 0             raw global phase coherence of u
// pmode = 1             global phase coherence of the periodic component of u 
// pmode = 2             global phase coherence of the 1/2,1/2-translated of u 
// pmode = 3 (default)   global phase coherence 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)
//
// Remarks: 
//   - this function also works for 1D signals (line or column vectors)
//   - the convergence is very slow (run the function several times to
//       have an idea of the precision attained for the value of N used)
//
// v1.0 (01/2012): initial Scilab version
// v1.1 (01/2014): added preprocessing modes
// v1.2 (01/2014): improved erfc estimate 
// v1.3 (02/2014): exact (up to machine precision) logerfc function

function gpc = global_phase_coherence(u,N,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, N = 100; end
if argn(2)<3, pmode = 3; end
if pmode==1 | pmode==3, u = perdecomp(u); end
if pmode==2 | pmode==3, u = dequant(u); end
// initialization for random phases
[ny,nx] = size(u)
fu = fft(u)
x = 1:nx; y = 1:ny; X = ones(y')*x; Y = y'*ones(x)
Xp = modulo(nx+1-X,nx)+1; Yp = modulo(ny+1-Y,ny)+1
i = (X-1)*ny+Y; ip = (Xp-1)*ny+Yp
A = find(ip>i)
B = find(ip==i)
phi = zeros(ny,nx)
// compute Total Variation of u
tvu = sum( abs(u(:,[2:nx,1])-u)+abs(u([2:ny,1],:)-u))
// MAIN LOOP (Monte-Carlo samples)
tv = zeros(1,N)
for i=1:N  
  // build random phases with Hermitian symmetry
  phi(A) = grand(A,"unf",0,2*%pi)
  phi(ip(A)) = -phi(A)
  phi(B) = %pi*grand(B,"uin",0,1)
  // compute tv of random phase image
  v = real(ifft(fu.*exp(%i*phi)))
  tv(i) = sum( abs(v(:,[2:nx,1])-v)+abs(v([2:ny,1],:)-v))
end
// estimate first and second order moments
mu = sum(tv)/N
var = sum((tv-mu).^2/(N-1))
// t is the estimate of ( E(TV)-tv )/sqrt(var(TV))
t = (mu-tvu)/sqrt(var)
// gpc = -log10(P(N(0,1)>t))
gpc = -logerfc(t/sqrt(2))/log(10)+log10(2)
endfunction
