<commit_before>import numpy as np
import nibabel as nib
from dipy.data import fetch_stanford_hardi, read_stanford_hardi, get_sphere
from dipy.reconst.shm import CsaOdfModel, normalize_data
from dipy.reconst.odf import peaks_from_model

fetch_stanford_hardi()
img, gtab = read_stanford_hardi()


data = img.get_data()
affine = img.get_affine()
print('data.shape (%d, %d, %d, %d)' % data.shape)


mask = data[..., 0] > 50

mask_small  = mask[20:50,55:85, 38:40]
data_small  = data[20:50,55:85, 38:40]

csamodel = CsaOdfModel(gtab, 4, smooth=0.006)
#csa_fit = csa_model.fit(data)

sphere = get_sphere('symmetric724')
csapeaks = peaks_from_model(model=csamodel,
                            data=data_small,
                            sphere=sphere,
                            relative_peak_threshold=.8,
                            min_separation_angle=45,
                            mask=mask_small,
                            return_odf=False,
                            normalize_peaks=True)

GFA = csapeaks.gfa

print('GFA.shape (%d, %d, %d)' % GFA.shape)
nib.save(nib.Nifti1Image(GFA.astype('float32'), affine), 'gfa.nii.gz')    



from dipy.data import get_sphere
sphere = get_sphere('symmetric724')

from dipy.viz import fvtk
r = fvtk.ren()

odfs = csamodel.fit(data_small[:,:,1:2]).odf(sphere)
fvtk.add(r, fvtk.sphere_funcs(odfs, sphere, colormap='jet'))
fvtk.show(r)
fvtk.clear(r)

indices = np.where(odfs < 0)
odfs[indices] = 0
r = fvtk.ren()
fvtk.add(r, fvtk.sphere_funcs(odfs, sphere, colormap='jet'))
fvtk.show(r)


# min-max normalization
csa_mm = (odfs - np.min(odfs, -1)[..., None]) / (np.max(odfs, -1) - np.min(odfs, -1))[..., None]
r = fvtk.ren()
fvtk.add(r, fvtk.sphere_funcs(odfs, sphere, colormap='jet'))
fvtk.show(r)
fvtk.clear(r)

<commit_msg>added csa_odf_example_gfa_sh.py to highlight the problem with the csa odf model sh coefficients and gfa computed from it<commit_after>import numpy as np
import nibabel as nib
from dipy.data import fetch_stanford_hardi, read_stanford_hardi, get_sphere
from dipy.reconst.shm import CsaOdfModel, QballModel, normalize_data
from dipy.reconst.odf import gfa, odf_remove_negative_values, minmax_normalize

fetch_stanford_hardi()
img, gtab = read_stanford_hardi()


data = img.get_data()
affine = img.get_affine()
print('data.shape (%d, %d, %d, %d)' % data.shape)


mask = data[..., 0] > 50

mask_small  = mask[20:50,55:85, 38:40]
data_small  = data[20:50,55:85, 38:40]

csamodel = CsaOdfModel(gtab, 4, smooth=0.006)
csa_fit = csamodel.fit(data_small)

sphere = get_sphere('symmetric724')
csa_odf = csa_fit.odf(sphere)
gfa_csa = gfa(csa_odf)

odfs = odf_remove_negative_values(csa_odf)
gfa_csa_wo_zeros = gfa(odfs)

csa_mm = minmax_normalize(odfs) 
gfa_csa_mm = gfa(csa_mm)

qballmodel = QballModel(gtab, 6, smooth=0.006)
qball_fit = qballmodel.fit(data_small)
qball_odf = qball_fit.odf(sphere)
gfa_qball = gfa(qball_odf)
gfa_qball_mm = gfa(minmax_normalize(qball_odf))


print 'Saving GFAs...'
nib.save(nib.Nifti1Image(gfa_qball.astype('float32'), affine), 'gfa.nii.gz')    
nib.save(nib.Nifti1Image(gfa_qball_mm.astype('float32'), affine), 'gfa_mm.nii.gz')    
nib.save(nib.Nifti1Image(gfa_csa.astype('float32'), affine), 'gfa_csa.nii.gz')    
nib.save(nib.Nifti1Image(gfa_csa_wo_zeros.astype('float32'), affine), 'gfa_csa_wo_neg.nii.gz')    
nib.save(nib.Nifti1Image(gfa_csa_mm.astype('float32'), affine), 'gfa_csa_mm.nii.gz')    



coeff = csa_fit._shm_coef
print 'Printing min-max of the zeroth order SH coeff of the CSA odf'
print np.min(coeff[:,:,:,0]),np.max(coeff[:,:,:,0])
gfa_sh = np.sqrt(1.0 - (csa_fit._shm_coef[:,:,:,0] ** 2 / ( np.sum(np.square(csa_fit._shm_coef), axis=3) ) ) )
gfa_sh[np.isnan(gfa_sh)] = 0
nib.save(nib.Nifti1Image(gfa_sh.astype('float32'), affine), 'gfa_sh_csa.nii.gz')    

coeff_qb = qball_fit._shm_coef
gfa_sh_qb = np.sqrt(1.0 - (coeff_qb[:,:,:,0] ** 2 / ( np.sum(np.square(coeff_qb), axis=3) ) ) )
gfa_sh_qb[np.isnan(gfa_sh_qb)] = 0
nib.save(nib.Nifti1Image(gfa_sh_qb.astype('float32'), affine), 'gfa_sh_qball.nii.gz')    


{'fraction_non_alphanumeric': 0.15308151093439365, 'fraction_numerical': 0.031312127236580514, 'mean_word_length': 4.296052631578948, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3390579', 'n_tokens_mistral': 1882, 'n_tokens_neox': 1766, 'n_words': 284}