[scikit-learn] using a mask for brain images
Alle Meije Wink
a.m.wink at gmail.com
Fri May 5 06:13:08 EDT 2017
Thanks for that Gael - I do know nilearn but in this case I did a
depth-first search on doing SVM on brain images and ended up here :)
The size of 'coef' is (1,205739), the size of mask[mask].T is (205739,)
Same number of elements, different storage layout(?).
Turned out that the numpy.ravel() function -pointed out by my colleague-
can solve that!
>>> wmap[mask]=np.ravel(coef)
bw
Alle Meije
On 5 May 2017 at 07:59, Gael Varoquaux <gael.varoquaux at normalesup.org>
wrote:
> Hi Alle,
>
> I think that what has changed between 2014 and today is that the
> coefficients coef are now a 2D array (number of hyperplanes x number of
> features). In your case, the first direction is of length one, so you
> could just do:
>
> coef = clf.coef_[0]
>
> and your script should work.
>
> The code of the Abraham paper cannot be updated, because papers are not
> living objects. However, we are maintaining a package that encodes all
> these patterns in higher-level construct: http://nilearn.github.io/
> It might be a good idea to use this package, as it is maintained and has
> quality assurance.
>
> Best,
>
> Gaël
>
> On Thu, May 04, 2017 at 05:52:27PM +0200, Alle Meije Wink wrote:
> > I have a script to classify MRI perfusion maps from healthy subjects and
> > patients. For the file IO and the classifier I have started with the
> example
> > code in Abraham et al 2014 [https://arxiv.org/pdf/1412.3919.pdf].
>
> > I use the same classifier as in the paper to produce a back-projected
> map of
> > classification weights, which I then want to 'unmask' like in the paper:
>
> > coef=clf.coef_
> > coef=featureselection.inverse_transform(coef)
>
> > and
>
> > map_name='weights_check.nii.gz'
> > wmap=np.zeros(mask.shape, dtype=X.dtype)
> > wmap[mask]=coef
> > img=nb.Nifti1Image(wmap,np.eye(4))
> > img.to_filename(map_name)
>
> > But the line "wmap[mask]=coef" throws an error "ValueError: boolean
> index array
> > should have 1 dimension". I tried the example code from the paper and
> that
> > works.
>
> > Is the 'coef' array of back-projected SVM weights in some way different
> than
> > the masked input image? Or am I doing something else wrong? The error
> suggests
> > that the mask array is the problem.
>
> > The complete script is attached.
>
> > Many thanks for your help!
> > Alle Meije
>
> > # -*- coding: utf-8 -*-
> > """
> > classify MR images from controls (CON) and patients (MCI)
>
> > """
>
> > import os
> > import numpy as np
> > import nibabel as nb
> > import sklearn as sl
>
> > from sklearn.feature_selection import f_classif
> > from sklearn import svm
>
> > def subj_lists( rootdir='/data/a.wink/MR', starts=['CON','MCI'] ):
>
> > # make and empty list for each patient group
> > slists=[]
> > for i in range(0,len(starts)):
> > slists.append([]);
>
> > # add subjects to the group based on
> > for root, dirs, files in os.walk(rootdir):
> > for fname in files:
> > for i in range(0,len(starts)):
> > if fname.startswith(starts[i]):
> > slists[i].append(os.path.join(root,fname))
>
> > print('finished building subject lists')
> > return slists
>
> > def mk_mask ( subj_lists=[] ):
>
> > # check whether mask can be loaded
> > mname='mask_check.nii.gz'
>
> > if os.path.isfile(mname):
>
> > msk=nb.load(mname).get_data()
>
> > # or must be built (values >20% in the subject images and >20% grey
> matter in the template)
> > else:
>
> > num_im=0;
> > for i in range(0,len(subj_lists)):
> > for fname in subj_lists[i]:
> > imdata=nb.load(fname).get_data()
> > if 'msk' not in locals():
> > msk=np.absolute(imdata)
> > else:
> > msk=msk+np.absolute(imdata)
> > num_im+=1
>
> > msk/=num_im
> > msk/=np.amax(msk)
> > msk[msk<.2]=0
> > msk[msk>0]=1
>
> > grey=nb.load('/usr/local/spm8/apriori/grey.nii').get_data()
> > grey[grey<.2]=0
> > grey[grey>0]=1
>
> > msk=msk*grey
>
> > img=nb.Nifti1Image(msk,np.eye(4))
> > img.to_filename(mname)
>
> > msk=msk.astype(bool)
>
> > print('finished building mask')
> > return msk, mname
>
> > def load_images( subj_lists=[], mask=[] ):
>
> > # load all the images, build a matrix X or subjects (rows) *
> inmask-voxels (columns)
> > num_im=0;
> > for i in range(0,len(subj_lists)):
> > for fname in subj_lists[i]:
> > imdata=nb.load(fname).get_data()
> > num_im+=1
> > print '\r' + str(num_im) +'\r',
> > if 'the_matrix' not in locals():
> > the_matrix=imdata[mask].T
> > else:
> > the_matrix=np.vstack((the_matrix,imdata[mask].T))
>
> > print('finished building matrix X')
> > return the_matrix
>
> > def main():
>
> > # get the input data
> > subjlists = subj_lists()
> > y=np.concatenate( ([-1 for _ in range(len(subjlists[0]))],
> > [ 1 for _ in range(len(subjlists[1]))]) )
> > mask, mas_name = mk_mask(subjlists)
> > X = load_images(subjlists,mask)
> > print "number of subjects, voxels: %d, %d" % X.shape
>
> > # select features
> > featureselection=sl.feature_selection.SelectKBest(f_classif, k=8000)
> > X_reduced=featureselection.fit_transform(X,y)
>
> > # classify
> > clf=svm.SVC(kernel='linear')
> > clf.fit(X_reduced,y)
>
> > # make discrimination map
> > coef=clf.coef_
> > coef=featureselection.inverse_transform(coef)
>
> > map_name='weights_check.nii.gz'
> > wmap=np.zeros(mask.shape, dtype=X.dtype)
> > wmap[mask]=coef
> > img=nb.Nifti1Image(wmap,np.eye(4))
> > img.to_filename(map_name)
>
> > if __name__ == "__main__":
> > main()
>
>
> > _______________________________________________
> > scikit-learn mailing list
> > scikit-learn at python.org
> > https://mail.python.org/mailman/listinfo/scikit-learn
>
>
> --
> Gael Varoquaux
> Researcher, INRIA Parietal
> NeuroSpin/CEA Saclay , Bat 145, 91191 Gif-sur-Yvette France
> Phone: ++ 33-1-69-08-79-68
> http://gael-varoquaux.info http://twitter.com/GaelVaroquaux
> _______________________________________________
> scikit-learn mailing list
> scikit-learn at python.org
> https://mail.python.org/mailman/listinfo/scikit-learn
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scikit-learn/attachments/20170505/e9145ba4/attachment-0001.html>
More information about the scikit-learn
mailing list