Thread Subject: Find gaussian humps in a 3D dataset

Hi all,

I have a 3D image (scalar intensity field) of particles in a fluid flow.

There may be many particles (1<N<1000), in a large domain (up to 400^3 voxels)

Each particle could be approximately represented by a gaussian hump in intensity, surrounding the centre position. Standard deviation of the gaussian might be ~ 2 voxels.

So the question is, how best to find the positions of the particles?

I suspect that an optimisation would be very computationally intensive. I can find local maxima very quickly to the accuracy of a single voxel; but I need to get sub-voxel accuracy - perhaps using some kind of polynomial fitting then a derivative-based approach.

Has anyone got any ideas?

Thanks,

Tom Clark

"Thomas Clark" <t.clarkremove_spam@cantab.net> wrote in message <galg9a$995$1@fred.mathworks.com>...
> Hi all,
>
> I have a 3D image (scalar intensity field) of particles in a fluid flow.
>
> There may be many particles (1<N<1000), in a large domain (up to 400^3 voxels)
>
> Each particle could be approximately represented by a gaussian hump in intensity, surrounding the centre position. Standard deviation of the gaussian might be ~ 2 voxels.
>
> So the question is, how best to find the positions of the particles?
>
> I suspect that an optimisation would be very computationally intensive. I can find local maxima very quickly to the accuracy of a single voxel; but I need to get sub-voxel accuracy - perhaps using some kind of polynomial fitting then a derivative-based approach.

No. Don't bother with a polynomial model.

As long as you are willing to assume there
is no near overlap between humps, find the
local peaks, then grab about 5 pixels in each
direction and fit that with a single gaussian
mode. Go on to the next local peak.

John

Thomas Clark wrote:
> Hi all,
>
> I have a 3D image (scalar intensity field) of particles in a fluid flow.
>
> There may be many particles (1<N<1000), in a large domain (up to 400^3 voxels)
>
> Each particle could be approximately represented by a gaussian hump in intensity, surrounding the centre position. Standard deviation of the gaussian might be ~ 2 voxels.
>
> So the question is, how best to find the positions of the particles?

Tom, it's possible that the GMDISTRIBUTION function in the Statistics Toolbox, which fits a Gaussian mixture distribution, might be helpful. I don't know enough about your application, though.

- Peter Perkins
  The MathWorks, Inc.

Thanks to you both.

John, I did have some fear that the polynomial approach wouldn't work... so thanks for confirming that. It's a shame, because 'windowing' in such a large number of particles could be computationally expensive. We'll see!

I'll implement both solutions and see which works best. I'll write back to let you know (might be a while, as I have yet to get any experimental data!).

Kind Regards

Tom Clark

On Sep 15, 10:45=A0am, "Thomas Clark" <t.clarkremove_s...@cantab.net>
wrote:
> Thanks to you both.
>
> John, I did have some fear that the polynomial approach wouldn't work... =
so thanks for confirming that. It's a shame, because 'windowing' in such a =
large number of particles could be computationally expensive. We'll see!
>
> I'll implement both solutions and see which works best. I'll write back t=
o let you know (might be a while, as I have yet to get any experimental dat=
a!).

I suggest you start with similated data. This is not an easy
problem because of the difficulty of separating overlapping Gaussians.

Hope this helps.

Greg

Greg Heath <heath@alumni.brown.edu> wrote:
>
> I suggest you start with similated data. This is not an easy
> problem because of the difficulty of separating overlapping Gaussians.
>

That certainly seems like a good idea to me. It might take quite a bit of
practice to tune your methods here.

An alternative method of detection might be to first normalise your data,
and then attempt to spot deviations from the norm.

Take your data set and apply a multi-dimensional smoothing to it, like a
moving average, but in higher dimensions. You can do this quite
efficiently using convolutions, or equivalently using fftn().

Once you have the smoothed function, subtract it from your original data
set, and you will be left with the deviations from the local average
levels. Now you can search the whole data set for places where the
deviation from the average is above your threshold level.

As Greg said, though, the overlapping Gaussians might cause you some
extra troubles. You may be able to locate these by noticing greater than
expected intensity, indicating there is more than one particle. I'm not
sure what the best method of fitting would be, though. It does sound like
a constrained optimisation problem.

You said that you have maybe 1000 particles. I take it that you also have
many data sets, not just the one?

--
Dr Tristram J. Scott
Energy Consultant

Greg,

Thanks, yes, I'll be starting with simulated data (gaussians plus noise).

One fortunate feature is that particles don't tend to clump together, and of course they can't be in the same place... so any two peaks will occur at least two standard deviations apart. Hopefully that'll help re. overlapping peaks.

Tristram,

Thanks for the additional idea, I'd not come across that before. I'll store it away for the future, but not implement that yet... as I see it, I'd still have to make some kind of fit to detect the exact location of the peak (thresholding could only ever give voxel-scale accuracy, rather than sub-voxel).

And yes, ~1000 particles per frame, ~1000 frames per test, ~100 tests to get a PhD. It'll need to be relatively quick :( !

Kind Regards

Tom Clark



tristram.scott@ntlworld.com (Tristram Scott) wrote in message <fwJzk.60139$6Y7.42807@newsfe29.ams2>...
> Greg Heath <heath@alumni.brown.edu> wrote:
> >
> > I suggest you start with similated data. This is not an easy
> > problem because of the difficulty of separating overlapping Gaussians.
> >
>
> That certainly seems like a good idea to me. It might take quite a bit of
> practice to tune your methods here.
>
> An alternative method of detection might be to first normalise your data,
> and then attempt to spot deviations from the norm.
>
> Take your data set and apply a multi-dimensional smoothing to it, like a
> moving average, but in higher dimensions. You can do this quite
> efficiently using convolutions, or equivalently using fftn().
>
> Once you have the smoothed function, subtract it from your original data
> set, and you will be left with the deviations from the local average
> levels. Now you can search the whole data set for places where the
> deviation from the average is above your threshold level.
>
> As Greg said, though, the overlapping Gaussians might cause you some
> extra troubles. You may be able to locate these by noticing greater than
> expected intensity, indicating there is more than one particle. I'm not
> sure what the best method of fitting would be, though. It does sound like
> a constrained optimisation problem.
>
> You said that you have maybe 1000 particles. I take it that you also have
> many data sets, not just the one?
>
> --
> Dr Tristram J. Scott
> Energy Consultant

Thomas Clark <t.clarkremove_spam@cantab.net> wrote:
> One fortunate feature is that particles don't tend to clump together, and
> of course they can't be in the same place... so any two peaks will occur at
> least two standard deviations apart. Hopefully that'll help re. overlapping
> peaks.

That is very useful.

>
> Thanks for the additional idea, I'd not come across that before. I'll
> store it away for the future, but not implement that yet... as I see it,
> I'd still have to make some kind of fit to detect the exact location of the
> peak (thresholding could only ever give voxel-scale accuracy, rather than
> sub-voxel).

I am assuming that your intensity data can be interpreted as a probability
density function for the Gaussian. Assuming that is the case, then your
discrete approximation to the PDF should be usable to find the mean by
looking at the moments of the PDF.

mean = \int (x.p) .dx

Presumably for a symmetric distribution of intensities you can do this in
each dimensions independently of the others?


> And yes, ~1000 particles per frame, ~1000 frames per test, ~100 tests to
> get a PhD. It'll need to be relatively quick :( !
>

Perhaps it is not so much that the process needs to be fast, but rather
that it needs to operate without manual intervention.

Drop me an email if you want more details on any of my vague ideas.

--
Dr Tristram J. Scott
Energy Consultant

Peter Perkins <Peter.PerkinsRemoveThis@mathworks.com> wrote in message <galiqa$6es$1@fred.mathworks.com>...
> Thomas Clark wrote:.....
> > I have a 3D image (scalar intensity field) of particles in a fluid flow. .....
> > Each particle could be approximately represented by a gaussian hump in intensity, surrounding the centre position. .....
> > So the question is, how best to find the positions of the particles?
>
> Tom, it's possible that the GMDISTRIBUTION function in the Statistics Toolbox, which fits a Gaussian mixture distribution, might be helpful. I don't know enough about your application, though.
>
> - Peter Perkins
> The MathWorks, Inc.
Hello,
    I have a question similar to one Thomas Clark asked (above): I have 3D scalar (intensity) field and I have to find particle centers by fitting mixed Gaussian distribution.
    gmdistribution.fit you suggested seems like good idea except for one problem: the input data are “events” or “hits” in XYZ space, which have no value. To convert each Xi,Yi,Zi voxel with intensity Ni to this form, I should add Ni identical [Xi Yi Zi ] rows to the input, to result in a huge matrix. Is there a better workaround or maybe other function similar to gmdistribution.fit but taking scalar data field as input?
  

Valentin Magidson wrote:

> I have a question similar to one Thomas Clark asked (above): I have 3D scalar (intensity) field and I have to find particle centers by fitting mixed Gaussian distribution.

Valentin, it sounds like you're fitting a (hyper)surface in the sense of a regression model, rather than fitting a probability density. The latter is what gmdistribution is intended for, not the former. You may want to look into using nlinfit instead, with a model function that is the sum of several Gaussians, but Gaussian distributions. You may find

<http://www.mathworks.com/products/statistics/demos.html?file=/products/demos/shipping/stats/cfitdfitdemo.html>

helpful in that distinction, though it discusses ting in only one dimension. You will likely need to initialize nlinfit with some reasonable starting values for the components. Such models are notorious for having many bad local minima of the sum of squares surface.

Hope this helps.

- Peter Perkins
  The MathWorks, Inc.