Page 1 of 51 Photoacoustic Image Reconstruction ENPH 479 Final Report Kevin Zhou Zongyi Zhang Sean Yang Project Sponsor: Dr. Shuo Tang Group Number: 1264 Engineering Physics The University of British Columbia Jan 7, 2013 Page 2 of 51 Executive Summary Our project?s main objective is to design and improve an accurate, reliable image reconstruction algorithm using the Delay and Sum method. The Fast Fourier Transform algorithm is also investigated to evaluate the efficiency and performance of our choice of algorithm. The developed MATLAB programs are to be able to take in the pre-beam-formed data acquired by the ultrasound system and reconstruct proper images of the planned seeds in the phantom. The developed algorithms should work to produce images with a satisfactory resolution for both homogeneous and inhomogeneous materials. Photoacoustic imaging is a fairly modern technique in the medical imaging filed, and it?s under rapid development. While the principle of photoacoustic effect and ultrasound wave production has been discovered more than a century ago, its potential application as a medical tool has really only come into development after the common use of ultrasound imaging and the availability of advanced ultrasound transducers. Our algorithms aim to reconstruct images using acquired data with a satisfactory resolution. Speed is also one consideration, as the algorithms have the potential to be integrated into a future system to achieve real-time imaging. Clarity of reconstruction judged from a qualitative standpoint, numerical analysis of signal to noise ratio (SNR), and some considerations of computation time were discussed in this document. Simulated data and experimental data, generated using phantoms constructed of metal seeds, were both used in analyzing the performance of our algorithms. Our findings suggest that while the computation time for both algorithms under identical situations is reasonably short for practical use, the FFT algorithm that we developed has shown significant SNR performance advantage over the Delay and Sum algorithm for homogeneous medium, indicating higher image quality. However, the FFT algorithm cannot reconstruct images using data obtained with inhomogeneous medium. Therefore, the Delay and Sum algorithm is a more viable image reconstruction algorithm to be used in various practical cases with medium not limited to homogenous materials. Page 3 of 51 Contents 1 Introduction..................................................................................................................... 7 1.1 Historical background and motivation ..................................................................... 7 1.2 Previous imaging modalities ..................................................................................... 8 1.3 Current state of Photoacoustic Imaging ................................................................... 8 1.4 Scope and limitations ............................................................................................... 9 1.5 Project Objectives ..................................................................................................... 9 2 Discussion ...................................................................................................................... 11 2.1 Photoacoustic Effect ............................................................................................... 11 2.2 Experiment Set-up .................................................................................................. 11 2.3 Image Reconstruction Methods ............................................................................. 13 2.4 Raw Ultrasound Data Spec ..................................................................................... 17 2.5 Delay and Sum algorithm ....................................................................................... 17 2.5.1 Delay and sum photoacoustic image reconstruction results .......................... 17 2.5.2 Test with simulation data ................................................................................ 17 2.5.3 Test with experimental data............................................................................ 20 2.6 Fast Fourier Transform (FFT) Algorithm ................................................................. 27 2.6.1 FFT Reconstruction Theory .............................................................................. 27 2.6.2 Test with Simulation Data ............................................................................... 28 2.6.3 Test with experimental data............................................................................ 29 3 Conclusions.................................................................................................................... 32 3.1 Results..................................................................................................................... 32 3.2 Challenges ............................................................................................................... 33 4 Project Deliverables ...................................................................................................... 35 4.1 List of Deliverables .................................................................................................. 35 4.2 Financial Summary ................................................................................................. 35 4.3 Ongoing Commitments by Team Members ........................................................... 35 5 Recommendations ........................................................................................................ 36 Page 4 of 51 6 Appendices .................................................................................................................... 37 6.1 Signal to Noise Ratio Matlab Code ......................................................................... 37 6.2 Delay and Sum Algorithm for Homogeneous Medium .......................................... 37 6.3 Delay and Sum Algorithm for Inhomogeneous Medium ....................................... 40 6.4 FFT Implementation Homogeneous ....................................................................... 44 6.5 Simulated Data Generation Code ........................................................................... 47 7 References ..................................................................................................................... 50 Page 5 of 51 List of Figures Figure 1 Experimental Setup ................................................................................................ 12 Figure 2 Ultrasound Data Acquisition Experiment Setup ..................................................... 12 Figure 3 The Delay and Sum algorithm ................................................................................ 14 Figure 4 The FFT Algorithm ................................................................................................... 15 Figure 5 The FFT algorithm for inhomogeneous medium .................................................... 16 Figure 6 Two seeds in a homogeneous medium - sensor data ............................................ 18 Figure 7 Two seeds in a homogeneous medium - reconstructed image ............................. 19 Figure 8 Two seeds in a homogeneous medium - actual positions of the seeds defined in the simulation ....................................................................................................................... 20 Figure 9 Three seeds in a homogeneous medium - sensor data ......................................... 21 Figure 10 Three seeds in a homogeneous medium - reconstructed image ......................... 22 Figure 11 Four seeds in a homogeneous medium - actual seeds positions ......................... 23 Figure 12 Four seeds in a homogeneous medium ? sensor data ......................................... 24 Figure 13 Four seeds in a homogeneous medium - reconstructed image ........................... 25 Figure 14 Nine seeds randomly inserted in an inhomogeneous medium - sensor data ..... 26 Figure 15 Nine seeds randomly inserted in an inhomogeneous medium - reconstructed image .................................................................................................................................... 27 Figure 16 Delay compensation for FFT Theory ..................................................................... 28 Figure 17 Two seeds in a homogeneous medium ................................................................ 29 Page 6 of 51 List of Tables Table 1 The Comparison of three algorithms ....................................................................... 16 Table 2 Raw Ultrasound Data Spec ...................................................................................... 17 Table 3 Algorithm SNR Performance Comparison ............................................................... 32 Table 4 Computation Time Comparison ............................................................................... 32 Page 7 of 51 1 Introduction This project, sponsored by Dr. Shuo Tang from the Electrical Engineering Department at the University of British Columbia, is to develop a working algorithm for image reconstruction by the means of photoacoustic imaging (PAI). This document discusses the goals, methods and deliverables of a medical imaging project concerned primarily with image reconstruction of ultrasound data generated through the relatively unused and lesser known imaging methodology of PAI. In this section, some background information on a medical imaging methodology, photoacoustic tomography, and the motivation behind our project is provided. Historical background of the technique is first discussed, and the current state of photoacoustic imaging and project objectives are also introduced. 1.1 Historical background and motivation Nowadays, scientific instruments have now been used extensively in the health care industry, in providing diagnostic analysis, treatments and many other benefits. Medical imaging of patients? body to improve the accuracy of diagnoses for human disease is becoming of critical importance. Photoacoustic imaging, also known as photoacoustic tomography, is a modern medical imaging technique under rapid development. While the principle of photoacoustic effect and ultrasound wave production has been discovered more than a century ago, its potential application as a medical tool has really only come into development after the common use of ultrasound imaging and the availability of advanced ultrasound transducers. It was only in the past few decades that lasers became capable of producing fast enough pulses to meet certain conditions with regard to thermal and acoustic confinement to make photoacoustic tomography a plausible method. The main motivation behind this project is to use photoacoustic tomography in the early detection of prostate cancer. Prostate cancer is one of the most commonly seen diagnosis in male cancer patients, therefore a wide range of research is being conducted in the medical community regarding diagnosis and treatment. Choice of diagnosis tests and the need for treatment in various cases and stages of prostate cancer is still an area of much disagreement, with the primary reason being specific unique biological qualities of the cancer. Page 8 of 51 Prostate brachytherapy is a popular prostate cancer treatment option that involves the permanent implantation of radioactive seeds into the prostate. However, contemporary brachytherapy procedure is limited by the lack of imaging systems that can provide real-time seed-position feedback. While many other imaging systems have been proposed, photoacoustic imaging has emerged as a potential ideal modality to address this need, since it could easily be incorporated into the current ultrasound system used in the operating rooms. Although still in its infancy, the initial results of photoacoustic imaging systems from a few research groups for the application of prostate brachytherapy seed localization are highly promising. 1.2 Previous imaging modalities Currently, the most commonly used imaging method for prostate cancer diagnosis is Transrectal Ultrasonography (TRUS), which is used primarily for determining the prostate gland volume rather than imaging of the cancer cells or differentiating healthy and cancerous tissue. The current TRUS technology alone has become generally accepted as insufficient to generate the seed position information due to its poor image quality. Magnetic resonance imaging (MRI) is also used for prostate cancer imaging, enabling a much more detailed visualization of prostate tissue. However, there has been evidence that post-biopsy recovery of tissue is worse with MRI than TRUS. Meanwhile, the cost (and resulting lack of availability) of MRI makes it somewhat prohibitive in comparison with ultrasound techniques. In addition, the resolution and contrast of photoacoustic tomography may be able to detect prostate cancer at a far earlier stage, which is almost always without initial symptoms. Current techniques, such as digital rectal examination (DRE) and serum level evaluation, are still incapable of detecting early stage cancers. A faster and more sensitive immediate diagnosis may be made with photoacoustic tomography. 1.3 Current state of Photoacoustic Imaging Photoacoustic imaging is currently in use for applications related to mining and geological applications, such as rock analysis. In this project, we are only concerned with its application in medical imaging. Authors such as Alexander Oraevsky (Fairway Medical Technologies) and Lihong Wang (Washington University, St. Louis) have extensively Page 9 of 51 experimented with photoacoustic medical applications and proven that it holds promise as a medical imaging modality. This technology is centered on a physical phenomenon known as the photoacoustic effect, in which electromagnetic energy, such as pulsed light from a laser source, is absorbed by a target material, converted into acoustic energy through rapid thermoelastic expansion, and can consequently be detected with ultrasound transducers. PAI would be effective in localizing seeds due to the strong contrast it generates between seeds and tissue at certain wavelengths of light. It would also be practical for BT since it takes advantage of already-existing ultrasound equipment in the OR, thus reducing any additional burden on the surgeon. PAI would also be cost-effective, since it would not require any large imaging equipment such as a CT or MRI scanner, but only the addition of a laser. Other current developments of photoacoustic imaging include focusing lasers so that only a small localized area is illuminated. Scanning over a series of points provides much greater detail and accuracy, but it clearly requires extensive scanning which may be time intensive. 1.4 Scope and limitations This project focuses on designing and improving the Delay and Sum algorithm. The Fast Fourier Transform algorithm is also investigated to evaluate the efficiency and performance of our choice of algorithm. Several methods of image improvement and actual testing results using both simulated data and acquired experiment data from Dr. Shuo Tang?s Microsystem and Nanotechnology Lab are also discussed. Both homogeneous medium and inhomogeneous medium are considered when developing the imaging algorithm. The physical apparatus used for the experiment, which has been developed and used by previous team in research on photoacoustic imaging, is covered. 1.5 Project Objectives The main objective of our project is to establish and test two accurate, reliable image reconstruction algorithms using the Delay and Sum method and FFT method. The Page 10 of 51 developed MATLAB programs are to be able to take in the pre-beam-formed data acquired by the ultrasound system and reconstruct proper images of the planned seeds in the phantom. When an initial working algorithm is developed, the second primary objective is to optimize the algorithm to produce images with a higher resolution. The testing of the performance of developed algorithms starts with a simplified case ? producing a photoacoustic image of the metal seeds in homogenous media. Once satisfactory results are obtained with using homogenous media, the project will then focus on testing using inhomogeneous media, which is a more likely case in practical applications. At the end of the project, a complete set of working image reconstruction algorithms will be delivered, along with this final project report discussing the design and performance of these algorithms in details. Page 11 of 51 2 Discussion 2.1 Photoacoustic Effect The photoacoustic effect is a conversion between light and acoustic waves due to absorption and localized thermal excitation. When rapid pulses of light are incident on a sample of matter, they can be absorbed and the resulting energy will then be radiated as heat. This heat causes detectable sound waves due to pressure variation in the surrounding medium. The acoustic wave is measured by an ultrasound probe and the raw data recorded can then be used for image reconstruction. Three promising algorithms for photoacoustic image reconstruction are: Delay and Sum, FFT and Time Reversal. Although Time Reversal has the potential to produce the highest image quality, it is currently not suitable for real-time image reconstruction due to its high computational cost (20 minutes computation time for one frame). Therefore, in our project, we will be focusing on developing the algorithms using Delay and Sum as well as FFT. 2.2 Experiment Set-up The experimental apparatus for data acquisition, although not a major concern for the purpose of this project, is briefly introduced in this section. The hardware of the photoacoustic imaging system consists of three major parts: the laser, the ultrasound and the data acquisition and processing system, as shown in the following illustration. Page 12 of 51 Figure 1 Experimental Setup Figure 2 Ultrasound Data Acquisition Experiment Setup Page 13 of 51 In the photoacoustic phenomenon, light is absorbed by a material and converted to heat, and the subsequent thermo elastic expansion generates an acoustic wave. The acoustic wave is measured by an ultrasound probe and the raw data recorded can then be used for image reconstruction. 2.3 Image Reconstruction Methods The image reconstruction algorithm will be implemented entirely in MATLAB. Delay and Sum algorithm will be our primary focus for this project. However, k-Wave Matlab Toolbox offers image reconstruction routines (FFT & Time reversal versions only) which can generate in theory more accurate photoacoustic images than Delay and Sum. The Delay and Sum algorithm starts from a number of acoustic measurements and combines them (Sum) taking into account the different acoustic paths (Delay). The hypothesis behind this algorithm is that a sound wave emitted at a certain time from a source is perceived by an observer, placed at a distance R from the source, after a time delay ?t=R/c, where c is the sound speed in the considered medium. The same signal undergoes an amplitude reduction which, supposing a spherical propagation model can be considered linear with the distance. Thus it is possible to reconstruct the original signal by applying a proper delay and amplitude correction. If the source position is unknown, a measurement method can be built to locate the source itself. Page 14 of 51 Figure 3 The Delay and Sum algorithm The Fast Fourier Transform algorithm is an alternative to the delay and sum. This method is fundamentally different from the delay and sum as it performs a fast Fourier transform to frequency domain on the time step axis and a transform to wavenumber domain in the sensor position axis. Then a mapping from wavenumber-frequency domain to wavenumber-wavenumber domain is applied based on the dispersion relationship of a plane wave in homogeneous medium. The final step is to do an inverse Fourier Transform from wavenumber-wavenumber space to get an x-y spatial domain image. These transforms and mapping are all done via the k-Wave toolbox in MATLAB. Page 15 of 51 Figure 4 The FFT Algorithm In the development of the FFT algorithm, we experimented with method to implement this algorithm for inhomogeneous medium. Since in the case of inhomogeneous medium, the wave speed (i.e. the wave number) would be different a direct FFT would not work. Since we would know the medium density distribution from the ultra sound data, we would be able to use this information in developing an inhomogeneous FFT reconstruction. In the simplest case, we would have 2 layer of different homogeneous medium. In this scenario, it would be possible to use an approximate FFT reconstruction by first reconstructing the pressure wave to the point prior to crossing the boundary between the two layers using a constant wave speed. Then we perform a second FFT with a different wave speed to reconstruct the point source. Page 16 of 51 Figure 5 The FFT algorithm for inhomogeneous medium This idea was never implemented because of the realization that the algorithm would only work for single point sources. If multiple sources are present, it would be impossible to reconstruct each wave to the media boundary simultaneously. We abandoned this development due to time limitations and decided to focus more on the more reliable delay and sum algorithm for inhomogeneous medium. Another reconstruction considered is time reversal approach. It numerically solves the wave equation in reverse order from the pressure recorded detectors outwards. It can be used to account for heterogeneous density and speed of sound in the medium. However, it is impossible to implement for real time imaging due to its ridiculous high computational time and cost, which takes approximately 20 minutes to process 1 frame. It is one of the more accurate reconstruction methods, but it?s not a suitable choice in actual clinical setting. Table 1 The Comparison of three algorithms Delay & Sum FFT Time Reversal Homogeneous Yes Yes Yes Inhomogeneous Yes No Yes Page 17 of 51 The performance of the algorithms is measured by signal-to-noise ratio (SNR). The signal strength of a reconstructed image is chosen to be the highest signal magnitude in the constructed image. All the pixels with magnitude less than half of the signal strength are considered as noise. The SNR is calculated by dividing the signal strength by the mean of the noise. Although there exists many other ways to define the SNR, this definition is effective for our purpose. 2.4 Raw Ultrasound Data Spec Table 2 Raw Ultrasound Data Spec Sampling Frequency 80 MHz Total Delay 24 ms Resolution 128 x 1560 Ultrasound Sensor Spacing 0.3048 mm Ultrasound Sensor Number (Number of Channels) 128 Sound speed in the Phantom 1.54E6 mm/s 2.5 Delay and Sum algorithm 2.5.1 Delay and sum photoacoustic image reconstruction results The delay and sum algorithm routine implemented for both homogeneous medium and inhomogeneous medium are tested with simulation data and actual measurement taken in Dr. Shuo Tang?s lab. 2.5.2 Test with simulation data 2.5.2.1 Two seeds in a homogeneous medium SNR: 4.5440 Page 18 of 51 Figure 6 Two seeds in a homogeneous medium - sensor data Page 19 of 51 Figure 7 Two seeds in a homogeneous medium - reconstructed image Page 20 of 51 Figure 8 Two seeds in a homogeneous medium - actual positions of the seeds defined in the simulation 2.5.3 Test with experimental data The actual position of the seeds was not provided by Dr. Shuo Tang?s lab. They believe the actual positions of the seeds are not important at this stage as they can always calibrate their system to generate the image with the correct seed position. The computation time on a regular recently laptop is about 1 second for homogeneous data and 2 seconds for inhomogeneous data. 2.5.3.1 Three seeds in a homogeneous medium SNR: 11.3646 Page 21 of 51 Figure 9 Three seeds in a homogeneous medium - sensor data Page 22 of 51 Figure 10 Three seeds in a homogeneous medium - reconstructed image 2.5.3.2 Four seeds in a homogeneous medium SNR: 6.9332 Page 23 of 51 Figure 11 Four seeds in a homogeneous medium - actual seeds positions Page 24 of 51 Figure 12 Four seeds in a homogeneous medium ? sensor data Page 25 of 51 Figure 13 Four seeds in a homogeneous medium - reconstructed image 2.5.3.3 Nine seeds randomly inserted in an inhomogeneous medium SNR: 3.5742 Page 26 of 51 Figure 14 Nine seeds randomly inserted in an inhomogeneous medium - sensor data Page 27 of 51 Figure 15 Nine seeds randomly inserted in an inhomogeneous medium - reconstructed image 2.6 Fast Fourier Transform (FFT) Algorithm 2.6.1 FFT Reconstruction Theory The Fast Fourier Transform (FFT) is an alternative to the delay and sum method. The advantages of FFT over delay and sum is that FFT gives a better reconstruction based on a signal-to-noise (SNR) merit as well as reconstruction speed. SNR is calculated in the same method as with the delay and sum. However, a drawback is that FFT is unable to reconstruct images in inhomogeneous medium. When implementing the FFT algorithm, we came across an issue with the delay between the laser pulse firing and the start time of the DAQ. Since the FFT algorithm does a Fourier transform from time domain to wave number domain, the delay will cause very big Page 28 of 51 inaccuracies in the reconstruction. We developed a simplistic yet effective solution to this problem by inserting artificial zero pressure data in front of the experimental data to offset the delay which solve the issue of the delay, but then the points of interests had to be extracted from the reconstructed image due to the high resolution generated by the extra data. Figure 16 Delay compensation for FFT Theory 2.6.2 Test with Simulation Data The FFT algorithm was implemented using a MATLAB toolbox, k-Wave, on homogenous medium and tested with simulated data and lab data. These cases are the same with the delay and sum method so we can compare the different methods. 2.6.2.1 Two seeds in a homogeneous medium The reconstruction was performed on the same set of data as in the Delay and Sum case, but FFT's SNR is a magnitude of order greater than Delay and Sum Page 29 of 51 SNR: 47.8960 Figure 17 Two seeds in a homogeneous medium 2.6.3 Test with experimental data The computation time on a regular recently laptop is about 0.5 second for experimental data. This is faster than Delay and Sum which takes about 1 second for the identical data and task. 2.6.3.1 Three seeds in a homogeneous medium SNR: 159.5380 Page 30 of 51 2.6.3.2 Four seeds in a homogeneous medium SNR: 127.8663 Page 31 of 51 Page 32 of 51 3 Conclusions 3.1 Results Successfully developing the photoacoustic image reconstruction algorithm detailed in this report would give our project sponsor a system that, when fully integrated with a medical ultrasound imaging machine, would be used as a platform for research regarding cancer detection and ultrasound imaging. This photoacoustic imaging system, with further development, could eventually be integrated with other imaging methods for improved medical surveillance. On a larger scale, this result of this project will be the starting point in the research and development of more efficient and effective ways of detecting not just prostate cancer, but potentially other cancers as well. The performance of the Delay and Sum algorithm and the FFT algorithm are assessed for different sets of data, and the results are listed and compared below. Table 3 Algorithm SNR Performance Comparison SNR Delay and Sum FFT Two seeds in a homogeneous medium (simulation data) 4.5440 47.8960 Three seeds in a homogeneous medium 11.3646 159.5380 Four seeds in a homogeneous medium 6.9332 127.8663 Table 4 Computation Time Comparison Delay and Sum FFT Computation time (second) 1 0.5 From the table above, we can see the FFT algorithm has significant SNR performance advantage over the Delay and Sum algorithm for homogeneous medium. The computation time for both algorithms under identical situations is reasonable for practical use. Page 33 of 51 However, the FFT algorithm cannot handle inhomogeneous medium. Therefore, the performance comparison for the inhomogeneous case cannot be done. The advantages and drawbacks of the Delay and Sum method are: Advantages o Simpler to implement o Can be modified to handle the inhomogeneous case o Fast Drawbacks o Very sensitive to input parameters (sensor sizes, sound speed, delay in the laser firing and data acquisition system, etc.). A reasonably good image quality can only be achieved when all the input parameters are very accurate (<1% error). o Sensitive to noise in raw sensor data The advantages and drawbacks of the FFT algorithm are: Advantages o Fast o Much better SNR performance under identical situations o Much more robust compared to the Delay and Sum Drawbacks o Cannot be modified to deal with inhomogeneous medium 3.2 Challenges The biggest challenge of this project is to improve the image quality of the reconstructed image for inhomogeneous medium while maintaining a reasonable computation time. The Delay and Sum algorithm is capable of producing reasonably good image quality at a very Page 34 of 51 high speed for homogeneous medium. However, for inhomogeneous medium, due to the nature of this algorithm, it can only deal with fairly simple medium structures and several simplifications have to be made in order to maintain the reasonable computation time. Thus the potential of using this algorithm for image reconstruction for inhomogeneous medium is very limited. Page 35 of 51 4 Project Deliverables Photoacoustic image reconstruction algorithms implemented in MATLAB will be delivered to the project sponsor at the end of the project along with supporting documents including documentations of the MATLAB code and this final report for the project. 4.1 List of Deliverables ? Algorithms capable of reconstructing 2D images of initial pressure ? The scripts included in the Appendices are capable of taking DAQ data acquired in the lab and producing recognizable images. These and additional scripts for signal to noise ratio calculation and simulation data generation will be delivered to the Sponsor in soft copy .m-files. ? Report of the results of the algorithm and quality of reconstruction This report serves as a discussion of these results, and includes recommendations. 4.2 Financial Summary There is no monetary cost in this project, as the development of the algorithms has been done using personal computers with MATLAB installed. The cost associated with the production of phantoms, data acquisition using Ultrasound devise and other experimental procedures are not considered here, as the scope of this project is limited to developing the image reconstruction algorithms. 4.3 Ongoing Commitments by Team Members There are no further tasks planned for the project at this moment. We will be available, however, to be consulted by future project groups on the work accomplished in this project. Page 36 of 51 5 Recommendations Both the Delay and Sum and FFT algorithms can produce satisfactory image quality at reasonably fast speed. The Fast Fourier Transform method has been shown to produce superior quality images to delay-and-sum, and is more robust. However, due to its inability to account for inhomogeneous medium, the Delay and Sum method is recommended. Time reversal can be further investigated as an imaging technique. Especially if real-time imaging is not a goal, time reversal still gives most promising reconstruction performance. Because of the wide variety of applications suitable for photoacoustics, our project could eventually lead to commercial products that adapt the method of photoacoustics for medical imaging. Since photoacoustic imaging is still in its early stages, our project will provide great contribute to the initial advancements in the field of medical imaging. Page 37 of 51 6 Appendices 6.1 Signal to Noise Ratio Matlab Code function [SNR] = signal_to_noise(img) signal=max(img(:)); theshold = signal/2; sum=0; count=0; imgg=img(:); for i=1:size(imgg,1) if imgg(i)<theshold sum=sum+imgg(i); count = count +1; end end mean=sum/count; SNR = signal/mean; end 6.2 Delay and Sum Algorithm for Homogeneous Medium %this script tests the delay and sum algorithm with simulated data %read raw data from sensors load sensorData; %plot the raw image figure; colormap(bone) imagesc(sqrt(abs(hilbert(sensorData)))) colormap(bone) title('Simulated Sensor Data') xlabel('Channel Index') ylabel('Sampling') drawnow; %pos_sensor is a 128 by 2 matrix that contains the x,y coordinates of each Page 38 of 51 %sensor in the ultrasound array. unit is in mm. L=-30; d_sensor = 0.46;%physical spacing between two adjacent sensors (mm) pos_sensor = zeros(128,2); pos_sensor(:,1)=(-d_sensor*(128-1)/2):d_sensor:(d_sensor*(128-1)/2); pos_sensor(:,2)=L; %cutoff angle for each sensor cutoff = 10; %reconstructed image resolution size_img = [300,300]; %spacing spacing = [0.2 0.2]; delay = 0; %reconstruct the image [img] = recon_ds_homogeneous(sensorData,pos_sensor,cutoff,size_img,spacing,delay); %Image display spacing_x=spacing(1); spacing_y=spacing(2); figure; colormap(hot) xax=linspace(-spacing_x*((size_img(2)-1)/2),spacing_x*((size_img(2)-1)/2),size_img(2)); yax=linspace(-spacing_y*((size_img(1)-1)/2),spacing_y*((size_img(1)-1)/2),size_img(1)); imagesc(yax,xax,img); axis square; xlabel('Lateral (mm)') ylabel('Axial (mm)') title('Delay and Sum Image Reconstruction') drawnow Page 39 of 51 %compute signal to noise ratio signal_noise_R = signal_to_noise(img) %Image Reconstruction Routine for homogeneous medium function [img] = recon_ds_homogeneous(sensorData, pos_sensor, cutoff, size_img, spacing, delay) %input: %sensorData: raw data file from the ultrasound probe. %cutoff: cutoff angle (degree) for the ultrasound probe %size_img: the size of Output image B; e.g. size_B=[300, 300] %spacing: the physical distance between two adjacent pixels in img (mm) %delay: the delay in the laser firing system (us) %output: %img is the reconstructed image %constants speed = 1.54e6;%acoustic speed in water (mm/s) s_rate = 80e6; %sampling rate of each sensor in Hz spacing_x = spacing(1); spacing_y = spacing(2); cutoff=cutoff/180*pi; ilength=size(sensorData,1); %acutal position coordinates for each pixel in img [coord_ix,coord_iy]=meshgrid((-spacing_x*(size_img(2)-1)/2):spacing_x:(spacing_x*(size_img(2)-1)/2),-(spacing_y*(size_img(1)-1)/2):spacing_y:(spacing_y*(size_img(1)-1)/2)); tic; %Image reconstruction img = zeros(size_img(1),size_img(2)); Page 40 of 51 for n = 1 : size_img(1) %image range in Y direction for m = 1 : size_img(2) %image range in X direction for j = 1 : 128 sep_x=coord_ix(n,m)-pos_sensor(j,1); sep_y=coord_iy(n,m)-pos_sensor(j,2); dis = sqrt(sep_x^2+sep_y^2); ind = ceil(s_rate * dis / (speed)-delay*1e-6*s_rate); deg=abs(atan(sep_x/sep_y)); if (ind <=ilength && ind >0 && deg<(cutoff)) img(n,m) = img(n,m) + sensorData(ind, j)*cos(deg); end end end end img =sqrt(abs(hilbert(img))); toc; 6.3 Delay and Sum Algorithm for Inhomogeneous Medium %this script tests the delay and sum algorithm with experimental data %read raw data from sensors filename='121121_mpshete5pvc'; path = strcat('C:\Users\Freebody\UBC\5.1\enpy479\src\',filename,'\'); filetag='CH000.daq'; fid = fopen([path, filetag],'r'); hdr=fread(fid,3,'int32'); numFrame=hdr(2); %number of frames contained in the data file ilength=hdr(3); %number of rows in each frame Page 41 of 51 [hdr, sensorData] = readDAQ(path, ones(1, 128), 3 , true); shallow=450; deep=1400; sensorData(1:shallow,:)=0; sensorData(deep:ilength,:)=0; sensorData=medfilt2(sensorData, [5 1]); %remove horizontal lines %plot the raw image figure; colormap(bone) imagesc(sqrt(abs(hilbert(sensorData)))) colormap(bone) title('Sensor Data') xlabel('Channel Index') ylabel('Sampling') drawnow; %pos_sensor is a 128 by 2 matrix that contains the x,y coordinates of each %sensor in the ultrasound array. unit is in mm. L=-58; d_sensor = 0.3048; %physical spacing between two adjacent sensors (mm) pos_sensor = zeros(128,2); pos_sensor(:,1)=(-d_sensor*(128-1)/2):d_sensor:(d_sensor*(128-1)/2); pos_sensor(:,2)=L; %cutoff angle for each sensor cutoff = 25; %reconstructed image resolution size_img = [300,300]; %spacing spacing = [0.08 0.08]; delay = 21; pos_layer = -48; % the y-coordinate of the boundary of two layers Page 42 of 51 waveSpeed = [1.54e6 1.3e6]; %[v1 v2] v2:=the speed of the sound in the layer next to the proble(mm/s) %reconstruct the image tic; [img] = recon_ds_inhomogeneous(sensorData,pos_sensor,pos_layer,waveSpeed,cutoff,size_img,spacing,delay); toc; %Image display spacing_x=spacing(1); spacing_y=spacing(2); figure; colormap(hot) xax=linspace(-spacing_x*((size_img(2)-1)/2),spacing_x*((size_img(2)-1)/2),size_img(2)); yax=linspace(-spacing_y*((size_img(1)-1)/2),spacing_y*((size_img(1)-1)/2),size_img(1)); imagesc(yax,xax,img); axis square; xlabel('Lateral (mm)') ylabel('Axial (mm)') title('Delay and Sum Image Reconstruction') drawnow %compute signal to noise ratio signal_noise_R = signal_to_noise(img) %Image Reconstruction Routine for inhomogeneous medium function [img] = recon_ds_inhomogeneous(sensorData,pos_sensor,pos_layer,waveSpeed,cutoff,size_img,spacing,delay) %input: %sensorData: raw data file from the ultrasound probe. %pos_sensor: the spatial coordinates of the sensors in the probe Page 43 of 51 %pos_layer: the y-coordinate of the boundary layer %wavespeed: [v1 v2] v2:=the speed of the sound in the layer next to the proble(mm/s) %cutoff: cutoff angle (degree) for the ultrasound probe %size_img: the size of Output image B; e.g. size_B=[300, 300] %spacing: the physical distance between two adjacent pixels in img (mm) %delay: the delay in the laser firing system (us) %output: %img is the reconstructed image %constants s_rate = 80e6; %sampling rate of each sensor in Hz spacing_x = spacing(1); spacing_y = spacing(2); cutoff=cutoff/180*pi; ilength=size(sensorData,1); %acutal position coordinates for each pixel in img [coord_ix,coord_iy]=meshgrid((-spacing_x*(size_img(2)-1)/2):spacing_x:(spacing_x*(size_img(2)-1)/2),-(spacing_y*(size_img(1)-1)/2):spacing_y:(spacing_y*(size_img(1)-1)/2)); %Image reconstruction img = zeros(size_img(1),size_img(2)); for n = 1 : size_img(1) %image range in Y direction for m = 1 : size_img(2) %image range in X direction for j = 1 : 128 sep_x=coord_ix(n,m)-pos_sensor(j,1); sep_y=coord_iy(n,m)-pos_sensor(j,2); dis = sqrt(sep_x^2+sep_y^2); if (coord_iy(n,m)<=pos_layer) Page 44 of 51 ind = ceil(s_rate * dis / (waveSpeed(2))-delay*1e-6*s_rate); else ind = ceil(s_rate * ((dis-pos_layer)/waveSpeed(1)+pos_layer/waveSpeed(2))-delay*1e-6*s_rate); end deg=abs(atan(sep_x/sep_y)); if (ind <=ilength && ind >0 && deg<(cutoff)) img(n,m) = img(n,m) + sensorData(ind, j)*cos(deg); end end end end img=medfilt2(img); img =sqrt(abs(hilbert(img))); 6.4 FFT Implementation Homogeneous % This function performs a FFT algorithm on image reconstruction on raw wave data for homogeneous medium function [p_xy_new] = Mod_FFT_reconstruction(filename, framenum, deep, delay) % (1) filename: string containing the name of the folder containing the % ultra sound data % (2) framenum: the frame number of the data set % (3) deep: number of time step of data to removve from the bottom % (4) delay: the time delay in data collection and firing of the laser % in microseconds Page 45 of 51 path = strcat('C:\Users\Kevin\Desktop\enphy479\',filename,'\'); % get raw data from readDAQ [hdr, RF] = readDAQ(path, ones(1, 128), framenum , true); [m_raw, n_raw] = size(RF); % define the properties of the propagation medium sound_speed = 1.54e3; % [m/s] % define physical parameters s_rate = 80e6; % sampling rate dx = 0.46/1000.0; % grid point spacing in the x direction [m] - physical distance between crystal on probe dy = sound_speed/s_rate; % grid point spacing in the y direction [m] - speed of sound divided by sampling rate % generates zero pressure data to compensate for the delay in laser % pulse and start of data collection RFdelay = zeros((s_rate*delay*1e-6),n_raw); % removes negatve space from bottom of raw image and reconstruct % pressure matrix with the compensated delay data % deep is bottom, for phantom case deep = 1500 RFnew = [RFdelay ; RF(1:(deep),:)]; [m_delay, n_delay] = size(RFdelay); [m, n] = size(RFnew); % create the computational grid kgrid = makeGrid((m-m_delay),dy,n,dx); % create the time array dt = 1/s_rate; kgrid.t_array = 0:dt:m*dt; % smooths out data (taken from Leo's code) for i=1:128 rf_temp(:,i)=gradient(RFnew(:,i),kgrid.t_array); Page 46 of 51 end gradRFf(:,:)=rf_temp; RFbf(:,:)=2*RFnew-2*gradRFf; sensor_data = RFbf; % reconstruct the initial pressure p_xy = kspaceLineRecon(sensor_data, dx, dt, sound_speed, 'PosCond', true); % extract the interested portion of the reconstruction p_xy_new = p_xy(m_delay:m,:); % plot raw data imagesc(sqrt(abs(hilbert(RF)))); colormap(bone); caxis auto; title(strcat('Raw',filename)); xlabel('Channel Index'); ylabel('Sampling'); figure; % plot the reconstructed initial pressure Img = imagesc(kgrid.x_vec*1e3,kgrid.y_vec*1e3,sqrt(abs(hilbert(p_xy)))); colormap(hot); ylabel('y-position [mm]'); xlabel('x-position [mm]'); colorbar; figure; % plot the extracted region or interest from initial pressure imagesc(kgrid.x_vec*1e3,kgrid.y_vec*1e3,sqrt(abs(hilbert(p_xy_new)))); colormap(hot); ylabel('y-position [mm]'); xlabel('x-position [mm]'); colorbar; Page 47 of 51 end 6.5 Simulated Data Generation Code % This function generates simulated pressure wave data % number of points/position/pressure intensity can be change in the body of the code function [sensor_data] = simulate_data(sampling_rate, noise, bandwidth) % (1) sampling_rate: the rate at which sampling occurs at the sensors, % vertical distance per cell calculated by 1/sampling_rate % % (2) noise: 'true' if want noise to be added to the signal % % (3) bandwidth: the bandwidth for the noise generated (smaller value = % greater noise) % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] % create the computational grid Nx = 3140; % number of grid points in the x (row) direction % 3140 Ny = 128; % number of grid points in the y (column) direction (sensors) dx = medium.sound_speed/sampling_rate; % grid point spacing in the x direction [m] total 0.05888m dy = 0.46/1000.0; % grid point spacing in the y direction [m] kgrid = makeGrid(Nx, dx, Ny, dy); %============================================================== % simulated data specification %============================================================== % create initial pressure distribution using makeDisc disc_magnitude = 5; % [au] Page 48 of 51 disc_x_pos = 1000; % [grid points] disc_y_pos = 70; % [grid points] disc_radius = 10; % [grid points] disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); disc_x_pos = 1500; % [grid points] disc_y_pos = 110; % [grid points] disc_radius = 10; % [grid points] disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius); source.p0 = disc_1 + disc_2; %============================================================== % smooth the initial pressure distribution and restore the magnitude source.p0 = smooth(kgrid, source.p0, true); % define a binary line sensor sensor.mask = zeros(Nx, Ny); sensor.mask(1, :) = 1; % create the time array dt = 1/sampling_rate; kgrid.t_array = 0:dt:dt*Nx; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor) %, input_args{:}); sensor_data = sensor_data'; % add Noise if (noise == 'true') sensor_data = addNoise(sensor_data, bandwidth); end % plot the initial pressure and sensor distribution figure; Page 49 of 51 imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, source.p0 + sensor.mask*disc_magnitude, [-disc_magnitude disc_magnitude]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; colorbar; % plot the simulated sensor data figure; imagesc(sensor_data, [-1, 1]); colormap(getColorMap); ylabel('Sensor Position'); xlabel('Time Step'); colorbar; end Page 50 of 51 7 References 1. Kuo N., Kang H., Song D.Y., Kang J.U., Boctor E.M.; Real-time photoacoustic imaging of prostate brachytherapy seeds using a clinical ultrasound system. Journal of Biomedical Optics, 17(6), 066005 (2012) 2. Cigada A., Ripamonti F., Vanali M., The delay and sum algorithm applied to microphone array measurements: Numerical analysis and experimental validation, Mechanical Systems and Signal Processing, Volume 21, Issue 6 (2007). 3. Wang, L. V., Multiscale Photoacoustic Microscopy and Computed Tomography. Nature Photonics, 3, 503-509 (2009). 4. Wang, L. V., Xu, M. Photoacoustic Imaging in Biomedicine. American Institute of Physics, 77, (2009). 5. Wang, L.V., Xu, M. Time-Domain Reconstruction for Thermoacoustic Tomography in a Spherical Geometry. IEEE Transactions on Medical Imaging, 21, 814-822 (2002). 6. Afnan, J., Tempany, C. M., Update on Prostate Imaging. Urology Clinic North America, 37, 23-25, (2010). 7. Treeby, B.E., Cox, B.T., k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. Journal of Biomedical Optics. 15, 021314-1-10, (2010). 8. Piras, A., Xia, W., Photoacoustic imaging of Breast Using the Twenty Photoacoustic Mammoscope: Present Status and Future Perspectives, IEEE Journal of Selected Topics in Quantum Electronics, 16, No.4, (2010). 9. Zhang, Hao F., Maslov, K., Stoica, G., & Wang, L. V., Functional Photoacoustic Microscopy for High-Resolution and Noninvasive In Vivo Imaging. Nature Biotechnology, 24, 848-851 (2006) 10. McDonald, F.Alan, Generalized Theory of the Photoacoustic Effect, Journal of Applied Physics, Vol 49, Issue 4, 2313-2322 11. Xu, M., Wang, L. V., Universal Back-Projection Algorithm for Photoacoustic Computed Tomography, Physical Review E, 71, 016706 (2005) 12. Wang, L. V., Prospects of Photoacoustic tomography, Medical Physics, 35(12), 5758-5767 (2008) Page 51 of 51 13. Yin, B. Z., Fast Photoacoustic Imaging System Based On 320-Element Linear Transducer Array, Physics in Medicine and Biology, 49, 1339 (2004)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Undergraduate Research /
- Photoacoustic image reconstruction
Open Collections
UBC Undergraduate Research
Photoacoustic image reconstruction Zhou, Kevin; Zhang, Zongyi; Yang, Sean Jan 7, 2013
pdf
Page Metadata
Item Metadata
Title | Photoacoustic image reconstruction |
Creator |
Zhou, Kevin Zhang, Zongyi Yang, Sean |
Date Issued | 2013-01-07 |
Description | Our project’s main objective is to design and improve an accurate, reliable image reconstruction algorithm using the Delay and Sum method. The Fast Fourier Transform algorithm is also investigated to evaluate the efficiency and performance of our choice of algorithm. The developed MATLAB programs are to be able to take in the pre-beam-formed data acquired by the ultrasound system and reconstruct proper images of the planned seeds in the phantom. The developed algorithms should work to produce images with a satisfactory resolution for both homogeneous and inhomogeneous materials. Photoacoustic imaging is a fairly modern technique in the medical imaging filed, and it’s under rapid development. While the principle of photoacoustic effect and ultrasound wave production has been discovered more than a century ago, its potential application as a medical tool has really only come into development after the common use of ultrasound imaging and the availability of advanced ultrasound transducers. Our algorithms aim to reconstruct images using acquired data with a satisfactory resolution. Speed is also one consideration, as the algorithms have the potential to be integrated into a future system to achieve real-time imaging. Clarity of reconstruction judged from a qualitative standpoint, numerical analysis of signal to noise ratio (SNR), and some considerations of computation time were discussed in this document. Simulated data and experimental data, generated using phantoms constructed of metal seeds, were both used in analyzing the performance of our algorithms. Our findings suggest that while the computation time for both algorithms under identical situations is reasonably short for practical use, the FFT algorithm that we developed has shown significant SNR performance advantage over the Delay and Sum algorithm for homogeneous medium, indicating higher image quality. However, the FFT algorithm cannot reconstruct images using data obtained with inhomogeneous medium. Therefore, the Delay and Sum algorithm is a more viable image reconstruction algorithm to be used in various practical cases with medium not limited to homogenous materials. |
Type |
Text |
Language | eng |
Series | University of British Columbia. ENPH 479 |
Date Available | 2013-11-29 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0074500 |
URI | http://hdl.handle.net/2429/45573 |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Peer Review Status | Unreviewed |
Scholarly Level | Undergraduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 52966-Zhou_K_et_al_ENPH_479_2013.pdf [ 1.42MB ]
- Metadata
- JSON: 52966-1.0074500.json
- JSON-LD: 52966-1.0074500-ld.json
- RDF/XML (Pretty): 52966-1.0074500-rdf.xml
- RDF/JSON: 52966-1.0074500-rdf.json
- Turtle: 52966-1.0074500-turtle.txt
- N-Triples: 52966-1.0074500-rdf-ntriples.txt
- Original Record: 52966-1.0074500-source.json
- Full Text
- 52966-1.0074500-fulltext.txt
- Citation
- 52966-1.0074500.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.52966.1-0074500/manifest