Hi,
I agree with almost, but not quite, everything here.
Certainly it is all about step 3, and that you can write things
equivalently as
matrices (no infinite support needed).
The one thing I don't see is why you are trying to define noise in the
image.
The noise is what happens prior to step 3, and once we've got our image
reconstructed, then everything is locked into that image. It is this
image that
is being convolved and this image that we want to recover with some
deconvolution process. Hence we want the exact same noise back as well.
That is why I was talking about rounding errors and machine precision,
as
this is where all the errors come from. With infinite precision, it
would
definitely be possible to deconvolve the image and get back the image
that you started with (which contained the locked-in noise).
Also, thinking about it a little more, for a Gaussian with a small
FWHM in
image space, the spectral representation actually does not get small
over the range of values needed (easily verified in matlab). So the
problem is actually reasonably well-posed. What becomes difficult is
(a) knowing *exactly* what the correct convolution kernel was, as any
mis-specification will result in nasty artifacts after the
deconvolution,
and (b) whether something non-trivial was done at the edges. Certainly
something will have been done at the edges, as the simple FFT-style
Gaussian convolution is not used in practice as it wraps around the
FOV due to the cyclic nature of the discrete Fourier transform. So it
may be implemented as a zero-padded image with FFT, where the
outside padding is thrown away once the convolution is done (which
is problematic for trying to invert this process). Alternatively, it
may
be renormalised by dividing by a function to take into account the
edge effects (e.g. the convolved, padded mask), or have a special
computation applied at the edges. If any of this is non-linear then
all bets are off. However, even if it is linear you generally end up in
a situation where you cannot represent this with FFTs and need to
use the matrix formulation. This then requires calculating a very
precise inverse of a matrix which is typically extremely poorly
conditioned, mainly due to the edge effects. If information honestly
has been lost (e.g. throwing away padded regions) then the matrix
will be rank-deficient and it is not possible to do the deconvolution
at all.
So, the bottom line isn't really different - I doubt you will be able to
perform an accurate and useful deconvolution. However, I do not
think it has anything to do with the SNR of the original image that is
formed (since the noise is locked in) but much more to do with the
accuracy that you can precisely specify and model the smoothing
process, and the ability to calculate a precise inverse given known
algorithms and machine-precision.
All the best,
Mark
On 7 Mar 2009, at 19:43, Souheil Inati wrote:
> 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
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>
>
|
