;*****************************************************
;Name:		qi_hyper.pro
;Author:		Emmanuel Christophe
;Contact:		e.christophe at melaneum.com
;Description:	Compute QI on hyperspectral images
;Version:		v1.1 - 2005/06
;
;Regroup all necessary function to easily compute quality
;criteria for hyperspectral images. The expression of the 
;different quality criteria used can be found in:
;
;[1] E. Christophe, D. Léger, and C. Mailhes. Comparison 
;and evaluation of quality criteria for hyperspectral imagery.
;In Image Quality and System Performance II , volume 5668., 
;p.204-213, SPIE, Jan. 2005.
; 		 
;
;qi_hyper(image1, image2, criteria)
;     -> return QI value between image1 and image2 for the criteria number
;        0  -> mse
;        1  -> mad
;        2  -> pmad
;        3  -> snr
;        4  -> psnr
;        5  -> Qlambda
;        6  -> mss
;        7  -> rmse
;        8  -> rrmse
;        9  -> msam
;        10 -> sid
;        11 -> pearson
;        12 -> mae
;        13 -> Qspac
;        14 -> Qmult
;        15 -> Fid1
;        16 -> Fid2
;        17 -> Fid3
;qiname(criteria)
;     -> return string containing the QI name for the criteria number
;        the final string is 7 characters long
;qiextname(criteria)
;     -> return string containing the extended QI name for the criteria 
;        number: usefull for automatic title of graph for example
;qinametex(criteria)
;     -> return string containing the latex QI name for the criteria 
;        number: usefull for generating latex tables
;qinumber(criteria)
;     -> return number corresponding to the string criteria
;mse_hyper(image1, image2, nowarning=nowarning)
;     -> return double value of the MSE between image1 and image2
;        nowarning is to avoid output inside function using it (RMSE...)
;mad_hyper(image1, image2)
;     -> return value of the MAD between image1 and image2
;        Warning, value is same type as image1 and image2, call with
;        double(image1) for coherence
;pmad_hyper(image1, image2)
;     -> return value of the PMAD between image1 and image2
;        Warning, value is at least float, call with
;        double(image1) for coherence
;snr_hyper(image1, image2)
;     -> return double value of the SNR between image1 and image2
;psnr_hyper(image1, image2)
;     -> return double value of the PSNR between image1 and image2
;qlambda_hyper(image1, image2, posX=posX, posY=posY)
;     -> return double value of the Qlamda criteria between image1 and image2
;        (posX,posY) optionnaly give the position where the minimum Qlambda occurs
;mss_hyper(image1, image2, posX=posX, posY=posY)
;     -> return double value of the MSS criteria between image1 and image2
;        (posX,posY) optionnaly give the position where the maximum MSS occurs
;rmse_hyper(image1, image2)
;     -> return double value of the RMSE between image1 and image2
;rrmse_hyper(image1, image2)
;     -> return double value of the RRMSE between image1 and image2
;        Warning: no use of mse_hyper, not adapted yet to 1D signal
;msam_hyper(image1, image2, posX=posX, posY=posY)
;     -> return double value of the maximum SAM between image1 and image2
;        (posX,posY) optionnaly give the position where the maximum SAM occurs
;sid_hyper(image1, image2, posX=posX, posY=posY)
;     -> return double value of the maximum SID between image1 and image2
;        (posX,posY) optionnaly give the position where the maximum SID occurs
;pearson_hyper(image1, image2, posX=posX, posY=posY)
;     -> return double value of the minimum Pearson correlation between image1 and image2
;        (posX,posY) optionnaly give the position where the minimum occurs
;mae_hyper(image1, image2)
;     -> return double value of the MAE between image1 and image2
;qspac_hyper(image1, image2, posL=posL)
;     -> return double value of the Qlamda criteria between image1 and image2
;        posL optionnaly give the plane where the minimum Qlambda occurs
;fid1_hyper(image1, image2)
;     -> return double value of the fidelity criteria between image1 and image2
;        this one compute criteria directly on the whole hypercube
;fid2_hyper(image1, image2, posX=posX, posY=posY)
;     -> return double value of the fidelity criteria between image1 and image2
;        this one compute the criteria for each spectral pixel and return the min
;        (posX,posY) optionnaly give the pixel where the minimum fidelity occurs
;fid3_hyper(image1, image2, posL=posL)
;     -> return double value of the fidelity criteria between image1 and image2
;        this one compute the criteria for each image plane and return the min
;        (posL) optionnaly give the plane where the minimum fidelity occurs


function qiname, criteria
   CASE criteria OF
      0:  return, 'mse    '
      1:  return, 'mad    '
      2:  return, 'pmad   '
      3:  return, 'snr    '
      4:  return, 'psnr   '
      5:  return, 'Qlambda'
      6:  return, 'mss    '
      7:  return, 'rmse   '
      8:  return, 'rrmse  '
      9:  return, 'msam   '
      10: return, 'sid    '
      11: return, 'pearson'
      12: return, 'mae    '
      13: return, 'Qspac  '
      14: return, 'Qmult  '
      15: return, 'Fid1   '
      16: return, 'Fid2   '
      17: return, 'Fid3   '
      else: return, 'Unknown'
   ENDCASE
end

function qiextname, criteria
   CASE criteria OF
      0:  return, 'Mean Square Error'
      1:  return, 'Maximum Absolute Difference'
      2:  return, 'Percentage Maximum Absolute Difference'
      3:  return, 'Signal to Noise Ratio'
      4:  return, 'Peak Signal to Noise Ratio'
      5:  return, 'Spectral Q'
      6:  return, 'Maximum Spectral Similarity'
      7:  return, 'Root Mean Square Error'
      8:  return, 'Relative Root Mean Square Error'
      9:  return, 'Maximum Spectral Angle'
      10: return, 'Maximum Spectral Information Divergence'
      11: return, 'Minimum Pearson Correlation'
      12: return, 'Mean Absolute Error'
      13: return, 'Spatial Q'
      14: return, 'Combined Q'
      15: return, 'Fidelity'
      16: return, 'Spectral Fidelity'
      17: return, 'Spatial Fidelity'
      else: return, 'Unknown'
   ENDCASE
end

function qinametex, criteria
  CASE criteria OF
   0:  return, 'MSE    '
   1:  return, 'MAD    '
   2:  return, 'PMAD   '
   3:  return, 'SNR    '
   4:  return, 'PSNR   '
   5:  return, '$Q_\lambda$'
   6:  return, 'MSS    '
   7:  return, 'RMSE   '
   8:  return, 'RRMSE  '
   9:  return, 'MSAM   '
   10: return, 'SID    '
   11: return, 'Pearson'
   12: return, 'MAE    '
   13: return, '$Q_{(x,y)}$  '
   14: return, '$Q_m$  '
   15: return, '$F$    '
   16: return, '$F_\lambda$'
   17: return, '$F_{(x,y)}$'
   else: return, 'Unknown'
  ENDCASE
end

;critical order for checking!
;mse is in rmse, checking the longest first...
function qinumber, criteria
   if (strpos(criteria,'rrmse') NE -1) then return, 8
   if (strpos(criteria,'rmse') NE -1) then return, 7
   if (strpos(criteria,'mse') NE -1) then return, 0

   if (strpos(criteria,'pmad') NE -1) then return, 2
   if (strpos(criteria,'mad') NE -1) then return, 1

   if (strpos(criteria,'psnr') NE -1) then return, 4
   if (strpos(criteria,'snr') NE -1) then return, 3

   if (strpos(criteria,'Qlambda') NE -1) then return, 5

   if (strpos(criteria,'mss') NE -1) then return, 6

   if (strpos(criteria,'msam') NE -1) then return, 9
   if (strpos(criteria,'sid') NE -1) then return, 10
   if (strpos(criteria,'pearson') NE -1) then return, 11

   if (strpos(criteria,'mae') NE -1) then return, 12

   if (strpos(criteria,'Qspac') NE -1) then return, 13
   if (strpos(criteria,'Qmult') NE -1) then return, 14

   if (strpos(criteria,'Fid1') NE -1) then return, 15
   if (strpos(criteria,'Fid2') NE -1) then return, 16
   if (strpos(criteria,'Fid3') NE -1) then return, 17

   return, -1
end

;*****************************************************
;**** Compute MSE between two hyperspectral cubes ****
;*****************************************************

;Function is valid for 1, 2, or 3 dimensionnal signal
;TODO, put everything in one step

function mse_hyper, image1, image2, nowarning=nowarning

   if n_elements(nowarning) EQ 0 then nowarning=0
   warning=~nowarning

   sizeImg=size(image1)

   if (warning) then print, 'Calculating MSE...'

   time=systime(1)

   ;compute the MSE

   ;Hyperspectral Case
   if (sizeImg[0] eq 3) then begin
      mse = total(double(image1-image2)^2,/double) / (double(sizeImg[1])*double(sizeImg[2])*double(sizeImg[3]))
   ;   print, '3D'

   ;Ordinary image
   endif else if (sizeImg[0] eq 2) then begin
      mse = total(double(image1-image2)^2,/double) / (double(sizeImg[1])*double(sizeImg[2]))
   ;   print, '2D'

   ;Ordinary signal
   endif else if (sizeImg[0] eq 1) then begin
      mse = total(double(image1-image2)^2,/double) / (double(sizeImg[1]))
   ;   print, '1D'
   endif

   time=systime(1)-time
   if (warning) then begin
      print, "t (mse): ",time
      print, "mse : ", mse
   endif

   return, mse

end


;************************************************
;*** Compute MAD between two hyperspectral cubes*
;************************************************
function mad_hyper, image1, image2

   sizeImg=size(image1)
   print, 'Calculating MAD...'

   time=systime(1)

   mad=max(abs(image1-image2))

   time=systime(1)-time
   print, "t (mad): ",time
   print, "mad : ", mad

   return, mad

end

;******************************************************
;***   Compute PMAD between two hyperspectral cubes  **
;***           Percentage MAD                        **
;******************************************************

function pmad_hyper, image1, image2

   sizeImg=size(image1)
   print, 'Calculating PMAD...'

   time=systime(1)

   pmad=max(abs(image1-image2)/(image1+1./12))

   time=systime(1)-time
   print, "t (pmad): ",time

   return, pmad

end


;******************************************************
;***  Compute SNR between two hyperspectral cubes  ****
;******************************************************
function snr_hyper, image1, image2

   sizeImg=size(image1)
   print, 'Calculating SNR...'
   snr=0.0D

   time=systime(1)
   mse=mse_hyper(image1, image2, nowarning=1)
   snr=10*alog10(variance(image1)/(mse+1/12.))

   time=systime(1)-time
   print, "t (snr): ",time

   return, snr
end


;******************************************************
;***  Compute PSNR between two hyperspectral cubes  ***
;******************************************************
function psnr_hyper, image1, image2

   sizeImg=size(image1)
   print, 'Calculating PSNR...'
   psnr=0.0D


   time=systime(1)

   psnr=10*alog10(max(image1)^2./(mse_hyper(image1, image2, nowarning=1)+1/12.))

   time=systime(1)-time
   print, "t (psnr): ",time

   return, psnr

end



;********************************************
;**** Compute Qlambda between all spectra ***
;****  and return the minimum             ***
;********************************************
;see below for the spatial version

function qlambda_hyper, image1, image2, posX=posX, posY=posY

   print, 'Calculating spectral Qlambda...'

   sizeImg=size(image1,/dimension)

   time=systime(1)
   print, memory()
   image1=reform(image1,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)

   ;print, memory() 

   ;im=[temporary(image1),temporary(image2)]

   ;print, memory() 

   ;mean=total(im,2,/double)/sizeImg[2]
   ;mean1=mean[0:sizeImg[0]*sizeImg[1]-1]
   ;mean2=mean[sizeImg[0]*sizeImg[1]:2*sizeImg[0]*sizeImg[1]-1]
   mean1=total(image1,2,/double)/sizeImg[2]
   mean2=total(image2,2,/double)/sizeImg[2]

   ;print, memory()

   ;imcent=im-replicate(1,1,sizeImg[2])##mean
   ;var=total(double(imcent)^2,2,/double)/(sizeImg[2]-1)

   var1=total(double(image1-replicate(1,1,sizeImg[2])##mean1)^2,2,/double)/(sizeImg[2]-1)
   var2=total(double(image2-replicate(1,1,sizeImg[2])##mean2)^2,2,/double)/(sizeImg[2]-1)

;   imcent1=imcent[0:sizeImg[0]*sizeImg[1]-1,*]
;   imcent2=imcent[sizeImg[0]*sizeImg[1]:2*sizeImg[0]*sizeImg[1]-1,*]

   ;print, 'apres variance centrage...'
   ;print, memory()

;   var1=var[0:sizeImg[0]*sizeImg[1]-1]
;   var2=var[sizeImg[0]*sizeImg[1]:2*sizeImg[0]*sizeImg[1]-1]
   

   ;covar=total(imcent1*imcent2,2,/double)/((sizeImg[2]-1)  );   *sqrt(var1*var2))
   ;for memory optimization...
   covar=total( $
                (image1-replicate(1,1,sizeImg[2])##mean1)$
                 *(image2-replicate(1,1,sizeImg[2])##mean2),2,/double)$
                /((sizeImg[2]-1))

   ;print, memory()

;   qlambda=4*covar*mean1*mean2/((var1+var2)*(mean1^2+mean2^2))

   qlambda=4*covar*mean1*mean2 $ 
             /((var1+var2)*(mean1^2+mean2^2))

   ;print, memory()

   mqlambda=min(qlambda, pos)

   posX=pos MOD sizeImg[0]
   posY=pos /sizeImg[0]
   ;Just a checking
   ;qlambda=reform(qlambda,sizeImg[0],sizeImg[1])

   image1=reform(image1,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)

   print, "t (Qlambda): ", systime(1)-time
   return, mqlambda

end

;************************************************
;*** Compute Maximum Spectral Similarity (MSS) **
;*********between two hyperspectral cubes  ******
;************************************************

function mss_hyper, image1, image2, posX=posX, posY=posY

   ; posX  -> abciss of the maximum
   ; posY  -> ordinate of the maximum


   print, 'Calculating Maximum SS...'


   sizeImg=size(image1, /dimensions)

   mss=0.0D
   mse=dblarr(sizeImg[0],sizeImg[1])
   ncorr2=dblarr(sizeImg[0],sizeImg[1])
   ss=dblarr(sizeImg[0],sizeImg[1])


   time=systime(1)

   ;compute the MSS

   mse=total(double(image1-image2)^2,3,/double)/sizeImg[2]

   image1=reform(image1,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)

   ;mean of image1 and image2 on wavelength direction
   mean1=total(image1,2,/double)/sizeImg[2]
   mean2=total(image2,2,/double)/sizeImg[2]

   ;variance of image1 and image2 on wavelength direction
   tmp1=image1-replicate(1,1,sizeImg[2])##mean1
   tmp2=image2-replicate(1,1,sizeImg[2])##mean2
   ;Beware, after change, wavelength is now on 2 for total...
   var1=total(double(tmp1)^2,2,/double)/(sizeImg[2]-1) ; without bias
   var2=total(double(tmp2)^2,2,/double)/(sizeImg[2]-1)


   corr2= total((double(tmp1)*double(tmp2))/(sizeImg[2]-1),2,/double)^2 / (var1 * var2)

   ;Back on normal order
   ncorr2=reform((1-corr2)^2,sizeImg[0],sizeImg[1])

   ss=sqrt(mse+ncorr2)
   mss=max(ss,pos)

   posX=pos MOD sizeImg[0]
   posY=pos /sizeImg[0]

   image1=reform(image1,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)

   time=systime(1)-time
   print, "t (mss): ",time

   return, mss

end

;******************************************************
;***   Compute RMSE between two hyperspectral cubes  **
;******************************************************
function rmse_hyper, image1, image2

   print, 'Calculating RMSE...'

   time=systime(1)

   rmse=sqrt(mse_hyper(image1,image2,nowarning=1))

   time=systime(1)-time
   print, "t (rmse): ",time

   return, rmse

end


;*********************************************************
;* Compute Relative RMSE between two hyperspectral cubes *
;*********************************************************
function rrmse_hyper, image1, image2

   sizeImg=size(image1)
   print, 'Calculating Relative RMSE...'

   time=systime(1)

   ;compute the Relative RMSE

   ;Hyperspectral Case
   if (sizeImg[0] eq 3) then begin
      rrmse = sqrt(total((double(image1-image2)/(image1+1./12))^2,/double) / (double(sizeImg[1])*double(sizeImg[2])*double(sizeImg[3])))
      ;rrmse = sqrt(total((float(image1-image2)/(image1+1./12))^2) / (float(sizeImg[1])*float(sizeImg[2])*float(sizeImg[3])))

   ;Normal image
   endif else if (sizeImg[0] eq 2) then begin
      rrmse = sqrt(total((double(image1-image2)/(image1+1./12))^2,/double) / (double(sizeImg[1])*double(sizeImg[2])))
   endif



   time=systime(1)-time
   print, "t (rrmse): ",time

   return, rrmse

end

;**************************************************
;** Compute Maximum Spectral Angle Mapper (MSAM) **
;******  between two hyperspectral cubes  *********
;**************************************************


function msam_hyper, image1, image2, posX=posX, posY=posY
   print, 'Calculating Maximum SAM...'


   sizeImg=size(image1, /dimensions)

   msam=0.0D
   num=dblarr(sizeImg[0],sizeImg[1])
   den=dblarr(sizeImg[0],sizeImg[1])
   sam=dblarr(sizeImg[0],sizeImg[1])

   ;TODO check the size between image1 and image2
   time=systime(1)

   ;compute the MSAD
   num=double(total(double(image1)*double(image2),3))
   den=sqrt(total(double(image1)^2,3)*total(double(image2)^2,3))
   sam= acos(num/den)
   msam=max(sam,pos)


   posX=pos MOD sizeImg[0]
   posY=pos /sizeImg[0]
   ;print, 'max X ', posX
   ;print, 'max Y ', posY

   time=systime(1)-time
   print, "t (msam): ",time

   ;r=display_image(sam[*,*],"SAM")

   return, msam

end

;***************************************************
;** Compute Spectral Information Divergence (SID) **
;******  between two hyperspectral cubes  **********
;***************************************************

function sid_hyper, image1, image2, posX=posX, posY=posY

   sizeImg=size(image1,/dimensions)

   print, 'Calculating Maximum SID...'

   sid=dblarr(sizeImg[0],sizeImg[1])

   time=systime(1)

   image1=reform(image1,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)

   imp1=double(image1)/rebin(total(image1,2),sizeImg[0]*sizeImg[1],sizeImg[2])
   imp2=double(image2)/rebin(total(image2,2),sizeImg[0]*sizeImg[1],sizeImg[2])

   impnz1=imp1*(imp1 GT 0)+0.000001*(imp1 LT 0.000001)
   impnz2=imp2*(imp2 GT 0)+0.000001*(imp2 LT 0.000001)

   tmp=alog(impnz1/impnz2)

   impnz1=0
   impnz2=0

   d12=total(imp1*tmp,2,/double)
   d21=total(-imp2*tmp,2,/double)

   sid=d12+d21

   sid=reform(sid,sizeImg[0],sizeImg[1],/overwrite)

   time=systime(1)-time
   print, "t (sid): ",time

   msid= max(sid,pos)

   image1=reform(image1,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)

   posX=pos MOD sizeImg[0]
   posY=pos /sizeImg[0]

   return, msid
end



;********************************************************
;**** Compute Pearson correlation between all spectra ***
;****  and return the minimum correlation             ***
;********************************************************

function pearson_hyper, image1, image2, posX=posX, posY=posY

   print, 'Calculating Pearson correlation...'

   sizeImg=size(image1,/dimension)

   time=systime(1)
   ;Beware, after change, wavelength is now on dimension 2 for total...
   ;mean of image1 and image2 on wavelength direction
   ;variance of image1 on wavelength direction


   image1=reform(image1,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)
   im=[image1,image2]

   mean=total(im,2,/double)/sizeImg[2]
   ;imcent=im-rebin(mean,2*sizeImg[0]*sizeImg[1],sizeImg[2])
   imcent=im-replicate(1,1,sizeImg[2])##mean

   ;Warning, not really the variance here: normalization factor (cf later)
   var=total(double(imcent)^2,2,/double);/(sizeImg[2]-1)

   imcent1=imcent[0:sizeImg[0]*sizeImg[1]-1,*]
   imcent2=imcent[sizeImg[0]*sizeImg[1]:2*sizeImg[0]*sizeImg[1]-1,*]
   imcent=0
   var1=var[0:sizeImg[0]*sizeImg[1]-1,*]
   var2=var[sizeImg[0]*sizeImg[1]:2*sizeImg[0]*sizeImg[1]-1,*]
   var=0
   pearson=total(imcent1*imcent2,2,/double)/(sqrt(var1*var2));/(sizeImg[2]-1)
   ;Removing the normalization factor of the formula as we've done in the 
   ;computation of the variance (cf up)

   mpearson=min(pearson, pos)

   posX=pos MOD sizeImg[0]
   posY=pos /sizeImg[0]
 
   image1=reform(image1,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)

   print, "t (pearson): ", systime(1)-time
   return, mpearson

end

;*************************************************
;*** Compute MAE between two hyperspectral cubes *
;*********** MAE= Mean Absolute error ************
;*************************************************

function mae_hyper, image1, image2

   sizeImg=size(image1)
   print, 'Calculating MAE...'
   mae=0.0D

   time=systime(1)

   ;compute the MAE

   ;Hyperspectral Case
   if (sizeImg[0] eq 3) then begin
   mae = total(abs(image1-image2)) / (double(sizeImg[1])*double(sizeImg[2])*double(sizeImg[3]))

   ;Normal image
   endif else if (sizeImg[0] eq 2) then begin
   mae = total(abs(image1-image2), /double) / (double(sizeImg[1])*double(sizeImg[2]))
   endif


   time=systime(1)-time
   print, "t (mae): ",time

   return, mae

end

;********************************************
;**** Compute Qspac between all plane   ***
;****  and return the minimum             ***
;********************************************
;see above for spectral version

function qspac_hyper, image1, image2, posL=posL

   print, 'Calculating spacial Qlambda...'

   sizeImg=size(image1,/dimension)

   time=systime(1)

   image1=reform(image1,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0]*sizeImg[1],sizeImg[2],/overwrite)
   ;im=[[image1],[image2]]

   mean1=total(image1,1)/(sizeImg[0]*sizeImg[1])
   imcent1=image1-replicate(1,sizeImg[0]*sizeImg[1],1)##mean1
   var1=total(double(imcent1)^2,1,/double)/(sizeImg[0]*sizeImg[1]-1)
   ;var1=total(float(imcent1)^2,1)/(sizeImg[0]*sizeImg[1]-1)

   mean2=total(image2,1)/(sizeImg[0]*sizeImg[1])
   imcent2=image2-replicate(1,sizeImg[0]*sizeImg[1],1)##mean2

   covar=total(imcent1*imcent2,1)/((sizeImg[0]*sizeImg[1]-1)  ) 
   imcent1=0

   var2=total(double(imcent2)^2,1,/double)/(sizeImg[0]*sizeImg[1]-1)
   ;var2=total(float(imcent2)^2,1)/(sizeImg[0]*sizeImg[1]-1)
   imcent2=0


   ;mean=total(im,1,/double)/(sizeImg[0]*sizeImg[1])
   ;mean=total(im,1)/(sizeImg[0]*sizeImg[1])
   ;mean1=mean[0:sizeImg[2]-1]
   ;mean2=mean[sizeImg[2]:2*sizeImg[2]-1]

   ;imcent=im-replicate(1,sizeImg[0]*sizeImg[1],1)##mean
   ;var=total(double(imcent)^2,1,/double)/(sizeImg[0]*sizeImg[1]-1)
   ;mean=0
   ;var=total(float(imcent)^2,1)/(sizeImg[0]*sizeImg[1]-1)

   ;imcent1=imcent[*,0:sizeImg[2]-1]
   ;imcent2=imcent[*,sizeImg[2]:2*sizeImg[2]-1]

   ;var1=var[0:sizeImg[2]-1]
   ;var2=var[sizeImg[2]:2*sizeImg[2]-1]
   ;var=0
   ;division par les variances explicite ?
   ;covar=total(imcent1*imcent2,1,/double)/((sizeImg[0]*sizeImg[1]-1)  );   *sqrt(var1*var2))
   ;covar=total(imcent1*imcent2,1)/((sizeImg[0]*sizeImg[1]-1)  )

   ;mean1=total(image1,1)/(sizeImg[0]*sizeImg[1])
   ;mean2=total(image2,1)/(sizeImg[0]*sizeImg[1])

   qlambda=4*covar*mean1*mean2/((var1+var2)*(mean1^2+mean2^2))

   mean1=0
   mean2=0
   var1=0
   var2=0

   ;Do not consider negative values ! 
   indextoconsider=where(qlambda GT 0)

   mqlambda=min(qlambda[indextoconsider], posL)

   image1=reform(image1,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)
   image2=reform(image2,sizeImg[0],sizeImg[1],sizeImg[2],/overwrite)

   print, "t (Qspacial): ", systime(1)-time
   return, mqlambda
   ;print, "WARNING, experimental" & return, qlambda

end

;********************************************
;**** Compute Qlambda and Qspac           ***
;****  and multiply result                ***
;********************************************
;combine both spectral and spacial version
function qmult_hyper, image1, image2
   print, 'Calculating Q combinaison...'
   return, qspac_hyper(image1, image2)*qlambda_hyper(image1,image2)
end

;********************************************
;**** Compute fidelity on the whole cube  ***
;********************************************
function fid1_hyper, image1, image2
   print, 'Calculating fid1...'   
   fid =1 - (total(double(image1-image2)^2,/double) / total(double(image1)^2,/double))

   return, fid

end
;********************************************
;**** Compute fidelity for each pixel     ***
;****   and return the min                ***
;********************************************
function fid2_hyper, image1, image2, posX=posX, posY=posY
   print, 'Calculating fid2...'   
   sizeImg=size(image1)
   nc=double(sizeImg[1])
   nl=double(sizeImg[2])
   nb=double(sizeImg[3])

   ;nc=float(sizeImg[1])
   ;nl=float(sizeImg[2])
   ;nb=float(sizeImg[3])

   fid =replicate(1,nc,nl) - (total((double(image1)-image2)^2,3,/double) / total(double(image1)^2,3,/double))
   ;fid =replicate(1,nc,nl) - (total((float(image1)-image2)^2,3) / total(float(image1)^2,3))

   mfid=min(fid, pos)
   

   posX=pos MOD nc
   posY=pos / nc

   ;print,  'Warning'  
   ;return, fid
  
   return, mfid
   
end

;********************************************
;**** Compute fidelity for each plane     ***
;****   and return the min                ***
;********************************************
function fid3_hyper, image1, image2, posL=posL
   print, 'Calculating fid3...'   
   sizeImg=size(image1)
   nc=double(sizeImg[1])
   nl=double(sizeImg[2])
   nb=double(sizeImg[3])
   image1=reform(image1,nc*nl,nb,/overwrite)
   image2=reform(image2,nc*nl,nb,/overwrite)

   image1=reform(image1,nc,nl,nb,/overwrite)
   image2=reform(image2,nc,nl,nb,/overwrite)   
   
   fid =replicate(1,nb) - (total(double(image1-image2)^2,1,/double) / total(double(image1)^2,1,/double))

   mfid=min(fid, posL)
   
   return, mfid
   
end




function qi_hyper, image1, image2, criteria

;image1   -> image to evaluate
;image2   -> original image (for bivariant evaluation)
;criteria -> choosen criteria

;0  -> mse
;1  -> mad
;2  -> pmad
;3  -> snr
;4  -> psnr
;5  -> Qlambda
;6  -> mss
;7  -> rmse
;8  -> rrmse
;9  -> msam
;10 -> sid
;11 -> pearson
;12 -> mae
;13 -> Qspac
;14 -> Qmult
;15 -> fid1
;16 -> fid2
;17 -> fid3

CASE criteria OF
   0: return, mse_hyper(double(image1), double(image2))
   1: return, mad_hyper(double(image1), double(image2))
   2: return, pmad_hyper(double(image1), double(image2))
   3: return, snr_hyper(double(image1), double(image2))
   4: return, psnr_hyper(double(image1), double(image2))
   5: return, qlambda_hyper(double(image1), double(image2))
   6: return, mss_hyper(double(image1), double(image2))
   7: return, rmse_hyper(double(image1), double(image2))
   8: return, rrmse_hyper(double(image1), double(image2))
   9: return, msam_hyper(double(image1), double(image2))
   10: return, sid_hyper(double(image1), double(image2))
   11: return, pearson_hyper(double(image1), double(image2))
   12: return, mae_hyper(double(image1),double(image2))
   13: return, qspac_hyper(double(image1), double(image2))
   14: return, qmult_hyper(double(image1), double(image2))
   15: return, fid1_hyper(double(image1), double(image2))
   16: return, fid2_hyper(double(image1), double(image2))
   17: return, fid3_hyper(double(image1), double(image2))

   else : return, -1
ENDCASE
end