/*********************************************************************
 *
 * Very simple code snippets to read/write nifti1 files
 * This code is placed in the public domain.
 *
 * If you are the type who doesn't want to use a file format unless
 * you can write your own i/o code in less than 30minutes, this
 * example is for you.
 *
 * This code does not deal with wrong-endian data, compressed data,
 * the new qform/sform orientation codes, parsing filenames, volume-
 * wise or timecourse-wise data access or any of a million other very useful
 * things that are in the niftilib i/o reference libraries.
 * We encourage people to use the niftilib reference library and send
 * feedback/suggestions, see http://niftilib.sourceforge.net/
 * But, if that is too much to tackle and you just want to jump in, this
 * code is a starting point.
 * This code was written for maximum readability, not for the greatest
 * coding style.
 *
 *
 * If you are already a little familiar with reading/writing Analyze
 * files of some flavor, and maybe even have some of your own code, here
 * are the most important things to be aware of in transitioning to nifti1:
 *
 * 1. nii vs .hdr/.img
 *      nifti1 datasets can be stored either in .hdr/.img pairs of files
 *      or in 1 .nii file.  In a .nii file the data will start at the byte
 *      specified by the vox_offset field, which will be 352 if no extensions
 *      have been added.  And, nifti1 really does like that magic field set
 *      to "n+1" for .nii and "ni1" for .img/.hdr
 *
 * 2. scaling
 *      nifti1 datasets can contain a scaling factor.  You need to check the
 *      scl_slope field and if that isn't 0, scale your data by
 *      Y * scl_slope  + scl_inter
 *
 * 3. extensions
 *      nifti1 datasets can have some "extension data" stuffed after the
 *      regular header.  You can just ignore it, but, be aware that a
 *      .hdr file may be longer than 348 bytes, and, in a .nii file
 *      you can't just jump to byte 352, you need to use the vox_offset
 *      field to get the start of the image data.
 *
 * 4. new datatypes
 *      nifti1 added a few new datatypes that were not in the Analyze 7.5
 *      format from which nifti1 is derived.  If you're just working with
 *      your own data this is not an issue but if you get a foreign nifti1
 *      file, be aware of exotic datatypes like DT_COMPLEX256 and mundane
 *      things like DT_UINT16.
 *
 * 5. other stuff
 *     nifti1 really does like the dim[0] field set to the number of
 *     dimensions of the dataset.  Other Analyze flavors might not
 *     have been so scrupulous about that.
 *     nifti1 has a bunch of other new fields such as intent codes,
 *     qform/sform, etc. but, if you just want to get your hands on
 *     the data blob you can ignore these.  Example use of these fields
 *     is in the niftilib reference libraries.
 *
 *
 *
 * To compile:
 * You need to put a copy of the nifti1.h header file in this directory.
 * It can be obtained from the NIFTI homepage  http://nifti.nimh.nih.gov/
 * or from the niftilib SourceForge site http://niftilib.sourceforge.net/
 *
 * cc -o nifti1_read_write nifti1_read_write.c
 *
 *
 * To run:
 * nifti1_read_write -w abc.nii abc.nii
 * nifti1_read_write -r abc.nii abc.nii
 *
 *
 * The read method is hardcoded to read float32 data.  To change
 * to your datatype, just change the line:
 * typedef float MY_DATATYPE;
 *
 * The write method is hardcoded to write float32 data.  To change
 * to your datatype, change the line:
 * typedef float MY_DATATYPE;
 * and change the lines:
 * hdr.datatype = NIFTI_TYPE_FLOAT32;
 * hdr.bitpix = 32;
 *
 *
 * Written by Kate Fissell, University of Pittsburgh, May 2005.
 *
 *********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "nifti1.h"

typedef float MY_DATATYPE;

#define MIN_HEADER_SIZE 348
#define NII_HEADER_SIZE 352


/**********************************************************************
 *
 * read_nifti_file
 *
 **********************************************************************/
static int read_nifti_file(char * const hdr_file, char * const data_file)
{
        nifti_1_header hdr;
        FILE *fp;
        int ret,i;
        double total;
        MY_DATATYPE *data=NULL;


/********** open and read header */
        fp = fopen(hdr_file,"r");
        if (fp == NULL) {
                fprintf(stderr, "\nError opening header file %s\n",hdr_file);
                exit(1);
        }
        ret = fread(&hdr, MIN_HEADER_SIZE, 1, fp);
        if (ret != 1) {
                fprintf(stderr, "\nError reading header file %s\n",hdr_file);
                exit(1);
        }
        fclose(fp);


/********** print a little header information */
        fprintf(stderr, "\n%s header information:",hdr_file);
        fprintf(stderr, "\nXYZT dimensions: %d %d %d %d",hdr.dim[1],hdr.dim[2],hdr.dim[3],hdr.dim[4]);
        fprintf(stderr, "\nDatatype code and bits/pixel: %d %d",hdr.datatype,hdr.bitpix);
        fprintf(stderr, "\nScaling slope and intercept: %.6f %.6f",hdr.scl_slope,hdr.scl_inter);
        fprintf(stderr, "\nByte offset to data in datafile: %ld",(long)(hdr.vox_offset));
        fprintf(stderr, "\n");


/********** open the datafile, jump to data offset */
        fp = fopen(data_file,"r");
        if (fp == NULL) {
                fprintf(stderr, "\nError opening data file %s\n",data_file);
                exit(1);
        }

        ret = fseek(fp, (long)(hdr.vox_offset), SEEK_SET);
        if (ret != 0) {
                fprintf(stderr, "\nError doing fseek() to %ld in data file %s\n",(long)(hdr.vox_offset), data_file);
                exit(1);
        }


/********** allocate buffer and read first 3D volume from data file */
        data = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * hdr.dim[1]*hdr.dim[2]*hdr.dim[3]);
        if (data == NULL) {
                fprintf(stderr, "\nError allocating data buffer for %s\n",data_file);
                exit(1);
        }
        ret = fread(data, sizeof(MY_DATATYPE), hdr.dim[1]*hdr.dim[2]*hdr.dim[3], fp);
        if (ret != hdr.dim[1]*hdr.dim[2]*hdr.dim[3]) {
                fprintf(stderr, "\nError reading volume 1 from %s (%d)\n",data_file,ret);
                exit(1);
        }
        fclose(fp);


/********** scale the data buffer  */
        if (hdr.scl_slope != 0) {
                for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
                        data[i] = (data[i] * hdr.scl_slope) + hdr.scl_inter;
        }


/********** print mean of data */
        total = 0;
        for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
                total += data[i];
        total /= (hdr.dim[1]*hdr.dim[2]*hdr.dim[3]);
        fprintf(stderr, "\nMean of volume 1 in %s is %.3f\n",data_file,total);


        return(0);
}


/**********************************************************************
 *
 * write_nifti_file
 *
 * write a sample nifti1 (.nii) data file
 * datatype is float32
 * XYZT size is 64x64x16x10
 * XYZ voxel size is 1mm
 * TR is 1500ms
 *
 **********************************************************************/
static int write_nifti_file(char * const hdr_file, char * const data_file)
{
        nifti_1_header hdr;
        nifti1_extender pad={0,0,0,0};
        FILE *fp;
        int ret,i;
        MY_DATATYPE *data=NULL;
        short do_nii;


/********** make sure user specified .hdr/.img or .nii/.nii */
        if ( (strlen(hdr_file) < 4) || (strlen(data_file) < 4) ) {
                fprintf(stderr, "\nError: write files must end with .hdr/.img or .nii/.nii extension\n");
                exit(1);
        }

        if ( (!strncmp(hdr_file+(strlen(hdr_file)-4), ".hdr",4)) &&
             (!strncmp(data_file+(strlen(data_file)-4), ".img",4)) ) {
                do_nii = 0;
        }
        else if ( (!strncmp(hdr_file+(strlen(hdr_file)-4), ".nii",4)) &&
                  (!strncmp(data_file+(strlen(data_file)-4), ".nii",4)) ) {
                do_nii = 1;
        }
        else {
                fprintf(stderr, "\nError: file(s) to be written must end with .hdr/.img or .nii/.nii extension\n");
                exit(1);
        }


/********** fill in the minimal default header fields */
        bzero((void *)&hdr, sizeof(hdr));
        hdr.sizeof_hdr = MIN_HEADER_SIZE;
        hdr.dim[0] = 4;
        hdr.dim[1] = 64;
        hdr.dim[2] = 64;
        hdr.dim[3] = 16;
        hdr.dim[4] = 10;
        hdr.datatype = NIFTI_TYPE_FLOAT32;
        hdr.bitpix = 32;
        hdr.pixdim[1] = 1.0;
        hdr.pixdim[2] = 1.0;
        hdr.pixdim[3] = 1.0;
        hdr.pixdim[4] = 1.5;
        if (do_nii)
                hdr.vox_offset = (float) NII_HEADER_SIZE;
        else
                hdr.vox_offset = (float)0;
        hdr.scl_slope = 100.0;
        hdr.xyzt_units = NIFTI_UNITS_MM | NIFTI_UNITS_SEC;
        if (do_nii)
                strncpy(hdr.magic, "n+1\0", 4);
        else
                strncpy(hdr.magic, "ni1\0", 4);


/********** allocate buffer and fill with dummy data  */
        data = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*hdr.dim[4]);
        if (data == NULL) {
                fprintf(stderr, "\nError allocating data buffer for %s\n",data_file);
                exit(1);
        }

        for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*hdr.dim[4]; i++)
                data[i] = i / hdr.scl_slope;


/********** write first 348 bytes of header   */
        fp = fopen(hdr_file,"w");
        if (fp == NULL) {
                fprintf(stderr, "\nError opening header file %s for write\n",hdr_file);
                exit(1);
        }
        ret = fwrite(&hdr, MIN_HEADER_SIZE, 1, fp);
        if (ret != 1) {
                fprintf(stderr, "\nError writing header file %s\n",hdr_file);
                exit(1);
        }


/********** if nii, write extender pad and image data   */
        if (do_nii == 1) {

                ret = fwrite(&pad, 4, 1, fp);
                if (ret != 1) {
                        fprintf(stderr, "\nError writing header file extension pad %s\n",hdr_file);
                        exit(1);
                }

                ret = fwrite(data, (size_t)(hdr.bitpix/8), hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*hdr.dim[4], fp);
                if (ret != hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*hdr.dim[4]) {
                        fprintf(stderr, "\nError writing data to %s\n",hdr_file);
                        exit(1);
                }

                fclose(fp);
        }


/********** if hdr/img, close .hdr and write image data to .img */
        else {

                fclose(fp);     /* close .hdr file */

                fp = fopen(data_file,"w");
                if (fp == NULL) {
                        fprintf(stderr, "\nError opening data file %s for write\n",data_file);
                        exit(1);
                }
                ret = fwrite(data, (size_t)(hdr.bitpix/8), hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*hdr.dim[4], fp);
                if (ret != hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*hdr.dim[4]) {
                        fprintf(stderr, "\nError writing data to %s\n",data_file);
                        exit(1);
                }

                fclose(fp);
        }


        free(data);
        return(0);
}

/********** The main program for this test case */
int main(int argc, char * argv[])
{

char *hdr_file, *data_file;
short do_read=0;
short do_write=0;


/********** process commandline parameters */
if (argc != 4) {
        fprintf(stderr, "\nUsage: %s <-r|-w> <header file> <data file>\n",argv[0]);
        exit(1);
}

if (!strncmp(argv[1],"-r",2))
        do_read=1;
else if (!strncmp(argv[1],"-w",2))
        do_write=1;
else {
        fprintf(stderr, "\nUsage: %s <-r|-w> <header file> <data file>\n",argv[0]);
        exit(1);
}

hdr_file = argv[2];
data_file = argv[3];



/********** do the simple read or write */
if (do_read)
        read_nifti_file(hdr_file, data_file);
else if (do_write)
        write_nifti_file(hdr_file, data_file);


exit(0);
}