Tsteifer: UWAGA! Usunięcie treści (strona pozostała pusta)!

2016-07-29T08:57:55Z

UWAGA! Usunięcie treści (strona pozostała pusta)!

Tsteifer: Utworzono nową stronę " # -- coding: UTF-8 -- ## (c) Tomasz Steifer 2016 ## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/Klasyczna_rekonstrukcja import numpy a..."

2016-07-29T08:55:34Z

Utworzono nową stronę "<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 a..."

Nowa strona

<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>

Rozwiązania - Historia wersji

Tsteifer: UWAGA! Usunięcie treści (strona pozostała pusta)!

Tsteifer: Utworzono nową stronę " # -*- coding: UTF-8 -*- ## (c) Tomasz Steifer 2016 ## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/Klasyczna_rekonstrukcja import numpy a..."

Tsteifer: Utworzono nową stronę " # -- coding: UTF-8 -- ## (c) Tomasz Steifer 2016 ## solutions to http://brain.fuw.edu.pl/edu/index.php/USG/Klasyczna_rekonstrukcja import numpy a..."