Revision 1d34d9c07e34cf0306445d49f7206994f5a79537 authored by James Scott on 19 July 2016, 08:40:22 UTC, committed by James Scott on 19 July 2016, 08:40:22 UTC
1 parent a7ab0a6
fdrBH95.cpp
#include<cstdio>
#include<cstdlib>
#include <algorithm>
#include<iostream>
#include<fstream>
#include<vector>
#include<string>
#include<cstring>
#include<sstream>//for convert int to string
#include<mex.h>
using namespace std;
void fdr_pvalue(double &q,int pN,vector<double> &p_value,vector<int> &signal_pvalue,int &R,double &threshold);
void mexFunction(int nlhs, mxArray *plhs[], int nrhs,const mxArray *prhs[])
{
int i;
if(nrhs!=4)
{
mexErrMsgTxt("Four inputs are required.");
}
//The first input
int Get_N1= mxGetN(prhs[0]);
if(mxGetM(prhs[0])!=1 || Get_N1!=4)
mexErrMsgTxt("Input argument 1 must be a row vector of size 4.");
double *input1=mxGetPr(prhs[0]);
for(i=0;i<3;i++)
{
if((input1[i]-(double((int)input1[i])))>0 || input1[i]<1)
mexErrMsgTxt("The first three elements of input argument 1 must be positive integers.");
}
if(!(input1[3]<1 && input1[3]>0))
mexErrMsgTxt("The last element of input argument 1 must be greater than 0 and less than 1.");
int LX=(int)input1[0];//dimension of the lattice on x-axis
int LY=(int)input1[1];//dimension of the lattice on y-axis
int LZ=(int)input1[2];//dimension of the lattice on z-axis
double q=input1[3];// control fdr at level q
int Num=LX*LY*LZ;//total voxels of the lattice on the image grid
//////////////////////////////////
//The 2nd to 4th inputs
char *file[3];// the directory of file,
/*file[0]: location file for the brain region of interest (ROI) on the LX*LY*LZ image grid
*file[1]: observed data (p-value) on the whole image grid including the brain part and the non-brain part
*file[2]: output file of signals for the ROI
*/
int filelen[3];//the length of directory
int status[3];
for(i=1;i<nrhs;i++)
{
// Input must be a string.
if (mxIsChar(prhs[i]) != 1)
mexErrMsgTxt("Input arguments from #3 must be a string.");
// Input must be a row vector.
if (mxGetM(prhs[i]) != 1)
mexErrMsgTxt("Input arguments from #3 must be a row vector.");
// Get the length of the input string.
int i_1=i-1;
filelen[i_1] = (mxGetM(prhs[i]) * mxGetN(prhs[i])) + 1;
/* Allocate memory for input and output strings. */
file[i_1] =(char *)mxCalloc(filelen[i_1], sizeof(char));
//Copy the string data from prhs[i] into a C++ string *file
status[i_1] = mxGetString(prhs[i], file[i_1], filelen[i_1]);
if (status[i_1] != 0)
mexWarnMsgTxt("Not enough space. String is truncated.");
}
vector<int> Loca;//Location of the ROI
ifstream ifs_loca(file[0]);
double temp_loca;
while(ifs_loca>>temp_loca)
Loca.push_back(int(temp_loca)-1);// The entry location of C++ array starts from 0, but vectorized image lattice location starts from 1
int PN=Loca.size();// the number of voxels in the ROI
vector<double> pvalue_all;//observed data (p-value) on the whole image grid including the brain part and the non-brain part
ifstream ifs_pvalue(file[1]);//
double temp_pvalue;
while(ifs_pvalue>>temp_pvalue)
pvalue_all.push_back(temp_pvalue);// read in y_ori
vector<double> pvalue(PN);//the p-values of the ROI
for(int pi=0;pi<PN;pi++)//read in just the p-value of the ROI
{
i=Loca[pi];
pvalue[pi]=pvalue_all[i];
}
vector<int> signals(PN);//The signals of the ROI,1=signal, 0=non-signal
double threshold;//threshold in the BH procedure
int R;//the number of signals (rejections)
fdr_pvalue(q,PN,pvalue,signals,R,threshold);
//output the signal file
ofstream signal_file;
signal_file.open(file[2]);
for(i=0;i<PN;i++)
{
signal_file<<signals[i]<<" ";
}
signal_file.close();
plhs[0]=mxCreateDoubleMatrix(1,3,mxREAL);
double *out;
out=mxGetPr(plhs[0]);
out[0]=R;
out[1]=PN;
out[2]=threshold;
}
void fdr_pvalue(double &q,int pN,vector<double> &p_value,vector<int> &signal_pvalue,int &R,double &threshold)
{
vector<double> sort_pvalue(p_value);
sort(sort_pvalue.begin(),sort_pvalue.end());
for(R=pN; R>0 && sort_pvalue[R-1]*pN>R*q; R--);
if(R!=0)
{
threshold=sort_pvalue[R-1];
int time=0;// for the tie pvalue problem
for(int j=0; j<pN && time<R; j++)
{
if(p_value[j]<=threshold)
{
signal_pvalue[j]=1;
time++;
}
}
}
else
{
R=0;
threshold=q/pN;
}
}
Computing file changes ...