#include <stdio.h>
#include <stdlib.h>

/*-- get the nx064 array definition --*/

#include "AlphaTab64.h"  /*-- auto-generated by script Alpha/malper64 --*/

/*---------------------------------------------------------------------------*/
/*! Function to indicate cluster significance level alpha:
     - csize = cluster size
     - nz    = # slices
     - fw    = FWHM ( pixels)
     - pv    = p-value per voxel

   Return value is
     -  -1  ==> unavailable
     -   0  ==> alpha < 0.01
     -   1  ==> alpha < 0.05
     -   2  ==> alpha < 0.10
     - 666  ==> alpha > 0.10
*//*-------------------------------------------------------------------------*/

int cluster_alphaindex_64( int csize , int nz , float fw , float pv )
{
   int iz , ja , kf , lp ;
   float        fkf ,flp , ccpp,ccpm,ccmp,ccmm , cc ;

   if( csize < 2 || nz < 1 || nz > NZMAX ) return -1 ;   /** bad inputs **/

        if( fw < 0.0f          ) return -1 ;
   else if( fw > fwhm[NFWHM-1] ) return -1 ;             /** bad input fw **/

        if( pv < pval[0]       ) pv = pval[0] ;
   else if( pv > pval[NPVAL-1] ) return -1 ;             /** bad input pv */

   iz = nz-1 ;  /* slice index in nx064 array*/

   /* find FWHM index so that fwhm[kf] >= fw >= fwhm[kf-1] */

  for( kf=1 ; kf < NFWHM && fw > fwhm[kf] ; kf++ ) ; /*nada */
  if( kf == NFWHM ) return -1 ;
  fkf = (fwhm[kf]-fw) / (fwhm[kf]-fwhm[kf-1]) ;  /* fractional step in kf */

  /* find pval index so that pval[lp] >= pv >= pval[lp-1] */

  for( lp=1 ; lp < NPVAL && pv > pval[lp] ; lp++ ) ; /*nada*/
  if( lp == NPVAL ) return -1 ;
  flp = (pval[lp]-pv) / (pval[lp]-pval[lp-1]) ;  /* fractional step in lp */

  /* scan different alpha levels to see if this cluster size is significant */

  for( ja=0 ; ja < NALPH ; ja++ ){
    ccpp = nx064[iz][ja][kf  ][lp  ] ;     /* bi-linearly interpolate in    */
    ccpm = nx064[iz][ja][kf  ][lp-1] ;     /* the nx064 table in fw and pv  */
    ccmp = nx064[iz][ja][kf-1][lp  ] ;     /* to find the critical cluster  */
    ccmm = nx064[iz][ja][kf-1][lp-1] ;     /* size at the ja-th alpha level */
    cc =  (flp*ccmm+(1.0f-flp)*ccmp)*fkf
        + (flp*ccpm+(1.0f-flp)*ccpp)*(1.0f-fkf) ;
    if( csize >= cc ) return ja ;          /* significance at this alpha */
  }
  return 666 ;  /* not significant */
}
