<jupyter_start><jupyter_text>test<jupyter_code>I.shape

directory1='/home/mahdi/python codes/sct_test/evaluation/crop_moco_mean/crop_moco_mean.nii'
I = sitk.ReadImage(directory1)
I = sitk.GetArrayFromImage(I)
I=I.astype(np.float32)
q=I[3,:,:]
q[27,2]=0
#q1=q.rotate(180, expand=True)
#q1=cv2.rotate(q,180)

plt.imshow(q,cmap='gray')
plt.show()
def GaussianFunction(x, sigma):
    if sigma == 0:
        return 0
    else:
        g = (1/math.sqrt(2*math.pi*sigma*sigma))*math.exp(-x*x)/(2*sigma*sigma)
    return g

# function returns the gaussian kernel using the GaussianFunction of size 3x3
def GaussianMask(sigma):
    g = []
    for i in range(-2, 3):#creating a gaussian kernel of size 3x3
        g1 = GaussianFunction(i,sigma)
        g2 = GaussianFunction(i-0.5, sigma)
        g3 = GaussianFunction(i+0.5, sigma)
        gaussian = (g1+g2+g3)/3
        g.append(gaussian)
    return g

sigma = 1.5
G = [] # Gaussian Kernel
G = GaussianMask(sigma)

def DownSample(I): 
    Ix = Iy  = []
    I = np.array(I)
    S = np.shape(I)      #shape of the image
    for i in range(S[0]):
        Ix.extend([signal.convolve(I[i,:],G,'same')])#convolution of the I[i] with G
    Ix = np.array(np.matrix(Ix))
    Iy = Ix[::2, ::2]#selects the alternate column and row
    return Iy

def UpSample(I):
    I = np.array(I)
    S = np.shape(I)

    Ix = np.zeros((S[0], 2*S[1]))#inserting alternate rows of zeros
    Ix[:, ::2] = I
    S1 = np.shape(Ix)
    Iy = np.zeros((2*S1[0], S1[1]))#inserting alternate columns of zeros
    Iy[::2, :] = Ix
    Ig = cv2.GaussianBlur(Iy, (5,5), 1.5, 1.5)#instead of using the user-defined gaussian function, I am using the Gaussian Blur functtion for double the size of gaussian kernel size
    return Ig
def LucasKanadeee(I1, I2,d,q):
    I1 = np.array(I1)
    I2 = np.array(I2)
    S  = np.shape(I1)
    
    smooth_kernel=5
    I1_smooth = cv2.GaussianBlur(I1,(smooth_kernel,smooth_kernel), 0)
    I2_smooth = cv2.GaussianBlur(I2 ,(smooth_kernel,smooth_kernel), 0)
    
    q=1/8
    
    Ix = signal.convolve2d(I1_smooth,[[-q,q ] , [-q,q]],'same') + signal.convolve2d(I2_smooth,[[-q,q],[-q,q]] ,'same')
    Iy = signal.convolve2d(I1_smooth,[[-q,-q],   [q,q]],'same') + signal.convolve2d(I2_smooth ,[[-q,-q],[q,q]]  ,'same')
    It = signal.convolve2d(I1_smooth,[[q,q  ]  , [q,q]],'same') + signal.convolve2d(I2_smooth ,[[-q,-q],[-q,-q]],'same')
    
    
    features = cv2.goodFeaturesToTrack(I1_smooth, 1000, q, d)
    features = np.int0(features)
      

    
 
    u = np.nan*np.ones(S)
    v = np.nan*np.ones(S)
    
    for l in features:
        j,i = l.ravel()
        IX = ([Ix[i-1,j-1],Ix[i,j-1],Ix[i+1,j+1],Ix[i-1,j],Ix[i,j],Ix[i+1,j],Ix[i-1,j+1],Ix[i,j+1],Ix[i+1,j-1]])
        IY = ([Iy[i-1,j-1],Iy[i,j-1],Iy[i+1,j+1],Iy[i-1,j],Iy[i,j],Iy[i+1,j],Iy[i-1,j+1],Iy[i,j+1],Iy[i+1,j-1]])
        IT = ([It[i-1,j-1],It[i,j-1],It[i+1,j+1],It[i-1,j],It[i,j],It[i+1,j],It[i-1,j+1],It[i,j+1],It[i+1,j-1]])
        #IX = ([Ix[i-1,j-1],Ix[i-1,j],Ix[i-1,j+1],Ix[i,j-1],Ix[i,j],Ix[i,j+1],Ix[i+1,j-1],Ix[i+1,j],Ix[i+1,j+1]])
        #IY = ([Ix[i-1,j-1],Ix[i-1,j],Ix[i-1,j+1],Ix[i,j-1],Ix[i,j],Ix[i,j+1],Ix[i+1,j-1],Ix[i+1,j],Ix[i+1,j+1]])
        #IT = ([Ix[i-1,j-1],Ix[i-1,j],Ix[i-1,j+1],Ix[i,j-1],Ix[i,j],Ix[i,j+1],Ix[i+1,j-1],Ix[i+1,j],Ix[i+1,j+1]])
        # Using the minimum least squares solution approach
        LK = (IX,IY)
        LK = matrix(LK)
        LK = array(LK)
       
        #Psedudo Inverse
     
        A2 = np.linalg.pinv(LK)
        IT = np.array(IT)
        b  = -IT
        (u[i,j],v[i,j]) = np.dot(A2.T,b) # we have the vectors with minimized square error
    
    return u,v
def LucasKanadeIterative(I1, I2, u1, v1,kernel):
   
    I1 = np.array(I1)
    I2 = np.array(I2)
    S  = np.shape(I1)
    u1 = np.round(u1)
    v1 = np.round(v1)
     
 
    u = np.nan*np.ones(S)
    v = np.nan*np.ones(S)

   
    for i in range(int(kernel/2),S[0]-int(kernel/2)):
        for j in range(int(kernel/2),S[1]-int(kernel/2)):
 
            I1new = I1[i-int(kernel/2):i+int(kernel/2)+1,j-int(kernel/2):j+int(kernel/2)+1]# picking 5x5 pixels at a time

            lr = (i-int(kernel/2))+v1[i,j]#Low Row Index
            hr = (i+int(kernel/2))+v1[i,j]#High Row Index
            lc = (j-int(kernel/2))+u1[i,j]#Low Column Index
            hc = (j+int(kernel/2))+u1[i,j]#High Column Index
            #window search and selecting the last window if it goes out of bounds
            if(lr < 0): 
                lr = 0
                hr = kernel-1
            if(lc < 0): 
                lc = 0
                hc = kernel-1
            if(hr > (len(I1[:,0]))-1): 
                lr = len(I1[:,0])-kernel
                hr = len(I1[:,0])-1                
            if(hc > (len(I1[0,:]))-1): 
                lc = len(I1[0,:])-kernel
                hc = len(I1[0,:])-1               
            if(np.isnan(lr)):
                lr = i-int(kernel/2)
                hr = i+int(kernel/2)
            if(np.isnan(lc)):
                lc = j-int(kernel/2)
                hc = j+int(kernel/2)
            #Selecting the same window for the second frame
            I2new = I2[int(lr):int((hr+1)),int(lc):int((hc+1))]
            # Now applying LK for each window of the 2 images
            IX = signal.convolve2d(I1new,[[-0.25,0.25],[-0.25,0.25]],'same') + signal.convolve2d(I2new,[[-0.25,0.25],[-0.25,0.25]],'same')
            IY = signal.convolve2d(I1new,[[-0.25,-0.25],[0.25,0.25]],'same') + signal.convolve2d(I2new,[[-0.25,-0.25],[0.25,0.25]],'same')
            IT = signal.convolve2d(I1new,[[0.25,0.25],[0.25,0.25]]  ,'same') + signal.convolve2d(I2new,[[-0.25,-0.25],[-0.25,-0.25]],'same')
            
            #if kernel>1:
            #IX = np.transpose(IX[1:kernel,1:kernel])
            #IY = np.transpose(IY[1:kernel,1:kernel])
            #IT = np.transpose(IT[1:kernel,1:kernel])

            
            IX = IX.ravel()
            IY = IY.ravel()
            IT = IT.ravel()

            LK = (IX,IY)
            LK = matrix(LK)
            LK = array(LK)
            
            #Psedudo Inverse

            A2 = np.linalg.pinv(LK)
            IT = np.array(IT)
            b  = -IT
            (u[i,j],v[i,j]) = np.dot(A2.T,b) # we have the vectors with minimized square error
              
          
    return u,v

def pyramid (Im1,Im2):

    #Im1=mask(I,158,5)
    #Im2=mask(I,159,5)

    #Im1=plt.imread(dir1)
    #Im2=plt.imread(dir2)

    #Im1=i1
    #Im2=j1


    I1 = np.array(Im1)
    I2 = np.array(Im2)
    S = np.shape(I1)

    pyramid1 = np.empty((S[0],S[1],3)) 
    pyramid2 = np.empty((S[0],S[1],3)) 
    pyramid1[:,:,0] = I1 #since the lowest level is the original image
    pyramid2[:,:,0] = I2 #since the lowest level is the original image

    #creating the pyramid by downsampling the original image
    for i in range(1,3):
        I1 = DownSample(I1)
        I2 = DownSample(I2)
        pyramid1[0:np.shape(I1)[0], 0:np.shape(I1)[1], i] = I1
        pyramid2[0:np.shape(I2)[0], 0:np.shape(I2)[1], i] = I2
    kernel=5
    #####################################################################################################################    

    k=0
    level=3
    I1new = pyramid1[0:int(len(pyramid1[:,0])/(2**(level-k-1))),0:int(len(pyramid1[0,:])/(2**(level-k-1))),level-k-1]
    I2new = pyramid2[0:int(len(pyramid2[:,0])/(2**(level-k-1))),0:int(len(pyramid2[0,:])/(2**(level-k-1))),level-k-1]

    #u,v=LucasKanadeee(mask(I,158,5),mask(I,159,5),1,.0001)
    u,v=LucasKanadeee(Im1,Im2,1,.0001)

    for i in range(0, 3):
        (u,v) = LucasKanadeIterative(I1new, I2new, u, v,kernel)



    u_l0 = u
    v_l0 = v
    I_l0 = I1new

    #u_l0[np.where(u_l0 == 0)] = nan
    #v_l0[np.where(v_l0 == 0)] = nan


    ####################################################################################################################
    #for level 1
    k = 1
    u = UpSample(u)+
    v = UpSample(v)
    I1new = pyramid1[0:int(len(pyramid1[:,0])/(2**(level-k-1))),0:int(len(pyramid1[0,:])/(2**(level-k-1))),level-k-1]
    I2new = pyramid2[0:int(len(pyramid2[:,0])/(2**(level-k-1))),0:int(len(pyramid2[0,:])/(2**(level-k-1))),level-k-1]
    (u,v) = LucasKanadeIterative(I1new, I2new, u, v,kernel)

    u_l1 = u
    v_l1 = v
    I_l1 = I1new
    #u_l1[np.where(u_l1 == 0)] = nan
    #v_l1[np.where(v_l1 == 0)] = nan

    ########################################################################################################
    k = 2
    u = UpSample(u)
    v = UpSample(v)
    I1new = pyramid1[0:int(len(pyramid1[:,0])/(2**(level-k-1))),0:int(len(pyramid1[0,:])/(2**(level-k-1))),level-k-1]
    I2new = pyramid2[0:int(len(pyramid2[:,0])/(2**(level-k-1))),0:int(len(pyramid2[0,:])/(2**(level-k-1))),level-k-1]
    #for i in range(0, 3):

    (u,v) = LucasKanadeIterative(I1new, I2new, u, v,kernel)

    S1=I1new.shape
    for i in range(S1[0]):
        for j in range(S1[1]):
             if norm((u[i,j],v[i,j]))<.05  or norm((u[i,j],v[i,j]))>2:
                    u[i,j]=0
                    v[i,j]=0
    u_l2 = u
    v_l2 = v
    I_l2 = I1new

    #u_l2[np.where(u_l2 == 0)] = nan
    #v_l2[np.where(v_l2 == 0)] = nan
    
    
    
    nr, nc = Im1.shape
    
    row_coords, col_coords = np.meshgrid(np.arange(nr), np.arange(nc),
                                         indexing='ij')
    im1_warp = warp(I1new, np.array([row_coords + v_l2, col_coords + u_l2]),
                       order=0) 
    return im1_warp
import nibabel as nib
import numpy as np 
from scipy import ndimage, misc
import time
import os
import subprocess
#####################
directory='/home/mahdi/python codes/sct_test/evaluation/crop before shift/shift after crop .nii'
I = sitk.ReadImage(directory)
I = sitk.GetArrayFromImage(I)
I=I.astype(np.float32)
#####################
start_time = time.time()
#==============================
img  = nib.load(directory)
img_mask_affine = img.affine
#################################
header = img.header
nb_img = header.get_data_shape()
nb_img_h = nb_img[0] #Hauteur
#################################
o=np.zeros((nb_img[0],nb_img[1],nb_img[2],nb_img[3]))
#o=np.zeros((110,110,nb_img[2],nb_img[3]))
kernel=5
for v in range(0,nb_img[3]):
     for s in range(0,nb_img[2]):
            #a=pyramid(I[v,s,:,:],I[159,s,:,:])
            a=pyramid(I[v,s,:,:],I[159,s,:,:])
            o[:, :, s,v] = a.T
o=o.astype(np.int16)

print("--- %s second ---" % (time.time() - start_time))
img_reg = nib.Nifti1Image(o, affine=img_mask_affine, header=header)
nib.save(img_reg,'/home/mahdi/python codes/motion result/shift after crop edit '+str(kernel))
subprocess.Popen(['fsleyes','/home/mahdi/python codes/motion result/shift after crop edit'+str(kernel)]) ## just change the output names.
395.1838300228119/60<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.24028629856850717, 'fraction_numerical': 0.05623721881390593, 'mean_word_length': 2.6845890410958906, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '221772', 'license': 'MIT', 'n_tokens_mistral': 4958, 'n_tokens_neox': 4664, 'n_words': 703, 'path': 'v7.ipynb', 'repo_name': 'mahdimplus/Optical-flow', 'size': '10758'}