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
