Fitting a curve to 3D data

Hi people, I have a problem with fitting 3D data. I have a 3D matrix (99x386x384) and each dimension of the matrix represents x, y and z coordinates and value inside the matrix s value at a certain point value(x,y,z). Basically, I have a cluster of data points. I want to fit that data to the 3D curve/volume. Does anyone know how to do it?

[KSSV]

Why you want to fit data to 3D volume?

[Bjorn Gustavsson]

So you have a scalar s that varies in 3-D? And you want some specific advice how to best represent that volume-distribution with a simple "few-parameter" function? What does the function look like? Is it something like a 3-D Gaussian? Smooth linear increase in an arbitrary gradient direction? Does s obey some physics-based conservation-laws, like a convection-diffusion dominated continuity-equation?

So many question - and I havent even started yet...

[Nikola Segedin]

@KSSV I'm a physicist who works in the hospital at Gamma knife (we treat bran tumours using ionizing radiation) and I need to obtain an analytica expression for a single shot on a gamma knife so I can work with it. A single shot ai approx spherical distribution where each point represents a value of the delivered dose to the tumour volume.

@Bjorn Gustavsson Yes, it is a scalar field. A function can be anything from Gauss to whatever, as long as it fits the data. The best fit I got for this type of distribution in 1D was a polynomial of 16th order. The distribution in 1D looks like Gauss but it is broader at the top so Gauss is not the best fit. In 3D it is approx. a sphere and the value of the point falls as you move from the centre of the distribution. There is no conservation law that this data would follow. It is just experimental data, obtained at the lab, we know that the distribution is the sphere in 3D, and in order to manipulate the data we need to obtain an analytical expression.

[Nikola Segedin]

@Bjorn Gustavsson An analytical expression would be a function in form f(x,y,z)

[Walter Roberson]

The best technique is going to depend a lot on the mathematical model that you are using.

You indicate that you have a cluster of data points. Is that along the lines that most of the locations in the array have an s value indicating "not present" and that the other ones form a shape that you want to model along with the s value? For example s==0 being unoccupied and all other values being meaningful?

Or is it a cuboid in which all s values are relevant?

[Nikola Segedin]

Most of the data points in the array have a value. Others are indicated as 0 and are not relevant. I can just remove the if necessary.

[Walter Roberson]

Degree 16 polynomials are almost always garbage for coordinates outside +/- 1. The high order terms grow so quickly that they overwhelm numeric accuracy considerations.

It would be more likely that a trig or exponential or Gaussian or spherical bessel was involved, with what you found being similar to a truncated Taylor expansion. (A truncated Taylor expansion of a trig function gets very bad by half a period but the original trig can still be quite meaningful.

[Nikola Segedin]

I can play later wit the fit function, for now I need an idea how to even start!

[Torsten]

I wonder why people always want analytical expressions, especially for something like a surface in 4d.

Use

sq = interp3(X,Y,Z,S,xq,yq,zq)

if you want a good approximation sq to your data in a query point (xq,yq,zq).

[Bjorn Gustavsson]

@Nikola Segedin if it looks like a Gaussian then tweak the Gaussian instead of jumping to polynomial fits. Perhaps something like this would be better:

Or some similar modifications...

Answer 1

[Bjorn Gustavsson]

If you need to "get an idea to start" you can start from here:

function err = your_3D_error_fcn(pars,x,y,z,S,sigmaS,idx_is_OK_linear,fit_fcn)

S_model = fit_fcn(pars,x,y,z);
err = sum((S(idx_is_OK_linear)-S_model(idx_is_OK_linear)).^2./sigmaS(idx_is_OK_linear).^2)
end

This function you could use to find the best fiting parameters for a model-function using fminsearch:

fit_G3D = @(pars,x,y,z) pars(1)*exp(-(x-pars(2)).^2/pars(3)^2).* ...
                                exp(-(x-pars(4)).^2/pars(5)^2).* ...
                                exp(-(x-pars(6)).^2/pars(7)^2);
% Guessing parameters for an initial 3-D Gaussian with peak at 1 centred at
% [x,y,z] = [2 3 4] with widths in all directions equal to 1/2
pars0 = [1, 2,1/2,3,1/2,4,1/2];
pars_best = fminsearch(@(pars) your_3D_error_fcn(pars,x,y,z,S,sigmaS,idx_is_OK_linear,fit_G3D),pars0);

Here you'll have to provide the 3-D coordinates of x, y, z and s as well as an array with the linear indices to the good points and the standard-deviation of s (if that doesn't apply just remove it from the error-function). The error-function you'll have to adjust to something that suits your problem.

HTH

[Nikola Segedin]

Thanks, I will try it.