USG/Rozwiązania

# -*- coding: UTF-8 -*-
## (c) Tomasz Steifer 2016
## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/Klasyczna_rekonstrukcja

import numpy as np
import pylab as py

import scipy.signal as ss

f0=5.5e6 # Transducer center frequency [Hz]
fs=50e6 # Sampling frequency [Hz]


c=1490. # Speed of sound [m/s]
pitch=0.00021 #Transducer's pitch
NT=192 # Number of full aperture elements
Nrec=64 #Number of receive events/lines
Ntr=64 # Transmit subaperture
Rf1 = 40*(1e-3) #depth of focal point

fname="usg1_bf_nitki.npy" #numpy file with RF data from Ultrasound Research Platform

RF=np.load(fname)
Mz=np.shape(RF)[1]
##reminder about data:
##first dimension -> receive transducers //Nrec
##second dimension -> time of acquisition ~ depth
##third dimension -> transmit event  // (NT - Ntr)

#First, we need to generate transmit delays - to be used in reconstruction

transmit_delays = -(np.arange(Ntr)-(float(Ntr)/2-0.5))*pitch
transmit_delays = (np.sqrt(transmit_delays**2+Rf1**2)-Rf1)/c*fs

#For sake of clarity, beamforming is performed in loops
image=np.zeros((NT-Ntr,Mz)) #
tmp=np.zeros((Mz,Ntr))
for line in range(NT-Ntr):
    for k in range(Nrec):
        tmp0=RF[line+k,transmit_delays[k]:Mz+transmit_delays[k],line]
        tmp[:len(tmp0),k]=tmp0
    tmp2=np.mean(tmp,axis=1)
    image[line,:len(tmp2)]=tmp2

#Some simple filters
b, a = ss.butter(10, 0.05, 'highpass')    
image = ss.filtfilt(b, a, image,axis=1)
b, a = ss.butter(10, 0.7, 'lowpass')    
image = ss.filtfilt(b, a, image,axis=1)


#We need to interpolate the image
xx=np.linspace(-(NT - Ntr)/2.*pitch,((NT - Ntr))/2.*pitch,(NT - Ntr)) #old x scale
zz=np.linspace(0,np.shape(image)[1]/2*(1./fs)*c,np.shape(image)[1]) #old depth scale
from scipy import interpolate
#f=interpolate.interp2d(xx,zz,image,kind='cubic')
f=interpolate.interp2d(xx,zz,image.transpose(),'cubic')

#new (natural) grid
step=c/(2*f0)
step=step/4
maxdepth=np.max(zz)
mindepth=0.005 #try to set minimal depth to 0 and see what happens
maxWidth=0.01
minWidth=0.01

X0 = np.arange(-minWidth,maxWidth,step)#.reshape((-1,1))
R0 = np.arange(mindepth,maxdepth,step)

#and interpolation
image=f(X0,R0)

#Lastly, envelope detection using Hilbert transform
image=np.abs(ss.hilbert(image,axis=0))

low=50 #lower bound for dB scale

#and log-conversion etc <- this may also be done by appropriate matplotlib methods
indices = image>0
indices2 =image<0
indices = indices+indices2
image=np.abs(image)/np.max(np.abs(image[indices]))
image[indices]=10*np.log(image[indices])

indices = indices <-low
image[indices]=-low
indices = image>=0
image[indices]=0

#Let's have a look:
py.imshow(image,cmap="gray")
py.show()

# -*- coding: UTF-8 -*-
## (c) Tomasz Steifer 2016
## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/PWI

import numpy as np
import pylab as py

import scipy.signal as ss

f0=5.5e6 # Transducer center frequency [Hz]
fs=50e6 # Sampling frequency [Hz]


c=1490. # Speed of sound [m/s]
pitch=0.00021 #Transducer's pitch
NT=192 # Number of full aperture elements
Nrec=192 #Number of receive events/lines
Ntr=192 # Transmit subaperture
Rf1 = 1000 #depth of focal point *-1
Theta = [ -1.74532925e-01 , -1.39626340e-01,  -1.04719755e-01,  -6.98131701e-02,
  -3.49065850e-02,  -2.77555756e-17,   3.49065850e-02 ,  6.98131701e-02,
   1.04719755e-01  , 1.39626340e-01  , 1.74532925e-01] #Transmit angles

fname="usg2_pwi_nitki.npy" #numpy file with RF data from Ultrasound Research Platform

RF=np.load(fname)
RF=np.transpose(RF,[1,0,2])
Mz=np.shape(RF)[1]
##reminder about data:
##first dimension -> receive transducers //Nrec
##second dimension -> time of acquisition ~ depth
##third dimension -> transmit event  // (NT - Ntr)

#First, we need to generate transmit delays - to be used in reconstruction
transmit_delays = np.zeros((len(Theta), Ntr))

for m in xrange(len(Theta)):
    for n in xrange(Ntr):         
        Xf = -1*Rf1*np.sin(Theta[m]) - ( n - (Ntr/2-0.5) )*pitch
        transmit_delays[m,n] = (np.sqrt(Xf**2 + (Rf1*np.cos(Theta[m]))**2) - Rf1)
    transmit_delays[m,:] = transmit_delays[m,:]

#Function for reconstruction
def pwi_shift(matrix,gridx,gridz,theta,fs,c,pitch,startdepth):
    #Reconstruction for one angle
    Nrec=np.shape(matrix)[1]
    out=np.zeros((len(gridz),len(gridx),Nrec))
    gridx,gridz = np.meshgrid(gridx,gridz,indexing='xy')
    shift_tr=(gridx*np.sin(theta)+gridz*np.cos(theta) )/c*fs
    for nrec in xrange(Nrec):
        Rec_delay=np.sqrt(gridz**2+(gridx-(-(float(Nrec)/2-1./2)*pitch + nrec*pitch))**2)/c
        shift_rec=Rec_delay*fs
        indx=(shift_tr+shift_rec+startdepth).astype(np.int)
        out[:,:,nrec]=matrix[:,nrec][indx]
    return out

def pwi(matrix,gridx,gridz,angles,fs,c,pitch):
    out=np.zeros((len(gridz),len(gridx),len(angles)))
    for angle in xrange(len(angles)):
        startdepth=1/2.*(np.max(np.abs(transmit_delays[angle,:]))*fs/c)+15 #start sample shift
        shifted=pwi_shift(matrix[:,:,angle],gridx,gridz,angles[angle],fs,c,pitch,startdepth)
        indices=np.abs(shifted)>0
        divs=np.ones(np.shape(shifted))
        divs[indices]=shifted[indices]/shifted[indices]
        divs=np.sum(divs,axis=2)
        out[:,:,angle]=np.sum(shifted,axis=2)/divs
    return out

#grid
step=c/(2*f0) 
step=step/4
maxdepth=(np.shape(RF)[0]*c/fs) *0.3 
mindepth=0.005 #try to set minimal depth to 0 and see what happens
maxWidth=0.01
minWidth=0.01

X0 = np.arange(-minWidth,maxWidth,step)#.reshape((-1,1))
R0 = np.arange(mindepth,maxdepth,step)

#Some simple filters
b, a = ss.butter(3, (f0-(0.5*f0))/(0.5*fs), 'highpass')    
RF = ss.filtfilt(b, a, RF,axis=0)
b, a = ss.butter(3, (f0+(0.5*f0))/(0.5*fs), 'lowpass')    
RF = ss.filtfilt(b, a, RF,axis=0)

image = np.mean(pwi(RF,X0,R0,Theta,fs,c,pitch),axis=2)


#Lastly, envelope detection using Hilbert transform
image=np.abs(ss.hilbert(image,axis=0))

low=50 #lower bound for dB scale

#and log-conversion etc <- this may also be done by appropriate matplotlib methods
indices = image>0
indices2 =image<0
indices = indices+indices2
image=np.abs(image)/np.max(np.abs(image[indices]))
image[indices]=10*np.log(image[indices])

indices = indices <-low
image[indices]=-low
indices = image>=0
image[indices]=0

#Let's have a look:
py.imshow(image,cmap="gray")
py.show()

## (c) Tomasz Steifer 2016
## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/Doppler
import numpy as np
import pylab as py
import scipy.signal as sig




def demodulate(array,f=5.5e6, ar = 3.079291762894534e-05 ):
    t=np.arange(np.shape(array)[0])*(2.*ar/1540)
    e=np.exp(-1j*2*np.pi*f*t)
    niuarray=np.zeros(np.shape(array))+0j
    for y in range(np.shape(array)[1]):
        for z in range(np.shape(array)[2]):
            niuarray[:,y,z]=array[:,y,z]*e
    return niuarray



def autocorr_freq(array):
    r=np.sum(array[1:]*np.conj(array[:-1]))
    return np.arctanh(r.imag/r.real)
    

def doppler(array):
    newarray=np.zeros(np.shape(array)[:2])
    for x in range(np.shape(array)[0]):
        for y in range(np.shape(array)[1]):
            newarray[x,y]=autocorr_freq(array[x,y,:])
    return newarray

def moda(array):
    return np.average(np.arange(len(array)),weights=array)



big=np.load("usg3_doppler_circle_10lh.npy")


big=big[:,:,:]

IQbig=demodulate(big)
import scipy.signal as sig
b, a = sig.butter(3, 5.5/9, 'lowpass')    
IQbig = sig.filtfilt(b, a, IQbig,axis=0)

low=60
Vmax=np.abs(big)
indices = Vmax>0
indices2= Vmax<0
indices= indices+indices2
Vmax=np.abs(Vmax)/np.max(np.abs(Vmax[indices]))
Vmax[indices] = 10*np.log(Vmax[indices])

indices = Vmax <-low
Vmax[indices]=-low
indices = Vmax >= 0
Vmax[indices]=np.min(Vmax)
RFbig=Vmax

print np.shape(RFbig)

py.subplot(2,1,1)
py.imshow(RFbig[:,:,0])
py.subplot(2,1,2)
flow=doppler(IQbig)
py.imshow(sig.medfilt2d(flow))
py.show()

## (c) Tomasz Steifer 2016
## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/Parametryczne
#Input: RF images taken from two different angles
def korelacje(ar1,ar2,d):
    a,b=np.shape(ar1)
    out=np.zeros(np.shape(ar1))
    for x in xrange(a):
        for y in xrange(d/2,b-d-1):
            out[x,y]=np.argmax(np.correlate(ar1[x,y:y+d],ar2[x,y:y+d],mode="same"))
    return out

## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/GPU
kernels="
        __kernel void funkcja(__global float *a, __global float *b, __global float *out,)
        {
        //zaczynamy od zadeklarowania zmiennych
        const int n = get_global_id(0); //indeks w pierwszym wymiarze danej jednostki roboczej
        const int m = get_global_id(1); //indeks w drugim wymiarze danej jednostki roboczej
        const int M = get_global_size(1); 
        const int nm = n*M + m; // indeks w tablicy wyjsciowej odpowiadajacy wartości o współrzędnych (n,m)
        a[nm] = ceil(sqrt(a[nm]*a[nm]));
        b[nm] = ceil(sqrt(b[nm]*b[nm]));
        int c;
        while ( a != 0 ) {
          c = a[nm]; a[nm] = b[nm]%a[nm];  b[nm] = c;
        }
        out[nm]=b[nm];
        }
"