As i've explained in a previous question: I have a dataset consisting of a large semi-random collection of points in three dimensional euclidian space. In this collection of points, i am trying to find the point that is closest to the area with the highest density of points.
As high performance mark answered;
the most straightforward thing to do would be to divide your subset of Euclidean space into lots of little unit volumes (voxels) and count how many points there are in each one. The voxel with the most points is where the density of points is at its highest. Perhaps initially dividing your space into 2 x 2 x 2 voxels, then choosing the voxel with most points and sub-dividing that in turn until your criteria are satisfied.
Mark suggested i use triplequad
for this, but this is not a function i am familiar with, or understand very well. Does anyone have any pointers on how i could go about using this function in Matlab for what i am trying to do?
For example, say i have a random normally distributed matrix A = randn([300,300,300])
, how could i use triplequad
to find the point i am looking for? Because as i understand currently, i also have to provide triplequad
with a function fun
when using it. Which function should that be for this problem?
Here's an answer which doesn't use triplequad
.
For the purposes of exposition I define an array of data like this:
A = rand([30,3])*10;
which gives me 30 points uniformly distributed in the box (0:10,0:10,0:10)
. Note that in this explanation a point in 3D space is represented by each row in A
. Now define a 3D array for the counts of points in each voxel:
counts = zeros(10,10,10)
Here I've chosen to have a 10x10x10 array of voxels, but this is just for convenience, it would be only a little more difficult to have chosen some other number of voxels in each dimension, and there don't have to be the same number of voxels along each axis. Then the code
for ix = 1:size(A,1)
counts(ceil(A(ix,1)),ceil(A(ix,2)),ceil(A(ix,3))) = counts(ceil(A(ix,1)),ceil(A(ix,2)),ceil(A(ix,3)))+1
end
will count up the number of points in each of the voxels in counts
.
EDIT
Unfortunately I have to do some work this afternoon and won't be able to get back to wrestling with the triplequad
solution until later. Hope this is OK in the meantime.