<commit_before>from scipy.io import loadmat
import numpy as np
#?loadmat

phantom=loadmat('/home/eg309/Desktop/phantom_test_data.mat',struct_as_record=True)
all=phantom['all']
b_table=phantom['b_table']
odf_vertices=phantom['odf_vertices']
odf_faces=phantom['odf_faces']
s = all[14,14,1]
l_values=np.sqrt(b_table[0]*0.01506)
tmp=np.tile(l_values,(3,1))
b_vector=b_table[1:4,:]*tmp
q2odf=np.sinc(np.dot(b_vector.T, odf_vertices) * 1.2/np.pi)
odf=np.dot(s,q2odf)






<commit_msg>Peak finding and QA done<commit_after>from scipy.io import loadmat
import numpy as np
#?loadmat

phantom=loadmat('/home/eg309/Desktop/phantom_test_data.mat',struct_as_record=True)
all=phantom['all']
b_table=phantom['b_table']
odf_vertices=phantom['odf_vertices']
odf_faces=phantom['odf_faces']

#s = all[14,14,1]

#Yeh et.al, IEEE TMI, 2010
#calculate the odf using GQI

scaling=np.sqrt(b_table[0]*0.01506) # 0.01506 = 6*D where D is the free
# water diffusion coefficient 
# l_values sqrt(6 D tau) D free water
# diffusio coefficiet and tau included in the b-value

tmp=np.tile(scaling,(3,1))

b_vector=b_table[1:4,:]*tmp

Lambda = 1.2 # smoothing parameter - diffusion sampling length

q2odf_params=np.sinc(np.dot(b_vector.T, odf_vertices) * Lambda/np.pi) # implements equation no. 9 from Yeh et.al.


def Q2odf(s,q2odf_params):


    odf=np.dot(s,q2odf_params)

    return odf

def peak_finding(odf,odf_faces,odf_vertices):

    #proton density already include from the scaling b_table[0][0] and s[0]


    #find local maxima

    peak=odf.copy()

    # where the smallest odf values in the vertices of a face remove the
    # two smallest vertices 

    for face in odf_faces.T:

        i, j, k = face

        check=np.array([odf[i],odf[j],odf[k]])

        zeroing=check.argsort()

        peak[face[zeroing[0]]]=0

        peak[face[zeroing[1]]]=0

    #for later testing expecting peak.max 794595.94774980657 and
    #np.where(peak>0) (array([166, 347]),)


    #we just need the first half of peak

    peak=peak[0:len(peak)/2]

    #find local maxima and give fiber orientation (inds) and magnitute
    #peaks in a descending order

    inds=np.where(peak>0)[0]

    pinds=np.argsort(peak[inds])
    
    peaks=peak[inds[pinds]][::-1]


    return peaks, inds[pinds][::-1]




'''

#s = all[14,14,1]

s = all[0,0,0]

#s = all[]
        
odf = Q2odf(s,q2odf_params)

peaks,inds=peak_finding(odf,odf_faces,odf_vertices)

'''



S=all

x,y,z,g=S.shape

S=S.reshape(x*y*z,g)

QA = np.zeros((x*y*z,5))

IN = np.zeros((x*y*z,5))

fwd = 0

#Calculate Quantitative Anisotropy and find the peaks and the indices
#for every voxel

for (i,s) in enumerate(S):

    odf = Q2odf(s,q2odf_params)

    peaks,inds=peak_finding(odf,odf_faces,odf_vertices)

    fwd=max(np.max(odf),fwd)

    peaks = peaks - np.min(odf)

    l=min(len(peaks),5)

    QA[i][:l] = peaks[:l]

    IN[i][:l] = inds[:l]

    
QA/=fwd

#If you just want to generate a volume then is enough to use the first QA[0]


#Generate streamlines














{'fraction_non_alphanumeric': 0.14442961974649768, 'fraction_numerical': 0.03835890593729153, 'mean_word_length': 4.100340136054422, 'pattern_counts': {'":': 0, '<': 3, '<?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': '1027476', 'n_tokens_mistral': 1370, 'n_tokens_neox': 1240, 'n_words': 253}