Sorry I didn't notice your comment about reverse engineering cifti, it is intended to be an open format, we do have documentation here:
http://www.nitrc.org/projects/cifti/ In particular, http://www.nitrc.org/plugins/mwiki/index.php/cifti:ConnectivityMatrixFileFormatsdetails how things work in cifti-1. It is intended to describe the format fully, but in case it doesn't, questions are welcome. Yes, in cifti-1, the nifti dimensions were stored backwards in the header, we are close to releasing a proposal for cifti-2 that corrects this and other design quirks, and adds additional features. As for the surface data that is not in the cifti file, we generally exclude the medial wall from cifti, and when changing it back to a metric (gifti functional) file, we put zeros where we don't have data (and optionally output the ROI of where we have data), because gifti doesn't have a provision for "no data here". The list of vertices in that "cifti[0][2][0][0].text" element (NodeIndices) you found says which vertices are included, and their order in the cifti file. There is also a special case (which we plan to remove in cifti-2) where if the list is empty and the IndexCount equals SurfaceNumberOfNodes, it means all vertices, in ascending order. However, this case should be fairly rare, since we exclude the medial wall. Tim On Mon, Dec 16, 2013 at 7:39 AM, basile pinsard < [email protected]> wrote: > Hi again, > > it really seems to be a matter of C/Fortran order storage issue, > reopening the timeseries data in C order seems to give me the right > data. > > In [120]: data = np.memmap(nii.get_filename(), > nii.get_data_dtype(), > shape=nii.shape[4:], > offset=nii.header.get_data_offset(), > order='C', > mode='r') > > In [121]: data.shape > Out[121]: (91282, 1200) > > In [122]: np.allclose(data[:29696],left_tss[:,idxs].T) > Out[122]: True > > > Cheers. > > Basile > > > On Sun, Dec 15, 2013 at 10:15 PM, basile pinsard > <[email protected]> wrote: > > Hi Matt, > > > > sorry to come back on this again but I am mainly trying to also > > understand the CIFTI format, as I am also interested by formats to be > > not only readable by closed source (even if freely distributed) binary > > but could be read by others software which is reason why I spent hours > > trying to reverse-engineer it (and finding bugs it seems but the > > reason you mentionned is not the one causing problems in my analysis). > > > > Now that I have each hemisphere timeseries stored in gifti separately > > I was able to compare to what is stored in the cifti. > > Loading the cifti in python with nibabel gave me an array of size > > (1,1,1,1,91282,1200) and from the gifti file I got the data from left > > cortex as time frame in darrays. Here is what I have done in ipython > > to understand why my analysis results are wrong: > > > > > > In [92]: > nii=nb.load('../REST/REST_fix/101915/MNINonLinear/Results/rfMRI_REST1_LR/rfMRI_REST1_LR_Atlas_hp2000_clean.dtseries.nii') > > > > In [93]: left_ts=nb.gifti.read('101915.left_ts.gii') > > > > In [94]: nii.shape > > Out[94]: (1, 1, 1, 1, 91282, 1200) > > > > In [95]: len(left_ts.darrays) > > Out[95]: 1200 > > > > In [96]: left_ts.darrays[0].data.shape > > Out[96]: (32492,) > > > > In [97]: from lxml import etree > > > > In [98]: cifti = etree.XML(nii.header.extensions[0].get_content()) > > > > In [99]: cifti[0][2][0].items() > > Out[99]: > > [('IndexOffset', '0'), > > ('IndexCount', '29696'), > > ('ModelType', 'CIFTI_MODEL_TYPE_SURFACE'), > > ('BrainStructure', 'CIFTI_STRUCTURE_CORTEX_LEFT'), > > ('SurfaceNumberOfNodes', '32492')] > > > > In [100]: cifti[0][2][0][0].text[:100] > > Out[100]: '0 1 2 3 4 5 6 8 9 10 11 12 13 14 15 16 17 18 19 20 68 69 70 > > 71 72 73 74 75 76 77 78 79 80 81 82 83 8' > > > > In [113]: idxs=np.array(cifti[0][2][0][0].text.split(' ')).astype(np.int > ) > > > > # so ok until now it seems that the 29696 first columns of the cifti > > file should be the one from the left cortex grayordinates. > > > > In [102]: left_tss = np.c_[[d.data for d in left_ts.darrays]] > > In [103]: left_tss.shape > > Out[103]: (1200, 32492) > > > > tss = np.squeeze(nii.get_data()) > > > > In [105]: tss.shape > > Out[105]: (91282, 1200) > > > > In [106]: np.allclose(tss[0],left_tss[:,0]) > > Out[106]: False > > > > In [109]: np.allclose(tss[:1200,0],left_tss[:,0]) > > Out[109]: True > > > > In [110]: np.allclose(tss[1200:2400,0],left_tss[:,1]) > > Out[110]: True > > > > #now I do compare the all the data from left cortex (hard reshaping > > the data) excluding those which are not stored in the cifti (how are > > there magically restored when doing cifti-separate is a mystery??) > > > > In [117]: > np.allclose(tss.T.reshape(tss.shape)[:29696],left_tss[:,idxs].T) > > Out[117]: True > > > > so forcing the slicing in the other direction gave me the right data, > > so it seems there is a problem either: > > - in the way the cifti library is using the nifti2 format to store the > > array dimensions > > - or there is a mismatch between row/column-major order storage of the > data. > > Do you have any insight on this problem? > > > > Sorry again but another question is whether the cifti format/nifti2 > > implementation used in HCP opened to public? > > If not, do you plan to release it to the community or it is intended > > to remain closed source? > > > > Many thanks for your help. > > > > Cheers. > > > > Basile Pinsard > > > > > > On Thu, Dec 5, 2013 at 4:35 PM, Glasser, Matthew > > <[email protected]> wrote: > >> So there is a bug in the .dlabel.nii files that left and right ROIs are > not > >> kept separate (because FreeSurfer gives them the same name in both > >> hemispheres). I have just fixed this for the next data release, but it > will > >> have to be worked around here by processing left and right hemispheres > >> separately. My suggestion for extracting timecourses from the ROIs is > to > >> use this approach: > >> > >> wb_command -cifti-create-label > >> > ${StudyFolder}/${Subject}/MNINonLinear/fsaverage_LR32k/${Subject}.L.aparc.a2009s.32k_fs_LR.dlabel.nii > >> -left-label > >> > ${StudyFolder}/${Subject}/MNINonLinear/fsaverage_LR32k/${Subject}.L.aparc.a2009s.32k_fs_LR.label.gii > >> -roi-left > >> > ${StudyFolder}/${Subject}/MNINonLinear/fsaverage_LR32k/${Subject}.L.atlasroi.32k_fs_LR.shape.gii > >> wb_command -cifti-create-label > >> > ${StudyFolder}/${Subject}/MNINonLinear/fsaverage_LR32k/${Subject}.R.aparc.a2009s.32k_fs_LR.dlabel.nii > >> -right-label > >> > ${StudyFolder}/${Subject}/MNINonLinear/fsaverage_LR32k/${Subject}.R.aparc.a2009s.32k_fs_LR.label.gii > >> -roi-right > >> > ${StudyFolder}/${Subject}/MNINonLinear/fsaverage_LR32k/${Subject}.R.atlasroi.32k_fs_LR.shape.gii > >> > >> For left and right dlabel files, parcellate the dense timeseries to get > >> average timecourses in each ROI: > >> wb_command -cifti-parcellate <dtseries> <dlabel> COLUMN <ptseries> > >> > >> To get the average ROI timeseries out of the parcellated timeseries > file: > >> wb_command -nifti-information <ptseries> -print–matrix > >> > >> The alternative, if you were more comfortable with GIFTI files would be > to > >> split the CIFTI file into components: > >> > >> wb_command -cifti-separate <dtseries> -metric CORTEX_LEFT <func.gii> > -metric > >> CORTEX_RIGHT <func.gii>, would would produce GIFTI timeseries at which > >> point, you could directly use the GIFTI label files. > >> > >> Peace, > >> > >> Matt. > >> > >> From: basile pinsard <[email protected]> > >> Date: Thursday, December 5, 2013 3:03 AM > >> To: Matt Glasser <[email protected]> > >> Cc: "[email protected]" <[email protected]> > >> Subject: Re: [HCP-Users] Q3 rfMRI > >> > >> Hi Matt, > >> > >> yes I mean homotopic connections. > >> The individual dense connectivity map in the image you sent is not > available > >> from the HCP DB, am I wrong? > >> Was this computed from dense timeseries ICA-FIX as provided online > using the > >> wb_command -cifti-correlation? Does it add global signal regression to > get > >> the negative parts? > >> > >> Maybe the way I extract the signal from cortical ROIS is wrong, but I > >> heavily rely on the cifti header to determine which rows belong to each > ROIs > >> so I used the "brain models" to get the IndexOffset and IndexCount of > the > >> different brain models, then used the > >> ??????.?.aparc.a2009s.32k_fs_LR.label.gii to get the ROIs in the surface > >> brain models using only the vertices listed in the > CIFTI_STRUCTURE_CORTEX_* > >> brain model NodeIndices xml element (which certainly exclude the filled > >> holes in surface and fs medial wall ROI). > >> for instance in the 101915 first scan LR , there are the following brain > >> models: > >> > >> model type structure > >> index count > >> CIFTI_MODEL_TYPE_SURFACE CIFTI_STRUCTURE_CORTEX_LEFT 0 > 29696 > >> CIFTI_MODEL_TYPE_SURFACE CIFTI_STRUCTURE_CORTEX_RIGHT 29696 > >> 29716 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_ACCUMBENS_LEFT > 59412 > >> 135 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_ACCUMBENS_RIGHT > 59547 > >> 140 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_AMYGDALA_LEFT > 59687 > >> 315 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_AMYGDALA_RIGHT > 60002 > >> 332 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_BRAIN_STEM 60334 > >> 3472 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_CAUDATE_LEFT 63806 > >> 728 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_CAUDATE_RIGHT > 64534 > >> 755 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_CEREBELLUM_LEFT > 65289 > >> 8709 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_CEREBELLUM_RIGHT 73998 > >> 9144 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_L > >> 83142 706 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_DIENCEPHALON_VENTRAL_R > >> 83848 712 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_HIPPOCAMPUS_LEFT 84560 > >> 764 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_HIPPOCAMPUS_RIGHT > 85324 > >> 795 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_PALLIDUM_LEFT > 86119 > >> 297 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_PALLIDUM_RIGHT > 86416 > >> 260 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_PUTAMEN_LEFT 86676 > >> 1060 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_PUTAMEN_RIGHT > 87736 > >> 1010 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_THALAMUS_LEFT > 88746 > >> 1288 > >> CIFTI_MODEL_TYPE_VOXELS CIFTI_STRUCTURE_THALAMUS_RIGHT > 90034 > >> 1248 > >> > >> but the problem might be in the mapping of ROIs from the > aparc.a2009s+aseg > >> gifti files. > >> > >> Sorry if it is not clear but I had to reverse engineer the stuff so I > may > >> have done it wrong and it might not be possible to do it without some > prior > >> hidden in the workbench. > >> > >> Thanks for you insight on this. > >> > >> Basile > >> > >> > >> On Thu, Dec 5, 2013 at 12:52 AM, Glasser, Matthew < > [email protected]> > >> wrote: > >>> > >>> I'm not sure what you mean by homologous connectivity. Does that mean > >>> connectivity between corresponding parts of the hemispheres (like is > shown > >>> in the attached png from subject 101915—seed vertex is white, > corresponding > >>> vertex is green)? > >>> > >>> These data are acquired with a very fast TR (0.720ms), which means that > >>> the movement parameters may pick up physiological motion effects in > addition > >>> to actual subject movements. In slower TR dat sets these are aliased > into > >>> the timeseries so you cannot see them. > >>> > >>> Peace, > >>> > >>> Matt. > >>> > >>> From: basile pinsard <[email protected]> > >>> Date: Wednesday, December 4, 2013 4:24 PM > >>> To: "[email protected]" <[email protected]> > >>> Subject: [HCP-Users] Q3 rfMRI > >>> > >>> Hi all, > >>> > >>> trying to process individual rfMRI from HCP, we found unexpected > results > >>> of connectivity: > >>> - almost no homologous connectivity (see attachment) either by taking > mean > >>> timeseries in aparc.a2009s+aseg rois(as store in fsaverage32k gifti > files > >>> for surfaces) or by computing the 30gigs voxelwise dense matrix zscore > and > >>> averaging in the same rois. > >>> We did it with and without regressing out mean gray timecourse from the > >>> same dense timeseries because it seems to influence largely the > >>> connectivity. > >>> We also did in on the base preprocessed, and ica-fixed data to check, > >>> results are quite similar. > >>> Will the pipeline used for preprocessing be released publicly at some > >>> point? > >>> > >>> - the motion estimated is quite shaky and unrealistic (for example the > >>> scans of subject 101915 enclosed), do you know which method was used > for > >>> realignment of the data? If this estimated motion is used both for > >>> interpolation of data and regression of nuisance signal, this has > certainly > >>> huge effect on the connectivity. > >>> > >>> Was the same data/preprocessing used for estimating the 40subject dense > >>> correlation matrix ? This has expected connectivity structure. > >>> > >>> Thanks for your help. > >>> Best. > >>> > >>> Basile Pinsard. > >>> > >>> _______________________________________________ > >>> HCP-Users mailing list > >>> [email protected] > >>> http://lists.humanconnectome.org/mailman/listinfo/hcp-users > >>> > >>> > >>> > >>> ________________________________ > >>> > >>> The materials in this message are private and may contain Protected > >>> Healthcare Information or other information of a sensitive nature. If > you > >>> are not the intended recipient, be advised that any unauthorized use, > >>> disclosure, copying or the taking of any action in reliance on the > contents > >>> of this information is strictly prohibited. If you have received this > email > >>> in error, please immediately notify the sender via telephone or return > mail. > >>> > >>> _______________________________________________ > >>> HCP-Users mailing list > >>> [email protected] > >>> http://lists.humanconnectome.org/mailman/listinfo/hcp-users > >> > >> > >> > >> > >> ________________________________ > >> > >> The materials in this message are private and may contain Protected > >> Healthcare Information or other information of a sensitive nature. If > you > >> are not the intended recipient, be advised that any unauthorized use, > >> disclosure, copying or the taking of any action in reliance on the > contents > >> of this information is strictly prohibited. If you have received this > email > >> in error, please immediately notify the sender via telephone or return > mail. > >> > >> _______________________________________________ > >> HCP-Users mailing list > >> [email protected] > >> http://lists.humanconnectome.org/mailman/listinfo/hcp-users > > _______________________________________________ > HCP-Users mailing list > [email protected] > http://lists.humanconnectome.org/mailman/listinfo/hcp-users > _______________________________________________ HCP-Users mailing list [email protected] http://lists.humanconnectome.org/mailman/listinfo/hcp-users
