/*
 * John Ashburner
 * Copyright (C) 1997-2022 Wellcome Centre for Human Neuroimaging
 */

/* MATLAB independent image access routines
 * Gateway to routines compiled from spm_vol_utils.c 
 */

#include <math.h>
#include <stdio.h>
#include "spm_vol_access.h"
#include "spm_datatypes.h"
#ifdef SPM_WIN32
#include <windows.h>
#endif

int get_datasize(int type)
{
    if (type == SPM_UNSIGNED_CHAR      || type == SPM_SIGNED_CHAR)          return(8);
    if (type == SPM_SIGNED_SHORT       || type == SPM_SIGNED_SHORT_S)       return(16);
    if (type == SPM_UNSIGNED_SHORT     || type == SPM_UNSIGNED_SHORT_S)     return(16);
    if (type == SPM_SIGNED_INT         || type == SPM_SIGNED_INT_S)         return(32);
    if (type == SPM_UNSIGNED_INT       || type == SPM_UNSIGNED_INT_S)       return(32);
    if (type == SPM_UNSIGNED_LONG_LONG || type == SPM_UNSIGNED_LONG_LONG_S) return(64);
    if (type == SPM_SIGNED_LONG_LONG   || type == SPM_SIGNED_LONG_LONG_S)   return(64);
    if (type == SPM_FLOAT              || type == SPM_FLOAT_S)              return(32);
    if (type == SPM_DOUBLE             || type == SPM_DOUBLE_S)             return(64);
    return(0);
}

int resample(int m, MAPTYPE *vol, double *out, double *x, double *y, double *z, int hold, double background)
{
    extern void resample_uchar(   int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_schar(   int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_short(   int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_ushort(  int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_int(     int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_uint(    int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_int64(   int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_uint64(  int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_float(   int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_double(  int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_short_s( int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_ushort_s(int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_int_s(   int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_uint_s(  int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_int64_s( int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_uint64_s(int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_float_s( int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_double_s(int,void**,double*,double*,double*,double*,int,int,int,int,double,double*,double*);

#ifdef _MSC_VER
    /* https://msdn.microsoft.com/en-us/library/windows/desktop/aa366801.aspx */
    __try
    {
#endif
    if (vol->dtype == SPM_UNSIGNED_CHAR)
         resample_uchar(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_CHAR)
         resample_schar(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_SHORT)
         resample_short(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_SHORT)
         resample_ushort(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_INT)
           resample_int(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_INT)
           resample_uint(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_FLOAT)
         resample_float(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_DOUBLE)
        resample_double(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_SHORT_S)
        resample_short_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_SHORT_S)
        resample_ushort_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_INT_S)
        resample_int_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_INT_S)
        resample_uint_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_FLOAT_S)
        resample_float_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_DOUBLE_S)
        resample_double_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_LONG_LONG)
        resample_int64(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_LONG_LONG)
        resample_uint64(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_LONG_LONG_S)
        resample_int64_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_LONG_LONG_S)
        resample_uint64_s(m,vol->data,out,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else
    {
        (void)fprintf(stderr,"%d: Unknown datatype.\n", vol->dtype);
        return(1);
    }
#ifdef _MSC_VER
    }
    __except(GetExceptionCode()==EXCEPTION_IN_PAGE_ERROR ?
        EXCEPTION_EXECUTE_HANDLER : EXCEPTION_CONTINUE_SEARCH)
    {
        (void)fprintf(stderr,"An exception occurred while accessing the data.");
        return(1);
    }
#endif
    return(0);
}

int resample_d(int m, MAPTYPE *vol, double *out, double *gradx, double *grady, double *gradz, double *x, double *y, double *z, int hold, double background)
{
    extern void resample_d_uchar(   int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_schar(   int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_short(   int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_ushort(  int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_int(     int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_uint(    int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_float(   int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_double(  int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_uint64(  int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_int64(   int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_short_s( int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_ushort_s(int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_int_s(   int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_uint_s(  int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_float_s( int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_double_s(int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_uint64_s(int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);
    extern void resample_d_int64_s( int,void**,double*,double*,double*,double*,double*,double*,double*,int,int,int,int,double,double*,double*);

#ifdef _MSC_VER
    /* https://msdn.microsoft.com/en-us/library/windows/desktop/aa366801.aspx */
    __try
    {
#endif
    if (vol->dtype == SPM_UNSIGNED_CHAR)
         resample_d_uchar(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_CHAR)
         resample_d_schar(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_SHORT)
         resample_d_short(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_SHORT)
         resample_d_ushort(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_INT)
           resample_d_int(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_INT)
           resample_d_uint(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_FLOAT)
         resample_d_float(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_DOUBLE)
        resample_d_double(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_SHORT_S)
        resample_d_short_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_SHORT_S)
        resample_d_ushort_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_INT_S)
        resample_d_int_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_INT_S)
        resample_d_uint_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_FLOAT_S)
        resample_d_float_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_DOUBLE_S)
        resample_d_double_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_LONG_LONG)
        resample_d_uint64(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_LONG_LONG)
        resample_d_int64(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_LONG_LONG_S)
        resample_d_uint64_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_LONG_LONG_S)
        resample_d_int64_s(m,vol->data,out,gradx,grady,gradz,x,y,z,vol->dim[0],vol->dim[1],vol->dim[2],
            hold, background, vol->scale,vol->offset);
    else
    {
        (void)fprintf(stderr,"%d: Unknown datatype.\n", vol->dtype);
        return(1);
    }
#ifdef _MSC_VER
    }
    __except(GetExceptionCode()==EXCEPTION_IN_PAGE_ERROR ?
        EXCEPTION_EXECUTE_HANDLER : EXCEPTION_CONTINUE_SEARCH)
    {
        (void)fprintf(stderr,"An exception occurred while accessing the data.");
        return(1);
    }
#endif
    return(0);
}

int slice(double *mat, double *image, int xdim1, int ydim1, MAPTYPE *vol, int hold, double background)
{
    extern int slice_uchar(   double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_schar(   double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_short(   double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_ushort(  double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_int(     double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_uint(    double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_float(   double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_double(  double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_uint64(  double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_int64(   double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_short_s( double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_ushort_s(double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_int_s(   double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_uint_s(  double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_float_s( double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_double_s(double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_uint64_s(double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
    extern int slice_int64_s( double *, double *, int, int, void **, int, int, int, int, double, double*, double *);
 
    int sts = 1;
#ifdef _MSC_VER
    /* https://msdn.microsoft.com/en-us/library/windows/desktop/aa366801.aspx */
    __try
    {
#endif
    if (vol->dtype == SPM_UNSIGNED_CHAR)
         sts = slice_uchar(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_CHAR)
         sts = slice_schar(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_SHORT)
         sts = slice_short(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_SHORT)
         sts = slice_ushort(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_INT)
         sts = slice_int(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_INT)
         sts = slice_uint(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_FLOAT)
        sts = slice_float(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_DOUBLE)
        sts = slice_double(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_SHORT_S)
        sts = slice_short_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2], 
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_SHORT_S)
        sts = slice_ushort_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2], 
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_INT_S)
        sts = slice_int_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2], 
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_INT_S)
        sts = slice_uint_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2], 
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_FLOAT_S)
        sts = slice_float_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_DOUBLE_S)
        sts = slice_double_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2], 
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_LONG_LONG)
        sts = slice_uint64(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_LONG_LONG)
        sts = slice_int64(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_UNSIGNED_LONG_LONG_S)
        sts = slice_uint64_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else if (vol->dtype == SPM_SIGNED_LONG_LONG_S)
        sts = slice_int64_s(mat, image, xdim1,ydim1, vol->data, vol->dim[0],vol->dim[1],vol->dim[2],
            hold,background, vol->scale,vol->offset);
    else
    {
        (void)fprintf(stderr,"%d: Unknown datatype.\n", vol->dtype);
        sts = 1;
    }
#ifdef _MSC_VER
    }
    __except(GetExceptionCode()==EXCEPTION_IN_PAGE_ERROR ?
        EXCEPTION_EXECUTE_HANDLER : EXCEPTION_CONTINUE_SEARCH)
    {
        (void)fprintf(stderr,"An exception occurred while accessing the data.");
        sts = 1;
    }
#endif
    return(sts);
}
