Hi all,

Okay, I'm jumping in late, but here's my point of view. I agree with
Jeff and Bruce. In general "deblurring" is a bad idea because the
problem is extremely ill-conditioned. But there are some cases where
it might be possible, and Mark's point is correct. Please excuse the
lengthy post, but I feel the need to be careful :-)

Let's restrict the conversation to smoothing of a previously acquired
rectilinear EPI image, i.e. the following sequence of steps occur:
1) Data are sampled by the scanner in time. The scanner performs an
analog operation - continuous integral of the instantaneous
magnetization across space. The signal is digitized at discrete time
points.

2) An image is reconstructed. Time is converted to k. A k-space
filter is applied (this is optional but is often applied to reduce
Gibbs in low-res images). 1-D recon in the frequency direction
(interpolation and FFT). Navigator correction is applied. FFT is
applied in the phase direction. Signal is combined across coils as
sqrt of sum of squares, or if 1 coil, absolute value.

[Aside: Okay, after all that physicist junk we have what Mark likes to
work with, namely an image :-) ]

3) Now we smooth the image, say with a gaussian. This can be done
either by convolution in the image domain or FFT, multiply by a
frequency domain filter, IFFT.

As far as I see it, the question is only about step 3. Can step 3 be
inverted if the convolution kernel is known. This is a VERY specific
question. I mean, it is related to, but is not exactly the same
problem as deblurring a Hubble telescope image say. And, importantly,
this has absolutely nothing to do with k-space in the MRI sense. You
had already lost the phase information, etc.. You took a noisy image
with some weird artifacts (distortion, blurring, yuck.) and blurred
it. That is all.

I want to be very precise about the statement "invert step 3". If
it's well conditioned, I want the inverse. If it's badly conditioned,
I want a pseudo-inverse type least squares solution, preferably in a
way that let's me control the noise amplification. I will work with
the (FFT, multiply, IFFT) description of convolution because I can
write the problem down in a clean way. Here's what we have (think
matlab):
    I_s = Finv* D * F * I_orig
where I_orig is the original image as a column vector, F is the
discrete fourier transform as a matrix, D is a diagonal matrix
(multiply each point in the frequency domain by a number), Finv is the
discrete inverse fourier transform.

I want to do this:
   I_est = Finv * Dinv * F * I_s
so that I_est is as close to I_orig as possible.

On the computer, you should work in double precision so that you don't
have to worry too much about rounding errors. In double precision on
the computer Finv*F is the identity to a whole lot of digits.

What about Dinv? D is diagonal, so I know how to invert it. If d_j
is the j-th element on the diagonal, the dinv_j = 1/d_j. This is bad
of course, because some of the d_j's are very small - i.e. out at high
spatial frequency, in the tail of the gaussian filter, we had
multiplied by a very small number. But we can do a controlled
inverse, either a pseudo-inverse i.e. cutting off at some singular
value, i.e. dinv_j = 0 for d < cutoff, or we can do Tikhonov
regularization, i.e. dinv_j = (dinv_j) ./ (dinv_j^2 + sigma^2).

The question is what is the structure of the noise in I_s, and how do
we pick a cutoff for the pseudo-inverse. Other than the noise
inherent in the scanner, there is the stupidity of storing images in
DICOM and restricting them to 12-bit integers, etc., etc., so, I
would argue that at best, you have three digits of precision in your
original image, i.e. your noise is of the order of a few percent of
your signal. So you should cut off (or roll off) the inverse
somewhere around there, say 0.1.

So here's a practical procedure
I_est = Finv * Dinv * F * I_s
Use FFT for F and IFFT for Finv. For Dinv, build a filter that looks
like this:
d(f) = exp(-a f.^2) / (exp(-2 a f.^2) + 0.01)
where `a` is the original filter width.

This will give you an image I_est which is close to I_orig. Not
exactly the same, but as close as you can get in a least squares sense
and not have the noise explode.

Does this make sense?

Cheers,
Souheil

PS - Of course, you should try this out with new data first. Go
through the smoothing and the "unsmoothing" and see what happens.


---------------------------------

Souheil Inati, PhD
Research Associate Professor
Center for Neural Science and Department of Psychology
Chief Physicist, NYU Center for Brain Imaging
New York University
4 Washington Place, Room 809
New York, N.Y., 10003-6621
Office: (212) 998-3741
Cell: (646) 522-7330
Fax: (212) 995-4011
Email: [log in to unmask]



On Mar 7, 2009, at 1:08 PM, Mark Jenkinson wrote:

> Hi,
>
> My understanding of this is that you are essentially working with a
> discrete
> signal and applying a discrete operation to it. You can either
> consider this
> as a discrete convolution operation in image space, or as a
> multiplication in
> the Fourier space, but not the continuous k-space, the space you get
> to with
> the discrete FFT. This does contain phase and magnitude information,
> but not necessarily the same as the original k-space, as they are
> not the
> same due to the magnitude operation and the discrete operations.
> However, it is perfectly equivalent to perform the convolution by
> using the
> FFT, multiplying and then doing the inverse FFT. Any description in
> terms
> of diffusion, or k-space are actually continuous analogues of this
> process
> and not quite what you are really doing by blurring the discrete,
> magnitude
> reconstructed signal.
>
> As for inverting the operation, the problem is simply one of machine
> precision. Since we are talking about an already acquired signal that
> has been convolved with a (discrete) Gaussian, then even noise is not
> an issue. So what really counts is the ability to represent the
> values
> sufficiently accurately, and the precision to which the FFT can be
> calculated. Because the FFT involves a large summation, much of
> which is partially canceling, it is sensitive enough that the
> inversion is
> usually not sufficiently accurate since the suppression and then
> enhancement of the high-frequency components in the image
> typically have rounding errors which affect the image sufficiently to
> cause problems.
>
> So I agree that in practice it probably will not give good enough
> results,
> but in theory there isn't really a loss of information in the
> discrete case
> except with respect to the machine precision.
>
> All the best,
> Mark
>
>
>
> On 7 Mar 2009, at 16:43, Kochunov, Peter wrote:
>
>> Thank you Bruce,
>> I think, your explanation is quite sensible. Obviously, Gaussian
>> filtering in the magnitude domain is not equivalent to the
>> frequency-space multiplication as half of the information (phase)
>> is lost. Not to mention that the fourier transform of the inverse-
>> gaussian (Wald) function is not analytically defined.
>> pk
>>
>> ________________________________
>>
>> From: FSL - FMRIB's Software Library on behalf of Bruce Fischl
>> Sent: Sat 3/7/2009 10:20 AM
>> To: [log in to unmask]
>> Subject: Re: [FSL] Actual implementation? [Re: Q: How to de-smooth
>> BOLD images, previously smoothed with a known kernel-width?]
>>
>>
>>
>> Hi Peter,
>>
>> I'm pretty sure it is. You can think about Gaussian convolution as
>> moving
>> intensities around the image since the diffusion equation obeys an
>> underlying conservation law. So any one voxel has its intensity
>> "come from"
>> other voxels in the image. Unfortunately there is nothing unique
>> about it,
>> so you don't know for example if the intensity came from a high
>> value voxel
>> far away or a moderately valued voxel nearby. The inverse diffusion
>> equation is fundamentally ill-conditioned. The k-space thing is
>> probably
>> only true in the limit of infinite support, etc..., but I haven't
>> thought
>> about it.
>>
>> cheers,
>> Bruce
>>
>>
>> On Sat, 7 Mar 2009, Kochunov, Peter wrote:
>>
>>> Bruce,
>>> Is that really the case? I mean, the k-space-domain operations
>>> equivalent
>> to convolution/deconvolution with the Gaussian function are
>> inversable?
>>> pk
>>>
>>>
>>> -----Original Message-----
>>> From: FSL - FMRIB's Software Library on behalf of Bruce Fischl
>>> Sent: Sat 3/7/2009 8:43 AM
>>> To: [log in to unmask]
>>> Subject: Re: [FSL] Actual implementation? [Re: Q: How to de-smooth
>>> BOLD images, previously smoothed with a known kernel-width?]
>>>
>>> Hi Raj,
>>>
>>> Gaussian blurring is the equivalent of running the diffusion
>>> equation for
>>> time proportional to sigma^2 (since the Gaussian is the Green's
>>> Function of
>>> it), which is not time-reversible. Information is irretrievably
>>> lost in
>>> diffusion, so I'm afraid the inversion isn't possible.
>>>
>>> sorry :<
>>>
>>> Bruce
>>>
>>> On Fri, 6 Mar 2009, Rajeev Raizada wrote:
>>>
>>>> On Fri, 6 Mar 2009 09:27:24 -0800, Michael T Rubens
>>>> <[log in to unmask]> wrote:
>>>>
>>>>> take FFT of smoothed image, divided by FFT of gaussian. the
>>>>> inverse FFT
>>>>> should be your unsmoothed data.
>>>>
>>>> Thanks...
>>>> But please see below... :-)
>>>>
>>>>> On Fri, Mar 6, 2009 at 5:12 AM, Rajeev Raizada <[log in to unmask]
>>>>> wrote:
>>>> [...]
>>>>>> Non-specific high-level exhortations to recast the smoothing
>>>>>> as a 3D Fourier filter and then to apply the inverse filter
>>>>>> are also welcome, but probably won't be quite as useful :-)
>>>>
>>>> I believe that the application of an inverse filter
>>>> may be easier said than done.
>>>> It appears that for Gaussian deblurring, the inverse is "ill-
>>>> conditioned",
>>>> e.g. http://ieeexplore.ieee.org/iel5/5992/26914/01196312.pdf
>>>>
>>>> Two additional complications:
>>>> 1. Apparently there are some analytical results for deblurring of
>>>> 2D discrete Gaussians,
>>>> but I don't know enough to know whether these hold in 3D as well.
>>>> 2. I believe that the 3D smoothing is actually done by a Gaussian
>>>> convolved
>>>> by a sinc function, not just a plain vanilla Gaussian.
>>>>
>>>> Does anyone have an actual implementation of such "de-smoothing",
>>>> as opposed to an "in principle" description of what it ought to
>>>> involve?
>>>> Googling for gaussian deblurring turns up a lot of hits for blind
>>>> deconvolution
>>>> and methods of counteracting noise.
>>>> However, in this case the deconvolution is not blind at all,
>>>> as we know that it was a gaussian kernel of FWHM 6mm,
>>>> and also there wasn't any noise in the blurring process.
>>>> So, in principle those two facts ought to make things easier, I
>>>> think?
>>>>
>>>> Any help greatly appreciated.
>>>> The more specific the better. :-)
>>>>
>>>> Raj
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>