PAT reconstruction in lossy heterogeneous medium

k-Wave User Forum » Photoacoustic and Thermoacoustic Tomography » Image Reconstruction

PAT reconstruction in lossy heterogeneous medium

(18 posts) (3 voices)

Started 13 years ago by huangchao
Latest reply from huangchao

huangchao

Member
Joined: Nov '10
Posts: 57

Hi,

For PAT reconstruction, I can get good reconstructed images in heterogeneous lossless medium or homogeneous lossy medium using time reversal algorithm, but in lossy heterogeneous medium I can not get good reconstructed images. From my understanding, one of reasons could be the data acquisition time.

On one hand, we know that time reversal requires the acquisition time to be long enough to get accurate reconstruction images in heterogeneous medium; on the other hand, the acquisition time could not be too long to avoid artifact trapping in heterogeneous medium.

However I am not sure if the dilemma of data acquisition time is really the culprit, because I can get good reconstructed images in lossless heterogeneous medium if the measured signals are long enough. Yet in lossy heterogeneous medium, I can not get good images no matter how long the signals are. (Although the image looks better when the signals were truncated according to Huygens' principle, it's still much worse than the image in lossless medium.)

So I want to know your opinion about the data acquisition time. Or is there any other reason that could deteriorate the reconstructed images in lossy heterogeneous medium?

Thanks,
Chao

Posted 13 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

Hi Chao,

Just to clarify, are you using simulated or experimental data? Also, are you working in 2D or 3D?

When you have a lossy medium there will necessarily be some information loss which will deteriorate your reconstructed image. This is for two related reasons. First, because the absorption alters the frequency spectrum of the photoacoustic waves as they propagate through the medium. Second, because for a given measurement bandwidth, the absorption will likely reduce the magnitude of some high frequency components below the noise floor.

If you are using k-Wave to perform the time reversal reconstructions you have two options.

(1) You can set your medium to be non-absorbing and accept that there will be some information loss (make sure the medium.alpha_coeff and medium.alpha_power parameters are not set otherwise you will absorb the waves twice!). This will be no worse than any of the other possible reconstruction approaches that do not compensate for acoustic attenuation.

(2) You can use attenuation compensation. In k-Wave this is done by reversing the sign of the absorption operator using medium.alpha_sign = [-1, 1], and then providing some regularisation using medium.alpha_filter. See here and here for more details, or see the Image Reconstruction With Compensation For Acoustic Absorption Example in k-Wave.

As for the particular combination of a lossy and heterogeneous medium, I can't think of any immediate answers. I will try and dig a little deeper and get back to you.

Brad.

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi Brad,

Thanks for you reply! Actually I am using "Image Reconstruction With Compensation For Acoustic Absorption Example" for simulation, and I only change the homogeneous medium into heterogeneous medium (nothing else is changed). I can not get good reconstructed images in this case, even I compensate the absorption.

However, I observe a phenomenon. That is the image looks better when the signals were truncated according to Huygens' principle than the one reconstructed from longer signals, but it's still much worse than the image in lossy homogeneous medium.

Chao

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi Brad and Ben,

I made some more observations of the reconstruction in 2D lossy heterogeneous medium, in which a annulus was used with different speed of sound (3 mm/us) and density (2 g/cm^3) from other parts of the medium (sos= 1.5 mm/us, rho=1 g/cm^3). And the absorption parameters were set to be alpha_power = 1.5, alpha_coeff = 3.

If I used the equation of state including both absorption and dispersion, I get the following image:

If I used stokes equation that only describes the absorption, the reconstructed image looks like:

It seems the latter one looks better, but it's still not good enough. I don't know why I can not get good images in the lossy heterogeneous medium.

Another thing is that both of these two images are reconstructed by filtering the abs and disp parameters, and the image obtained without regularization of absorption parameters is just blank. Because the inverse crime was committed when reconstructing, I don't know what's wrong without regularization.

Thanks a lot!

Best,
Chao

Posted 13 years ago #
bencox

Developer
Joined: Mar '10
Posts: 292

Hi Chao,

Maybe this is right. If the success of the heterogeneous lossless reconstruction is dependent on recombining multiply reflected waves, then it will be much worse in the absorbing case as those waves will have travelled much further and been attenuated much more than the unreflected, direct, arrivals. When the signal has been attenuated to below the noise floor there is no way to get it back, so a blurry image is the best you can expect.

You could try a smaller difference between the sound speeds and densities, perhaps change them by 10% from the background to start with and see if it works well in that case. The contribution from multiply reflected waves will be less, so you should see a better image if my argument is correct.

Also, if you reduce the absorption (but not to zero) does it get better?

I assume you are you including the heterogeneity (the annulus) both in the simulation of the pressure data and in the time reveral imaging. Is that right?

Kind regards,

Ben

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi Dr. Cox,

Now even in the homogeneous case, if I did not regularize the absorption parameters (no filter), I can not get good images, it looks like:

Because I used inverse crime, I don't understand why I can not get good images if I did not filtered the absorption parameters, while I can get good images if I filtered it.

Also, it's right that I included the heterogeneity in both forward and time reversal simulations. I will do the other simulations you mentioned in your reply later, and get back to you.

Thanks a lot!

Best,
Chao

Posted 13 years ago #
bencox

Developer
Joined: Mar '10
Posts: 292

Hi Chao,

When simulating the data, the absorption reduces the amplitude of the various frequency components. For the highest frequencies this can reduce their amplitude to below the noise floor. When time reversing the data without using a filter (but with the absorption still turned on and not set to zero) the code amplifies each frequency component to try to correct for the reduction that occured during the simulation. Unfortunately, when the highest frequency components have been reduced so much they are in the noise, all that happens is that the noise is amplified, which is what you see in your image. It is therefore necessary to use a low-pass filter on the time-reversal step to prevent this from happening. It is a classic case of needing regularization for the inverse problem because there is a loss of information in the direct (forward) problem.

See our paper
Photoacoustic tomography in absorbing acoustic media using time reversal for more details.

Regards,

Ben

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi Dr. Cox,

Actually, I used inverse crime, which means there is no noise in the pressure data. I don't know if the filter is still necessary even in this case.

Thanks,
Chao

Posted 13 years ago #
bencox

Developer
Joined: Mar '10
Posts: 292

Chao,

There will always be some noise when running computations on a computer, even if only rounding errors. Even noise at a level as low as that could be enough to cause your problem when amplified exponentially during the time reversal.

Ben

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi Dr. Cox and Dr. Treeby,

Thanks for the explanation. It's very helpful, but I still have some questions about the absorption parameters defined in the k-wave code, because it looks different from the pseudo-code described in your paper (2010 inverse problem), which is shown as below:

However, in the code, the absorption parameters are defined as:

absorb_nabla1 = (kgrid.k.*sinc(c_max*dt*kgrid.k/2)).^(medium.alpha_power-2);
absorb_nabla1(isinf(absorb_nabla1)) = 0;
absorb_nabla1 = ifftshift(absorb_nabla1);
absorb_tau = -2*medium.alpha_coeff*c_max^(medium.alpha_power - 1);
absorb_param = absorb_tau.*absorb_nabla1.*rho0;

And the equation of state is:

p = c.^2.*( (rhox + rhoz) + real(ifft2( absorb_param.*(duxdx_k + duzdz_k) + dispers_param.*fft2(rhox + rhoz) )) );

My questions about the above code are:

1. Why 'absorb_tau', 'absorb_nabla1', and 'rho0' can be combined as one variable 'absorb_param' if they are separated by Fourier and inverse Fourier operators?

2. I am confused about which space the 'absorb_param' belongs to. From this line "absorb_nabla1 = ifftshift(absorb_nabla1)" and "absorb_param = absorb_tau.*absorb_nabla1.*rho0", it seems absorb_param belongs to spatial domain, while from 'ifft2( absorb_param.*(duxdx_k + duzdz_k) + dispers_param.*fft2(rhox + rhoz) ))', it should be in k-space domain.

3. Why to use c_max to compute absorb_tau, while

4. If my understanding is right, in the absorption term,

computes the derivative
,
and according to the equation of continuity
,
I think there should be a minus sign before

Posted 13 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

Hi Chao,

Thanks for your questions and for digging into the code to try and understand them. It seems you have discovered a well hidden bug. The problem is that the absorption terms were coded and validated based on homogeneous media, and the current extension to heterogeneous media in B.0.3 is not correct. To answer your questions in more detail:

1. In the homogeneous case, absorb_tau and rho0 are scalar variables and can be shifted between the Fourier transform integrals. This is not true in the heterogeneous case and thus the code needs some modification to be correct.

2. In the homogeneous case, absorb_param is a k-space variable (absorb_nabla is a k-space variable while absorb_tau and rho0 are scalars). In the heterogeneous case, these parameters will need to be used separately as you point out.

3. In the homogeneous case, c_max = c_0 so there is no difference. In the heterogeneous case, the maximum value was used so absorb_tau remained a scalar (this is consistent with the k-space sinc term in which the maximum sound speed is used to ensure stability).

4. The code is correct, however, you are right about a missing minus sign compared to the equations in the paper. The problem stems from a typographical error in the original JASA paper explaining the lossy operators (Eqs. 32 and 33). In these equations, the absorption and dispersion terms were given the wrong sign (Eq. 33 is also missing a set of brackets). There is an erratum for this paper being printed in JASA next month. Unfortunately, we copied the same equations into the Inverse Problems paper so these are also incorrect. In Eqs. 7 and 14, the two plus signs should be minus signs. Given the substitution you mention, there is also a missing minus sign in Eq. 15a.

I'll try and code up some bug fixes for the heterogeneous absorbing case in the coming days. Thanks for your persistence with trying to get to the bottom of this! Any more questions, please let us know.

Kind regards,

Brad.

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi Dr. Treeby,

Thank you for the clarifications, but I still have some questions:

1. I agree that absorb_nabla is a k-space veriable, in this case I think the line "absorb_nabla1 = ifftshift(absorb_nabla1)" is not necessary.

2. Is it required that absorb_tau needs to be a scalar? If not, I think c_max can be replace by c(x).

3. Another important aspect I want to consult you is about the k-space adjustment "sinc(c*k*dt/2)". If I understanding is right, in k-space method, the accuracy and efficiency of calculation of spatial derivative is improved by using Fourier transform, and the stability and accuracy of temporal derivative is improved by the k-space adjustment. Therefore, in my opinion, the k-space adjustment is not necessary in the line "absorb_nabla1 = (kgrid.k.*sinc(c_max*dt*kgrid.k/2)).^(medium.alpha_power-2)" since the equation of state does not involve temporal derivatives.

Based on the above 3 points, I modified the the corresponding parts of your delicate code as below:

1. Use "absorb_nabla1 = (kgrid.k).^(medium.alpha_power-2)" to replace "absorb_nabla1 = (kgrid.k.*sinc(c_max*dt*kgrid.k/2)).^(medium.alpha_power-2)"

2. Comment the line "absorb_nabla1 = ifftshift(absorb_nabla1)"

3. Use "absorb_tau = -2*medium.alpha_coeff*c.^(medium.alpha_power - 1)" to replace "absorb_tau = -2*medium.alpha_coeff*c_max.^(medium.alpha_power - 1)"

4. The equation of state becomes to " p = c.^2.*( (rhox + rhoz) + absorb_tau.*real(ifft2( absorb_nabla1.*fft2( rho0.*( real(ifft2(duxdx_k)) + real(ifft2(duzdz_k)) ) ) )) )" (Note: the dispersion term is not included here)

Based on these changes, I got the following reconstructed images using annulus heterogeneity:

It looks better than the image I got before (see the post above), and in homogeneous media, using the modified code will also give us better images as shown below (both exclude the dispersion term):
original:

modified:

However, if I include and modify the dispersion term in the same way, I get image like this for heterogeneous medium:

I don't know what's wrong with the dispersion term, and hope you give me some suggestions.

BTW, if the erratum for the JASA paper is being printed, could you send me a copy?

Thanks again!

Best,
Chao

Posted 13 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

Hi Chao,

Thanks again for your post. I am working on some bug fixes, but there are still a couple of issues I need to resolve. To answer your questions:

1. The fftshift and ifftshift functions just reorder the inputs between having the zero frequency component in the middle (the way the kgrid.k matrix is defined) to having the zero frequency component in the top left (the way MATLAB returns the output of the FFT). In the code, instead of always reordering the output of the FFT operations, we just reorder any k-space variables before we start. So this line of code is still needed.

2. If you use absorb_tau explicitly, there is no reason it needs to be a scalar.

3. I agree with your general line of thinking, however, in practice it seems to make very little difference. I think our reasoning for leaving the k-space operator in was that the accuracy of the absorption operators seemed to be slightly improved. I agree that in either case a more rigorous argument is required.

My only other comment is to be careful of basing your measure of "correctness" of the absorption terms on the accuracy of an inverse problem example. It might be better to use a simple forward example in which you extract the absorption and dispersion in a homogeneous region of a heterogeneous medium (Similar to the Modelling Power Law Absorption Example).

You can grab the erratum here.

Kind regards,

Brad.

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi, Dr. Treeby,

Thank you for your suggestion, and now I have tested my modified code in a 1D heterogeneous medium (sound speed, density and absorption coefficient (alpha_0) are all heterogeneous). And I compared the simulated absorption and dispersion with the ones predicted by frequency power law and KK relation. The codes I use are as following:

Nx = 512;
dx = 0.1e-3;
kgrid = makeGrid(Nx, dx);

medium.sound_speed = 1500*ones(Nx,1);
medium.sound_speed(200:300) = 3000;
medium.density = 1000*ones(Nx,1);
medium.density(200:300) = 2000;
medium.alpha_coeff = zeros(Nx,1);
medium.alpha_power = 1.5;

[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed,0.05);
kgrid.t_array = (0:13999)*dt;

source.p0 = zeros(Nx, 1);
source.p0(50) = 1;

sensor.mask = zeros(Nx, 1);
sensor.mask(100) = 1;
sensor.mask(400) = 1;

pw = my_kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotSim', false);

medium.alpha_coeff(200:300) = 1;

ps = my_kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotSim', false);

[f, asw, psw] = spectrum(pw(10501:13500), 1/dt);
[f, ass, pss] = spectrum(ps(10501:13500), 1/dt);

pf=log(asw(1:25)*ass(13)./(ass(1:25)*asw(13)))*100;

beta = db2neper(1, 1.5)*(2*pi*10^6)^1.5;

f=(0:24)*0.2; % MHz

figure;plot(pf,'o');hold on;plot(beta*(f.^1.5-2.4^1.5))

uwpsw=unwrap(psw);
uwpss=unwrap(pss);

qf=(uwpss(2:25)-uwpsw(2:25))./(2*pi*f(2:25)*0.01)-(uwpss(13)-uwpsw(13))./(2*pi*f(13)*0.01);

figure;plot(qf,'o');hold on;plot(-beta/2/pi*tan(1.5*pi/2)*(f(2:25).^0.5-2.4^0.5))

Here I used the method proposed by Ping He (you can grap it here: goo.gl/zrPDx) to calculate absorption and dispersion function 'pf' and 'qf', which are used to compare with the resutls computed from power law and KK relation. And the results are shown below:

Best,
Chao

Posted 13 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

Hi Chao,

Thanks for your post, we certainly appreciate your help in identifying and resolving this bug. It looks like you have had some success with your modifications. We plan to revert to explicitly using the absorption equations as given in the inverse problems paper (without the k-space adjustment) for the next release. This should be available by the middle of September.

Thanks again,

Brad.

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

Hi Dr. Treeby,

Actually I also explicitly use the absorption equations without k-space adjustment, but I notice that Dr. Cox mentioned in another thread in this forum that, quote, "absorption is included in the wavenumber domain (because it prevents us having to store the field at previous time-steps and so is efficient) which necessarily implies a homogeneous absorption value everywhere".

I don't understand why the absorption value has to be constant if the absorption is implemented in the wavenumber domain.

Thanks,
Chao

Posted 13 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

Hi Chao,

This comment is in reference to the power law exponent y which is used in the two fractional Laplacian operators. Because these operators are applied in the spatial frequency domain using the Fourier collocation spectral method, y must necessarily be homogeneous in x-space. If you can think of a way around this restriction, please let us know!

Kind regards,

Brad.

Posted 13 years ago #
huangchao

Member
Joined: Nov '10
Posts: 57

I see. I thought the absorption value means the power law pre-factor alpha_0, that's why I am confused. Actually, I didn't change the power law exponent y. I just change alpha_0 to be varying in x-space in my modified code.

Best,
Chao

Posted 13 years ago #