|
|
| Linia 1: |
Linia 1: |
| − | <source lang=python>
| |
| − | # -*- 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()
| |
| − |
| |
| − | </source>
| |
| − |
| |
| − | <source lang=python>
| |
| − | # -*- 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()
| |
| − |
| |
| − |
| |
| − | </source>
| |
| − |
| |
| − | <source lang=python>
| |
| − | ## (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()
| |
| − |
| |
| − | </source>
| |
| − |
| |
| − | <source lang=python>
| |
| − | ## (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
| |
| − | </source>
| |
| − |
| |
| − | <source lang=python>
| |
| − | ## 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];
| |
| − | }
| |
| − | "
| |
| − | </source>
| |
