roi_preprocess_05202018.cpp

/*********************************************************************
 *
 * This code is designed for preprocessing ROI timeseries and is based
 * on the open access data from humanconnectome.org. This code does
 * not perform motion correction, detrending, or any of the other 
 * preprocessing steps carried out by the HCP team.
 * 
 * 
 * 
 * 1. update_goodvoxels
 * 		The HCP data come with a file called "goodvoxels.nii" which
 * 		is a binary brain mask with additional voxels removed based		
 * 		on local coefficient of variation (inverse of SNR, std/mean).
 * 		This file is useful because it removes voxels with a low SNR,
 * 		subject to the surrounding voxels SNR. This provides similar  
 * 		voxel removal rates in areas of signal dropout as well as in
 * 		areas of normal BOLD signal. Since the olfactory regions are
 * 		located in areas of high signal dropout, this method is
 * 		particularly useful for our interests. HOWEVER, some voxels
 * 		are still included that have negative values for certain TRs.
 * 		Thus this function searches for remaining negative values and
 * 		provides an updated mask with a "0" indicating the voxels to
 * 		be excluded.
 * 
 *		Output: goodvoxels_new.csv
 *		Return: goodv_new, pointer to a binary int array (same as output)
 *
 * 2. check_snr
 * 		This function is included for checking the SNR and activity
 * 		of individual voxels. It prints out two files, one with the SNR 
 *      and one with the mean of each voxels timeseries. These two can 
 *		be used to calculate each voxels variance or standard deviation.
 *		It also returns the nuimber of TRs in the data file.
 * 
 * 		Output: mean_data.csv & snr_data.csv
 *		Return: hdr.dim[4], int of the 4th dimension (number of TRs)
 * 
 * 3. roi_extract
 * 		This is the most crucial function in this program. It reads in
 * 		the atlas and fmri data and uses this information to construct
 * 		the ROI timeseries. The end result will be a csv file with
 * 		each column corresponding to a single roi. (i.e. timepoints by
 * 		roi matrix).
 * 
 * 		Output: roi_data.csv
 *		Return: max, int maximum value from roi atlas. This should be the
 *				number of rois. If it isn't, renumber the atlas in MATLAB
 * 
 * 4. meancenter_col
 * 		Perfroms mean centering and matrix transposition. This function
 * 		relys upon the returned values from check_snr and roi_extract to
 * 		setup the data matrix (data[max][TRs]).
 * 
 * 5. voxel_summary
 * 		Prints out voxel removal information to a csv file. Examine in
 * 		Excel.
 * 
 * 
 * 
 * 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/
 * 
 * g++ roi_preprocess_0520201.cpp -o roi_preprocess_05202018
 * 
 * Run:
 * ./a.out rfmri.nii atlas.nii goodvoxels.nii
 *
 * Written by Thomas Campbell Arnold, FSU 2016
 *
 *********************************************************************/
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <fstream>
#include <iostream>
#include "nifti1.h"

typedef float MY_DATATYPE;

#define MIN_HEADER_SIZE 348
#define NII_HEADER_SIZE 352

using namespace std;

// PROTOTYPES
int * update_goodvoxels(char *data_file, char *goodv_file, char *basal_file);
int check_snr(char *data_file);
int roi_extract(char *data_file, char *atlas_file, int *goodv_update);
int meancenter_col(int max, int trs);
int voxel_summary(int max, char* data_file, char *atlas_file, char *goodv_file, int *goodv_update);

int main(int argc, char *argv[]) 
{
	
	char *data_file, *atlas_file, *goodv_file, *basal_file;
	
	// check for proper number of commandline arguements 
	if (argc != 5) {
		fprintf(stderr, "\nUsage: %s  <header file> <data file> <atlas.nii>\n",argv[0]);
        exit(1);
    }
	
	// get filenames from arguement
	data_file = argv[1];
	atlas_file = argv[2];
	goodv_file = argv[3];
	basal_file = argv[4];

	// set a time parameter
	time_t tstart, tend;
	tstart = time(0);
	
	// call functions
	int *goodv_update = update_goodvoxels(data_file, goodv_file, basal_file);
	int trs = check_snr(data_file);
	fprintf(stderr,"\nnumber of TRs = %d\n",trs);
	int max = roi_extract(data_file, atlas_file, goodv_update);
	fprintf(stderr,"\nNumber of Brain Regions = %d\n",max);
	meancenter_col(max,trs);
	fprintf(stderr,"\nfinished mean centering timeseries\n");
	voxel_summary(max, data_file, atlas_file, goodv_file, goodv_update);
	fprintf(stderr,"\nThe code has finished running. Check output.\n");

	// call back time parameter
	tend = time(0);
	fprintf(stderr, "\nIt took %.3f second(s) to run this code.\n",difftime(tend, tstart));
	
	exit(0);
}

/***************************** FIRST FUNCTION *****************************
*
*
*
*/

int * update_goodvoxels(char *data_file, char *goodv_file, char *basal_file)
{
	nifti_1_header hdr;
	nifti_1_header goodv_hdr;
	nifti_1_header basal_hdr;
	FILE *f;
	int ret, i, t, neg_data, neg_goodv;
	MY_DATATYPE *data=NULL, *goodv_data=NULL, *basal_data=NULL;
	
	// initialize negative counts at 0
	neg_data = 0;
	neg_goodv = 0;
	
	// FIRST LOAD IN GOODVOXELS.NII INFORMATION
	
	// open hdr file
	f = fopen(goodv_file,"r");
	if (f == NULL) {
        fprintf(stderr, "\nError opening header file %s\n",goodv_file);
        exit(1);
    }

    // read hdr information
    ret = fread(&goodv_hdr, MIN_HEADER_SIZE, 1, f);
    if (ret != 1) {
        fprintf(stderr, "\nError reading header file %s\n",goodv_file);
        exit(1);
    }
    fclose(f);
    
    // print hdr information to screen
    fprintf(stderr, "\n%s header information:",goodv_file);
    fprintf(stderr, "\nXYZT dimensions: %d %d %d %d",goodv_hdr.dim[1],goodv_hdr.dim[2],goodv_hdr.dim[3],goodv_hdr.dim[4]);
    fprintf(stderr, "\nDatatype code and bits/pixel: %d %d",goodv_hdr.datatype,goodv_hdr.bitpix);
    fprintf(stderr, "\nScaling slope and intercept: %.6f %.6f",goodv_hdr.scl_slope,goodv_hdr.scl_inter);
    fprintf(stderr, "\nByte offset to data in datafile: %ld",(long)(goodv_hdr.vox_offset));
    fprintf(stderr, "\n");
    
    // open data file jump to vox_offset before beginning to read
    f = fopen(goodv_file,"r");
    if (f == NULL) {
        fprintf(stderr, "\nError opening data file %s\n",goodv_file);
        exit(1);
	}
	ret = fseek(f, (long)(goodv_hdr.vox_offset), SEEK_SET);
	
	// allocate memory for holding mean
	//float goodv_data[goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]];
	//memset(goodv_data,0,goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]*sizeof(float));
	
	// allocate space to hold the data
	goodv_data = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]);
	if (goodv_data == NULL) {
		fprintf(stderr, "\nError allocating data buffer for %s\n",data_file);
		exit(1);
	}
	
	ret = fread(goodv_data, sizeof(MY_DATATYPE), goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3], f);
	if (ret != goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]) {
		fprintf(stderr, "\nError reading volume 1 from %s (%d)\n",goodv_data,ret);
		exit(1);
	}
	
	// close data file
	fclose(f);
	
	// LOAD IN BASAL GANGLION MASK
	
	// open hdr file
	f = fopen(basal_file,"r");
	if (f == NULL) {
        fprintf(stderr, "\nError opening header file %s\n",basal_file);
        exit(1);
    }

    // read hdr information
    ret = fread(&basal_hdr, MIN_HEADER_SIZE, 1, f);
    if (ret != 1) {
        fprintf(stderr, "\nError reading header file %s\n",basal_file);
        exit(1);
    }
    fclose(f);
    
    // print hdr information to screen
    fprintf(stderr, "\n%s header information:",basal_file);
    fprintf(stderr, "\nXYZT dimensions: %d %d %d %d",basal_hdr.dim[1],basal_hdr.dim[2],basal_hdr.dim[3],basal_hdr.dim[4]);
    fprintf(stderr, "\nDatatype code and bits/pixel: %d %d",basal_hdr.datatype,basal_hdr.bitpix);
    fprintf(stderr, "\nScaling slope and intercept: %.6f %.6f",basal_hdr.scl_slope,basal_hdr.scl_inter);
    fprintf(stderr, "\nByte offset to data in datafile: %ld",(long)(basal_hdr.vox_offset));
    fprintf(stderr, "\n");
    
    // open data file jump to vox_offset before beginning to read
    f = fopen(basal_file,"r");
    if (f == NULL) {
        fprintf(stderr, "\nError opening data file %s\n",basal_file);
        exit(1);
	}
	ret = fseek(f, (long)(basal_hdr.vox_offset), SEEK_SET);
	
	// allocate memory for holding mean
	//float goodv_data[goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]];
	//memset(goodv_data,0,goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]*sizeof(float));
	
	// allocate space to hold the data
	basal_data = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * basal_hdr.dim[1]*basal_hdr.dim[2]*basal_hdr.dim[3]);
	if (basal_data == NULL) {
		fprintf(stderr, "\nError allocating data buffer for %s\n",basal_file);
		exit(1);
	}
	
	ret = fread(basal_data, sizeof(MY_DATATYPE), basal_hdr.dim[1]*basal_hdr.dim[2]*basal_hdr.dim[3], f);
	if (ret != basal_hdr.dim[1]*basal_hdr.dim[2]*basal_hdr.dim[3]) {
		fprintf(stderr, "\nError reading volume 1 from %s (%d)\n",basal_data,ret);
		exit(1);
	}
	
	// close data file
	fclose(f);
	
	// BEGIN LOADING RFMRI DATA
	
	// open hdr file
	f = fopen(data_file,"r");
	if (f == NULL) {
        fprintf(stderr, "\nError opening header file %s\n",data_file);
        exit(1);
    }

    // read hdr information
    ret = fread(&hdr, MIN_HEADER_SIZE, 1, f);
    if (ret != 1) {
        fprintf(stderr, "\nError reading header file %s\n",data_file);
        exit(1);
    }
    fclose(f);
    
    // print hdr information to screen
    fprintf(stderr, "\n%s header information:",data_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 data file jump to vox_offset before beginning to read
    f = fopen(data_file,"r");
    if (f == NULL) {
        fprintf(stderr, "\nError opening data file %s\n",data_file);
        exit(1);
	}
	ret = fseek(f, (long)(hdr.vox_offset), SEEK_SET);
	
	// allocate memory for holding mean
	//float data[hdr.dim[1]*hdr.dim[2]*hdr.dim[3]];
	//memset(data,0,hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*sizeof(float));
	
	// loop through each of the TRs (1:hdr.dim[4])
	for(t=1; t<=hdr.dim[4]; t++){
		
		// allocate space to hold the data
		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], f);
		if (ret != hdr.dim[1]*hdr.dim[2]*hdr.dim[3]) {
			fprintf(stderr, "\nError reading volume %d from %s (%d)\n",t,data_file,ret);
			exit(1);
		}
		
		// add data from new TR to the sum
		for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
		{
			// if data value is below zero add to count
			if (data[i] < 0){
				neg_data++;
			}
			
			
			// change goodv_data if it is a basal voxel
			if (goodv_data[i] == 0){
				//REMOVED BASAL AFFECT if (basal_data[i] == 1){
					goodv_data[i] = 0;
				//}
			}
			
			// multiply data by goodvoxels
			data[i] = data[i] * goodv_data[i];
			
			/* REMOVED: removed any voxel that had a value below zero (not valid after DPARSFA)
			// if data value is below zero add to count (after goodvoxels)
			if (data[i] < 0){
				neg_goodv++;
				goodv_data[i] = 0;
			}*/
		}
		
	}
	
	// close data file
	fclose(f);

	// print results
	fprintf(stderr, "\nTotal instances of negative TRs in raw data = %d\n",neg_data);
	fprintf(stderr, "\nInstances of negative TRs after goodvoxels correction = %d\n",neg_goodv);
	
	// open file to print out new goodvoxels if needed
	if (neg_goodv > 0){
		f = fopen("goodvoxels_new.csv","w");
		for (i=0; i<goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]; i++){
				fprintf(f,"%d,",goodv_data[i]);
		}
		fclose(f);
	}
	
	// convert to int
	static int goodv_new[1000000];
	// static does not work becuase it is defined by a variable, this should be changed.
	//static int goodv_new[goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]]; 
	memset(goodv_new,0,goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]*sizeof(int));
	for (i=0; i<goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]; i++)
	{
		goodv_new[i] = goodv_data[i];
	}
	

	delete[] data;
	delete[] goodv_data;
	return goodv_new;
};
 
/***************************** SECOND FUNCTION *****************************
*
*
*
*/
 
int check_snr(char *data_file)
{
	nifti_1_header hdr;
	FILE *f, *f2;
	int ret, i, t;
	MY_DATATYPE *data=NULL, *mean_data=NULL, *var_data=NULL;
	
	// open hdr file
	f = fopen(data_file,"r");
	if (f == NULL) {
        fprintf(stderr, "\nError opening header file %s\n",data_file);
        exit(1);
    }

    // read hdr information
    ret = fread(&hdr, MIN_HEADER_SIZE, 1, f);
    if (ret != 1) {
        fprintf(stderr, "\nError reading header file %s\n",data_file);
        exit(1);
    }
    fclose(f);
    
    // open data file jump to vox_offset before beginning to read
    f = fopen(data_file,"r");
    if (f == NULL) {
        fprintf(stderr, "\nError opening data file %s\n",data_file);
        exit(1);
	}
	ret = fseek(f, (long)(hdr.vox_offset), SEEK_SET);
	
	// allocate memory for holding mean
	mean_data = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * hdr.dim[1]*hdr.dim[2]*hdr.dim[3]);
	//float mean_data[hdr.dim[1]*hdr.dim[2]*hdr.dim[3]];
	memset(mean_data,0,hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*sizeof(MY_DATATYPE));

	// loop through each of the TRs (1:hdr.dim[4])
	for(t=1; t<=hdr.dim[4]; t++){


		// allocate space to hold the data
		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], f);
		if (ret != hdr.dim[1]*hdr.dim[2]*hdr.dim[3]) {
			fprintf(stderr, "\nError reading volume %d from %s (%d)\n",t,data_file,ret);
			exit(1);
		}
		
		// add data from new TR to the sum
		for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
		{
			mean_data[i] += data[i];
		}
		
	}

	// divide the sum of all TR intensities by total number of TTs
	for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
	{
		mean_data[i] /= hdr.dim[4];
	}
	

	// close data file
	fclose(f);
	
	// BEGIN CALCULATION OF VARIANCE
	
	// open data file jump to vox_offset before beginning to read
    f = fopen(data_file,"r");
    if (f == NULL) {
        fprintf(stderr, "\nError opening data file %s\n",data_file);
        exit(1);
	}
	ret = fseek(f, (long)(hdr.vox_offset), SEEK_SET);

	// allocate memory for holding mean
	//float var_data[hdr.dim[1]*hdr.dim[2]*hdr.dim[3]];
	//memset(var_data,0,hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*sizeof(float));
	var_data = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * hdr.dim[1]*hdr.dim[2]*hdr.dim[3]);
	memset(var_data,0,hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*sizeof(MY_DATATYPE));

	// loop through each of the TRs (1:hdr.dim[4])
	for(t=1; t<=hdr.dim[4]; t++){
		
		// allocate space to hold the data
		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], f);
		if (ret != hdr.dim[1]*hdr.dim[2]*hdr.dim[3]) {
			fprintf(stderr, "\nError reading volume %d from %s (%d)\n",t,data_file,ret);
			exit(1);
		}
		
		// add data from new TR to the sum
		for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
		{
			var_data[i] += (data[i] - mean_data[i])*(data[i] - mean_data[i]);
		}
		
	}
	
	// divide the sum of all TR intensities by total number of TTs
	for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
	{
		var_data[i] /= hdr.dim[4];
	}
	
	// close data file
	fclose(f);

	// variables to count voxels w/ mean > 1 and calculate mean SNR
	int n = 0;
	float tot = 0;
	
	// open file for writing
	f = fopen("mean_data.csv", "w");
	f2 = fopen("snr_data.csv", "w");
	
	// calculate SNR: divide mean by sqrt of variance (std)
	for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
	{
		if (mean_data[i] > 1){
			fprintf(f,"%f,",mean_data[i]); // print mean before calculating snr
			mean_data[i] /= sqrt(var_data[i]);
			fprintf(f2,"%f,",mean_data[i]); // print snr (note: mean_data is altered)
			n++;
			tot += mean_data[i];
		}else{
			mean_data[i] = 0;
			fprintf(f,"%f,",mean_data[i]); // NOTE: This does not print actual mean or SNR, it prints 0.
			fprintf(f2,"%f,",mean_data[i]); //      Mean/SNR could be in below 1 or even negative.
		}
		
	}
	
	// close output files
	fclose(f);
	fclose(f2);
	
	// print out voxel count and mean SNR to screen
	tot /= n;
	fprintf(stderr, "\nTotal voxels before goodvoxels = %d\n\nmean SNR = %f\n",n,tot);
	

	delete[] data;
	delete[] mean_data;
	delete[] var_data;
	return hdr.dim[4];
};
 
/***************************** THIRD FUNCTION *****************************
*
*
*
*/
 
int roi_extract(char *data_file, char *atlas_file, int *goodv_update){
	nifti_1_header hdr;
	nifti_1_header aal_hdr;
	FILE *fp,*f,*aal;
	int ret,i,t,max;
	MY_DATATYPE *data=NULL, *mean_intensity=NULL;
	double *roi_data=NULL;
	int *roi_array=NULL;
	
	/**********************************************************************
	 *
	 * Process Atlas
	 *
	**********************************************************************/
	
	/********** open and read header */
	aal = fopen(atlas_file,"r");
	if (aal == NULL) {
	        fprintf(stderr, "\nError opening header file %s\n",atlas_file);
	        exit(1);
	}
	
	ret = fread(&aal_hdr, MIN_HEADER_SIZE, 1, aal);
	if (ret != 1) {
	        fprintf(stderr, "\nError reading header file %s\n",atlas_file);
	        exit(1);
	}
	fclose(aal);
	
	/********** print a little header information */
	fprintf(stderr, "\n%s header information:",atlas_file);
	fprintf(stderr, "\nXYZT dimensions: %d %d %d %d",aal_hdr.dim[1],aal_hdr.dim[2],aal_hdr.dim[3],aal_hdr.dim[4]);
	fprintf(stderr, "\nDatatype code and bits/pixel: %d %d",aal_hdr.datatype,aal_hdr.bitpix);
	fprintf(stderr, "\nScaling slope and intercept: %.6f %.6f",aal_hdr.scl_slope,aal_hdr.scl_inter);
	fprintf(stderr, "\nByte offset to data in datafile: %ld",(long)(aal_hdr.vox_offset));
	fprintf(stderr, "\n");
	
	/*** check if there are multiple timepoints in atlas ***/
	if (aal_hdr.dim[4] >= 2){
			fprintf(stderr, "\nError - multiple timepoints in atlas file: %s\n",atlas_file);
	        exit(1);
	}
	
	/********** open the datafile, jump to data offset */
	aal = fopen(atlas_file,"r");
	if (aal == NULL) {
	        fprintf(stderr, "\nError opening data file %s\n",atlas_file);
	        exit(1);
	}
	ret = fseek(aal, (long)(aal_hdr.vox_offset), SEEK_SET);
	
	/********** allocate buffer and read first 3D volume from data file */
	roi_data = (double *) malloc(sizeof(double) * aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]);
	if (roi_data == NULL) {
	        fprintf(stderr, "\nError allocating data buffer for %s\n",atlas_file);
	        exit(1);
	}
	ret = fread(roi_data, sizeof(double), aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3], aal);
	if (ret != aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]) {
	        fprintf(stderr, "\nError reading volume 1 from %s (%d)\n",atlas_file,ret);
	        exit(1);
	}
	
	/********** scale the data buffer  
	if (aal_hdr.scl_slope != 0) {
	        for (i=0; i<aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]; i++)
	                roi_data[i] = (roi_data[i] * aal_hdr.scl_slope) + aal_hdr.scl_inter;
	}*/
	
	/************ get max value of roi****************/
	max = roi_data[0]; // declare max is first value in data array
	for (i=0; i<aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]; i++)
	{
	 if(max < roi_data[i]){
	 	max = roi_data[i];
	 }
	}
	
	int n = 0;
	/*** convert from float to int ***/
	// multiply by goodvoxels to remove voxels from atlas
	//int roi_array[aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]];
	//memset(roi_array,0,aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]*sizeof(int));
	roi_array = (int *) malloc(sizeof(int) * aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]);
	memset(roi_array,0,aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]*sizeof(int));


	for (i=0; i<aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]; i++)
	{
		roi_array[i] = roi_data[i] * goodv_update[i];
		if (goodv_update[i] == 1){
			n++;
		}
	}
	
	fprintf(stderr, "\nTotal voxels after goodvoxels update = %d\n",n);
	
	fclose(aal);
	
	/**********************************************************************
	 *
	 * Process rfMRI data
	 *
	**********************************************************************/
	
	/********** open and read header */
	fp = fopen(data_file,"r");
	if (fp == NULL) {
	        fprintf(stderr, "\nError opening header file %s\n",data_file);
	        exit(1);
	}
	
	ret = fread(&hdr, MIN_HEADER_SIZE, 1, fp);
	if (ret != 1) {
	        fprintf(stderr, "\nError reading header file %s\n",data_file);
	        exit(1);
	}
	fclose(fp);
	
	//float mean_intensity[hdr.dim[1]*hdr.dim[2]*hdr.dim[3]];
	//memset(mean_intensity,0,hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*sizeof(float));
	mean_intensity = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * hdr.dim[1]*hdr.dim[2]*hdr.dim[3]);
	memset(mean_intensity,0,hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*sizeof(MY_DATATYPE));

	f = fopen("roi_data.csv", "w"); /**open CSV file**/
	fp = fopen(data_file,"r");
	for(t=1; t<=hdr.dim[4]; t++)
	{
	/********** open the datafile, jump to data offset */
	if (fp == NULL) {
	        fprintf(stderr, "\nError opening data file %s\n",data_file);
	        exit(1);
	}
	if (t == 1){
		ret = fseek(fp, NII_HEADER_SIZE, SEEK_SET);
	}
	
	/********** 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);
	}
	
	/********** 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;
	}*/
	
	/*float fmri_array[hdr.dim[1]*hdr.dim[2]*hdr.dim[3]];
	memset(fmri_array,0,hdr.dim[1]*hdr.dim[2]*hdr.dim[3]*sizeof(float));
	for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++)
	{
		fmri_array[i] = data[i];
	}*/

	 /***** roi count *****/
	 int roi_count[max+1];
	 memset(roi_count,0,(max+1)*sizeof(int));


	/*** print mean roi timeseries of data ***/
	float total[max+1];
	memset(total,0,(max+1)*sizeof(float));
	for (i=0; i<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; i++){
		if(data[i] >= 1){
	        total[roi_array[i]] = total[roi_array[i]] + data[i];
	        roi_count[roi_array[i]]++;
	    }
	}


	
	/*** divide by number of voxels included in roi ***/
	for (i=1; i<max+1; i++){
	        total[i] /= roi_count[i];
	}
	
	
	/*************** write CSV file **********
	added by TCA 11-20-2015 *****************/
	
	for (i=1; i<max+1; i++)
	{
		if (i<max){
			fprintf(f,"%f,",total[i]);
		}else{
			fprintf(f,"%f",total[i]);
		}
	}
	
	fprintf(f,"\n");
	
	}
	
	 fclose(fp);
	 fclose(f);
	 
	delete[] data;
	delete[] mean_intensity;
	delete[] roi_data;
	delete[] roi_array;
	return max;

};
 
/***************************** FOURTH FUNCTION *****************************
*
*
*
*/
 
int meancenter_col(int max, int trs){

	FILE *f,*fp;
	char delim;
	int i, j, n, val;
	float temp;

	// allocate memory for data
	float *data[max];
	for (i=0; i<max; i++){
	data[i] = (float *)malloc(trs * sizeof(float));
	}

	// open csv file
	f = fopen("roi_data.csv","r");

	// used to index matrix, start at zero
	i=0;
	j=0;
	n=0;
	
	/***************************************************************************
	 *
	 * The data must be read in from a .csv file and then transposed. The 
	 * output from nifti_simu.c is in 1200x605 (the last column is zeros). This
	 * needs to be transposed into a 605(rois) x 1200(timepoints) matrix. The
	 * result of this will be called "data"
	 *
	***************************************************************************/
	//fscanf(f,"%i",&temp);
	while(val = fscanf(f,"%f,",&temp)){
		if(i <= max-1){
			// is int value, put in array
			data[i][j] = temp;
			//fprintf(stderr,"\ni=%i j=%i\n",i,j); // used to check proper i,j output
			if(i>=max-1){
				i=0;
				j++;
				if(j>=trs){break;}
			}else{
				i++;
			}
		}else{
			break;
		}
	}

	/**** remove the mean ******/
	float tp_sum;
	for(i=0; i<max; i++){
		tp_sum = 0;
		for(j=0; j<trs; j++){
			tp_sum += data[i][j];
		}
		tp_sum /= trs;
		for(j=0; j<trs; j++){
			data[i][j] -= tp_sum;
		}
	}
	
	/***************************************************************************
	 *
	 * This section can be used to verify output file is properly transposed
	 *
	 *
	 ***************************************************************************/
	// open csv file to write
	fp = fopen("mc.csv","w");			
	for(j=0; j<trs; j++){
		for(i=0; i<max; i++){
			if (i<max-1){
				fprintf(fp,"%f,",data[i][j]);
			}else{
				fprintf(fp,"%f",data[i][j]);
			}
			
		}
		fprintf(fp,"\n",data[i][j]);
	}
	
	
	// close csv file to write
	fclose(f);
	fclose(fp);
	
    return 0;
};


int voxel_summary(int max, char* data_file, char *atlas_file, char *goodv_file, int *goodv_update){
	
	nifti_1_header goodv_hdr;
	nifti_1_header aal_hdr;
	nifti_1_header hdr;
	FILE *f,*aal;
	int ret,i,j;
	MY_DATATYPE *data=NULL, *goodv_data=NULL;
	double *roi_data=NULL;
	int *goodv_array=NULL, *roi_array=NULL;
	
	// load in atlas information
	/********** open and read header */
	aal = fopen(atlas_file,"r");
	if (aal == NULL) {
	        fprintf(stderr, "\nError opening header file %s\n",atlas_file);
	        exit(1);
	}
	
	ret = fread(&aal_hdr, MIN_HEADER_SIZE, 1, aal);
	if (ret != 1) {
	        fprintf(stderr, "\nError reading header file %s\n",atlas_file);
	        exit(1);
	}
	fclose(aal);
	
	/********** open the datafile, jump to data offset */
	aal = fopen(atlas_file,"r");
	if (aal == NULL) {
	        fprintf(stderr, "\nError opening data file %s\n",atlas_file);
	        exit(1);
	}
	ret = fseek(aal, (long)(aal_hdr.vox_offset), SEEK_SET);
	
	/********** allocate buffer and read first 3D volume from data file */
	roi_data = (double *) malloc(sizeof(double) * aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]);
	if (roi_data == NULL) {
	        fprintf(stderr, "\nError allocating data buffer for %s\n",atlas_file);
	        exit(1);
	}
	ret = fread(roi_data, sizeof(double), aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3], aal);
	if (ret != aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]) {
	        fprintf(stderr, "\nError reading volume 1 from %s (%d)\n",atlas_file,ret);
	        exit(1);
	}
	
	// close data file
	fclose(aal);
	
	/*** convert from float to int ***/
	//int roi_array[aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]];
	//memset(roi_array,0,aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]*sizeof(int));
	roi_array = (int *) malloc(sizeof(int) * aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]);
	memset(roi_array,0,aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]*sizeof(int));

	for (i=0; i<aal_hdr.dim[1]*aal_hdr.dim[2]*aal_hdr.dim[3]; i++)
	{
		roi_array[i] = roi_data[i];
	}
	
	// load in goodvoxels.nii information
	
	// open hdr file
	f = fopen(goodv_file,"r");
	if (f == NULL) {
        fprintf(stderr, "\nError opening header file %s\n",goodv_file);
        exit(1);
    }

    // read hdr information
    ret = fread(&goodv_hdr, MIN_HEADER_SIZE, 1, f);
    if (ret != 1) {
        fprintf(stderr, "\nError reading header file %s\n",goodv_file);
        exit(1);
    }
    fclose(f);
    
    // open data file jump to vox_offset before beginning to read
    f = fopen(goodv_file,"r");
    if (f == NULL) {
        fprintf(stderr, "\nError opening data file %s\n",goodv_file);
        exit(1);
	}
	ret = fseek(f, (long)(goodv_hdr.vox_offset), SEEK_SET);
	
	// allocate space to hold the data
	goodv_data = (MY_DATATYPE *) malloc(sizeof(MY_DATATYPE) * goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]);
	if (goodv_data == NULL) {
		fprintf(stderr, "\nError allocating data buffer for %s\n",data_file);
		exit(1);
	}
	
	ret = fread(goodv_data, sizeof(MY_DATATYPE), goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3], f);
	if (ret != goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]) {
		fprintf(stderr, "\nError reading volume 1 from %s (%d)\n",goodv_data,ret);
		exit(1);
	}
	
	// close data file
	fclose(f);
	
	/*** convert from float to int ***/
	//int goodv_array[goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]];
	//memset(goodv_array,0,goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]*sizeof(int));
	goodv_array = (int *) malloc(sizeof(int) * goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]);
	memset(goodv_array,0,goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]*sizeof(int));

	for (i=0; i<goodv_hdr.dim[1]*goodv_hdr.dim[2]*goodv_hdr.dim[3]; i++)
	{
		goodv_array[i] = goodv_data[i];
	}
	
	// load in data information, just first slice

	// open hdr file
	f = fopen(data_file,"r");
	if (f == NULL) {
        fprintf(stderr, "\nError opening header file %s\n",data_file);
        exit(1);
    }

    // read hdr information
    ret = fread(&hdr, MIN_HEADER_SIZE, 1, f);
    if (ret != 1) {
        fprintf(stderr, "\nError reading header file %s\n",data_file);
        exit(1);
    }
    fclose(f);
    
    // open data file jump to vox_offset before beginning to read
    f = fopen(data_file,"r");
    if (f == NULL) {
        fprintf(stderr, "\nError opening data file %s\n",data_file);
        exit(1);
	}
	ret = fseek(f, (long)(hdr.vox_offset), SEEK_SET);
	
	// allocate space to hold the data
	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], f);
	if (ret != hdr.dim[1]*hdr.dim[2]*hdr.dim[3]) {
		fprintf(stderr, "\nError reading volume 1 from %s (%d)\n",data,ret);
		exit(1);
	}
	
	// close data file
	fclose(f);
	
	// open file for writing
	f = fopen("voxels_summary.csv","w");
	
	// count voxels for each condition
	int roi_count[max+1];
	memset(roi_count,0,(max+1)*sizeof(int));
	int goodv_count[max+1];
	memset(goodv_count,0,(max+1)*sizeof(int));
	int goodv_new_count[max+1];
	memset(goodv_new_count,0,(max+1)*sizeof(int));
	
	for (j=0; j<hdr.dim[1]*hdr.dim[2]*hdr.dim[3]; j++)
	{
		
		if(data[j] >= 1)
		{
			roi_count[roi_array[j]]++;
			goodv_count[roi_array[j]*goodv_array[j]]++;
			goodv_new_count[roi_array[j]*goodv_update[j]]++;
		};
		
	};
	
	for (i=0; i<max+1; i++)
	{
		fprintf(f,"%i,%i,%i,%i,%f,%f\n",i,roi_count[i],goodv_count[i],goodv_new_count[i],((goodv_count[i]*1.0)/(roi_count[i]*1.0))*100,((goodv_new_count[i]*1.0)/(roi_count[i]*1.0))*100);
	}
	
	// close output file
	fclose(f);
	
	return 0;
};

/*********************************************************************
 *
 * Helpful comments from parent package. Note the original code was designed
 * to be compiled in c, not c++. The links near the bottom can supply you
 * with the original code, which is brief. 
 *
 * 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.
 *
 *********************************************************************/