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
[email protected]
https://mail.python.org/mailman/listinfo/scikit-learn
