Dear FSL experts
I have looked into the cope, varcope, tstat, and sigmasquareds images in a feat stats folder.
Everything correspond to each other well except the sigmasquareds.nii.gz.
However, it should be the source of how each varcope are being generated.
I have looked into the source code (glimGls.cc, glimGls.h) which is listed at the end for reference, and the math calculating them seems sound and logical to me.
Form the code, one is easy to spot that the variable "scale" (=con*inv_xx*con.t()) is the one that transforms the sigmasqureds to each particular varcope. The "scale" is calculated in the for loop for contrasts.
Therefore, one is to except that when you use fslmaths to divide varcope by sigmasquareds, although different varcope will generate different scale value, but within a particular varcope/sigmasquareds, the resulting image should be consist of voxels with the same scalar value.
I have tested with several different feat folders I have, but each time it seems that in a particular varcope, each of the voxel has different "scale" factor.
Moreover, the resulting varcope/sigmasquareds image actually still resembles a brain.
Because the cope, varcope and tstat corresponding with each other well, and the code for their calculation also seemed correct.
I am wondering whether this is caused by how the sigmasqureds.nii.gz image was saved that cause this strange behavior?
(In the .h file, sigmaSquareds was declared as a RowVector, and cope, varcope were declared as a matrix.)
Or if varcope/sigmasquareds is not supposed to be the same "scale" number in every voxel in the first place, could anyone kindly explain to me why?
Thank you very much for your patience and kind assitance!
Best wishes,
Chen-Chia Lan
PhD student, Imperial College London, UK
Attending physician, Department of Psychiatry, Taichung Veterans General Hospital, Taiwan
Source code from glimGls.cc:
void GlimGls::setData(const ColumnVector& y, const Matrix& x, const int ind, const Matrix& tContrasts, const Matrix& fContrasts)
{
// compute b
Matrix inv_xx = pinv(x.t()*x);
b.Column(ind) = inv_xx*x.t()*y;
// compute r
Matrix R = I-x*inv_xx*x.t();
r = R*y;
sigmaSquareds(ind) = (r.t()*r).AsScalar()/R.Trace();
for ( int tContrast=1;tContrast <= tContrasts.Nrows();tContrast++) {
Matrix con=tContrasts.Row(tContrast);
copes(tContrast,ind)=(con*b.Column(ind)).AsScalar();
double scale((con*inv_xx*con.t()).AsScalar());
if ( scale <= 0 && con.MaximumAbsoluteValue() > 0 )
cerr << "Neff Error" << endl;
varcopes(tContrast,ind)=sigmaSquareds(ind)*scale;
}
for ( int fContrast=1;fContrast <= fContrasts.Nrows();fContrast++) {
Matrix fullFContrast( 0, tContrasts.Ncols() );
for (int tContrast=1; tContrast<=fContrasts.Ncols() ; tContrast++ )
if (fContrasts(fContrast,tContrast)==1) fullFContrast &= tContrasts.Row(tContrast);
fDof[fContrast-1]=fullFContrast.Nrows();
fstats(fContrast,ind) = (b.Column(ind).t()*fullFContrast.t()*(fullFContrast*inv_xx*fullFContrast.t()*sigmaSquareds(ind)).i()*fullFContrast*b.Column(ind)).AsScalar()/(float)fullFContrast.Nrows();
}
// set corrections SetCorrection(inv_xx, ind);
}
void GlimGls::Save(const volume<float>& mask, volume4D<float>& saveVolume, fslSurface<float, unsigned int>& saveSurface, const string& saveMode, const float reftdim)
{
// Need to save b, sigmaSquareds, corrections and dof
Log& logger = LogSingleton::getInstance();
for(int i = 1; i <= numParams; i++) //Beta
saveData(logger.getDir() + "/pe" + num2str(i),b.Row(i),saveVolume,mask,true,false,-1,true,NIFTI_INTENT_ESTIMATE,saveSurface,saveMode);
ColumnVector zstat(numTS);
for ( int tContrast=1;tContrast <= nContrasts;tContrast++) {
saveData(logger.getDir() + "/cope" + num2str(tContrast),copes.Row(tContrast),saveVolume,mask,true,false,-1,true,NIFTI_INTENT_ESTIMATE,saveSurface,saveMode);
saveData(logger.getDir() + "/varcope" + num2str(tContrast),varcopes.Row(tContrast),saveVolume,mask,true,false,-1,true,NIFTI_INTENT_ESTIMATE,saveSurface,saveMode);
RowVector tstat(SD(copes.Row(tContrast),sqrt(varcopes.Row(tContrast))));
T2z::ComputeZStats(varcopes.Row(tContrast).AsColumn(), copes.Row(tContrast).AsColumn(), dof, zstat);
saveData(logger.getDir() + "/tstat" + num2str(tContrast),tstat,saveVolume,mask,true,false,-1,true,NIFTI_INTENT_TTEST,saveSurface,saveMode);
saveData(logger.getDir() + "/zstat" + num2str(tContrast),zstat.AsRow(),saveVolume,mask,true,false,-1,true,NIFTI_INTENT_ZSCORE,saveSurface,saveMode);
}
for ( int fContrast=1;fContrast <= nFContrasts;fContrast++) {
saveData(logger.getDir() + "/fstat" + num2str(fContrast),fstats.Row(fContrast),saveVolume,mask,true,false,-1,true,NIFTI_INTENT_FTEST,saveSurface,saveMode);
F2z::ComputeFStats(fstats.Row(fContrast).AsColumn(), fDof[fContrast-1], dof, zstat);
saveData(logger.getDir() + "/zfstat" + num2str(fContrast),zstat.AsRow(),saveVolume,mask,true,false,-1,true,NIFTI_INTENT_ZSCORE,saveSurface,saveMode);
}
saveData(logger.getDir() + "/sigmasquareds",sigmaSquareds,saveVolume,mask,true,false,-1,false,-1,saveSurface,saveMode);
ColumnVector dofVec(1);
dofVec = dof;
write_ascii_matrix(logger.appendDir("dof"), dofVec);
}
|