How to add two different surface curves in a single plot?

I have a bell shaped curve in 3D (generated from curve fitter app) as shown in the figures.

Both are the same curve the only difference is one is without point data and the other has point data.

Problem: I want to add another surface to the same plot for the same data which is at z=1 and parallel to x- and y-axis. Also, is it possible to find the values of major and minor axes of the ellipse formed from the intersection of both the surfaces. Or in other words, the values between the intersection of x and z values of both the curves and the y and z values of both the curves.

I hope I am able to express myself clearly.

Thank you

Answer 1

[Star Strider]

Try this —

‘I want to add another surface to the same plot for the same data which is at z=1 and parallel to x- and y-axis.’

[X,Y] = ndgrid(-3:0.1:3);
f  = @(x,y) exp(-(x.^2+(2*y).^2)*0.5);
Z = f(X,Y)*3;

figure
surf(X, Y, Z)
hold on
surf(X, Y, ones(size(Z)), 'FaceColor','r', 'FaceAlpha',0.5, 'EdgeColor','none')
hold off
colormap(turbo)

‘Also, is it possible to find the values of major and minor axes of the ellipse formed from the intersection of both the surfaces.’

figure
[c,h] = contour(X, Y, Z, [1 1]);
axis('equal')
grid
Ax = gca;
Ax.XAxisLocation = 'origin';
Ax.YAxisLocation = 'origin';

elpsfcn = @(b,xy) xy(1,:).^2/b(1)^2 + xy(2,:).^2/b(2)^2 - b(3);
opts = optimoptions('fminunc', 'MaxFunctionEvaluations', 5E+3, 'MaxIterations',1E+4);
[B, fv] = fminunc(@(b) norm(elpsfcn(b,c(:,2:end))), rand(3,1), opts)

Local minimum possible.

fminunc stopped because it cannot decrease the objective function
along the current search direction.

B = 3×1

52.8204 26.4510 0.0008

fv = 1.2653e-05

fprintf('Semimajor Axis = %.4f\nSemiminor Axis = %.4f\nConstant       = %.4f\n', B)

Semimajor Axis = 52.8204
Semiminor Axis = 26.4510
Constant       = 0.0008

text(-2.5, 2.5, sprintf('$\\frac{x^2}{%.2f^2} + \\frac{y^2}{%.2f^2} = %.4f$',B), 'Interpreter','latex', 'FontSize',16)

.

[Deepshikha Deo]

Thanks a lot for taking your time and answering the question with an example. I tried the sample code and it works well but I am having the trouble with replacing it with my data, I hope you can help me.

I have used curve fitter app and exported the output to the workspace. Thus a sfit of name 'fittedmodel' has been created in the workspace.

I replaced sixth line of your code (which is surf(X, Y, Z) ) to plot(fittedmodel,[x,y],z)

This created the following image

After that I followed the code as it is and it showed an error that stated:

Error using surf

Z must be a matrix, not a scalar or vector.

In an another trial, I added the function in attempts to create a matrix as follows:

f=@(x,y) a0*exp(-(((x-x0).*cos(phi)+(y-y0).*sin(phi)).^2/0.5*ax^2)-(((y-y0).*cos(phi)-(x-x0).*sin(phi)).^2/0.5*ay^2))

z=f(X,Y)*3;

But it also showed an error stating:

Unrecognized function or variable 'a0'.

Error in @(x,y)a0*exp(-(((x-x0).*cos(phi)+(y-y0).*sin(phi)).^2/0.5*ax^2)-(((y-y0).*cos(phi)-(x-x0).*sin(phi)).^2/0.5*ay^2))

I request you to help me with the problem. I have also attached the zip file which contains the data in excel and saved sfit file.

Thank you

Regards

[Star Strider]

I finally figured out how to open the ‘.sfit’ file, however getting information from it is beyond hopeless. It would be better to provide the code used to create it, since that may be an option. Anyway, I just gave up on getting the information from it, since that requires spelunking through multiple levels of the structure it creates.

(I don’t have the Curve Fitting Toolbox since I don’t need it for what I do, so I have very little experiece with it. The Statistics etc,, Optimization and Global Optimization Toolboxes do everything I need.)

It would be best to post the code used to create it, since I can work with that.

Please explain what the variables are in:

f=@(x,y) a0*exp(-(((x-x0).*cos(phi)+(y-y0).*sin(phi)).^2/0.5*ax^2)-(((y-y0).*cos(phi)-(x-x0).*sin(phi)).^2/0.5*ay^2))

Which are the parameters and which are the data?

How do I use the data in the .xlsx file? Specifically, what is ‘grid_code’ and does it have any significance in this problem?

In the interim, this works as far as I can take it —

Uz = unzip('BG_Apr03.zip')

Uz = 1×2 cell array
    {'BG_Apr03/data.xlsx'}    {'BG_Apr03/FP_Apr2003.sfit'}

% LD = load(Uz{2}, '-mat')
% sf = LD.savedSession
% stuff1 = sf.CurveFittingToolboxVersion
% stuff2 = sf.AllFitdevsAndConfigs
% stuff2{1}
% stuff2{2}



T1 = readtable(Uz{1})

T1 = 1260×4 table
      X         Y       grid_code       Z    
    ______    ______    _________    ________

    77.544    13.132     295.35      -0.40916
    77.552    13.132     295.33      -0.43155
     77.56    13.132     295.09      -0.67391
    77.568    13.132     294.68       -1.0863
    77.576    13.132     294.58       -1.1887
    77.632    13.132     294.53       -1.2553
    77.639    13.132     294.34       -1.4476
    77.647    13.132     294.19          -1.6
    77.655    13.132     294.39       -1.4024
    77.528    13.124     295.65      -0.11861
    77.536    13.124     295.52      -0.25099
    77.544    13.124     295.45      -0.32334
    77.552    13.124      295.4      -0.37573
     77.56    13.124     295.16      -0.61809
    77.568    13.124      294.7       -1.0805
    77.576    13.124     294.33       -1.4529

[Xd,Xu] = bounds(T1.X);
[Yd,Yu] = bounds(T1.Y);
[Zd,Zu] = bounds(T1.Z);

xv = linspace(Xu, Xd, 100);
yv = linspace(Yd, Yu, 100);
[Xm,Ym] = ndgrid(xv, yv);
ZF = scatteredInterpolant(T1.X, T1.Y, T1.Z);

Zm = ZF(Xm, Ym);

figure
surfc(Xm, Ym, Zm)
colormap(turbo)
xlabel('X')
ylabel('Y')
zlabel('Z')
title('File Data (Interpolated Surface)')

EDIT — (19 Mar 2024 at 02:05)

Added —

% % % b(1) = a0,  b(2) = phi,  b(3) = x0,  b(4) = y0,  b(5) = ax,  b(6) = ay 
% % % xy = [X Y]
f = @(b, xy) b(1)*exp(-(((xy(:,1)-b(3)).*cos(b(2))+(xy(:,2)-b(4)).*sin(b(2))).^2./0.5.*b(5).^2)-(((xy(:,2)-b(4)).*cos(b(2))-(xy(:,1)-b(3)).*sin(b(2))).^2./0.5*b(6).^2));

xy = [T1.X, T1.Y];

opts = optimoptions('lsqcurvefit', 'MaxFunctionEvaluations', 1E+4, 'MaxIterations',1E+4);
[B,rn] = lsqcurvefit(f, [4; 2; 77; 13; 1; 1], xy, T1.Z, [], [], opts);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fprintf('a0  = %8.4f\nphi = %8.4f\nx0  = %8.4f\ny0  = %8.4f\nax  = %8.4f\nay  = %8.4f\n', B)

a0  =   3.5603
phi =  -7.1186
x0  =  77.5578
y0  =  12.9724
ax  =   8.5711
ay  =   9.2511

fprintf('\nResidual Norm = %8.4f\n',rn)

Residual Norm = 423.1419

XY = cat(3,Xm, Ym);
Fmtx = @(b, xy) b(1)*exp(-(((xy(:,:,1)-b(3)).*cos(b(2))+(xy(:,:,2)-b(4)).*sin(b(2))).^2./0.5.*b(5).^2)-(((xy(:,:,2)-b(4)).*cos(b(2))-(xy(:,:,1)-b(3)).*sin(b(2))).^2./0.5*b(6).^2));

Zest = Fmtx(B,XY);

figure
surfc(Xm, Ym, Zm, 'FaceAlpha',0.7)
hold on
surf(Xm, Ym, Zest, 'FaceColor',[1 0 1], 'FaceAlpha',0.7)
hold off
colormap(turbo)
xlabel('X')
ylabel('Y')
zlabel('Z')
title('File Data (Interpolated Surface) With Estimated Surface')

figure
surf(Xm, Ym, Zest)
hold on
surf(Xm, Ym, ones(size(Zest)), 'FaceColor','r', 'FaceAlpha',0.5, 'EdgeColor','none')
hold off
colormap(turbo)
xlabel('X')
ylabel('Y')
zlabel('Z')
title('File Data (Estimated Surface) With Intersecting Plane At Z=1')

figure
[c,h] = contour(Xm, Ym, Zest, [1 1]);
axis('equal')
grid
Ax = gca;
xlabel('X')
ylabel('Y')
title('Contour With Fitted Function')
% Ax.XAxisLocation = 'origin';
% Ax.YAxisLocation = 'origin';

elpsfcn = @(b,xy) xy(1,:).^2/b(1)^2 + xy(2,:).^2/b(2)^2 - b(3);
opts = optimoptions('fminunc', 'MaxFunctionEvaluations', 5E+3, 'MaxIterations',1E+4);
[B, fv] = fminunc(@(b) norm(elpsfcn(b,c(:,2:end))), rand(3,1), opts)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

B = 3×1

104.5118 37.1843 0.6724

fv = 0.0210

fprintf('Semimajor Axis = %8.4f\nSemiminor Axis = %8.4f\nConstant       = %8.4f\n', B)

Semimajor Axis = 104.5118
Semiminor Axis =  37.1843
Constant       =   0.6724

text(min(xlim)+diff(xlim)*0.4, min(ylim)+diff(ylim)-0.05, sprintf('$\\frac{x^2}{%.2f^2} + \\frac{y^2}{%.2f^2} = %.4f$',B), 'Interpreter','latex', 'FontSize',16)

figure
[c,h] = contour(Xm, Ym, Zest, [1 1]);
axis('equal')
grid
Ax = gca;
xlabel('X')
ylabel('Y')
title('Contour With Fitted Function')

elpsfcn = @(b,xy) (xy(1,:)-b(4)).^2/b(1)^2 + (xy(2,:)-b(5)).^2/b(2)^2 - b(3);
opts = optimoptions('fminunc', 'MaxFunctionEvaluations', 5E+3, 'MaxIterations',1E+4);
[B, fv] = fminunc(@(b) norm(elpsfcn(b,c(:,2:end))), [rand(3,1); 77; 13], opts)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

B = 5×1

1.7444 1.7627 0.0026 77.5579 12.9724

fv = 0.0021

fprintf('Semimajor Axis = %8.4f\nSemiminor Axis = %8.4f\nConstant       = %8.4f\n', B(1:3))

Semimajor Axis =   1.7444
Semiminor Axis =   1.7627
Constant       =   0.0026

text(min(xlim)+diff(xlim)*0.4, min(ylim)+diff(ylim)-0.05, sprintf('$\\frac{x^2}{%.2f^2} + \\frac{y^2}{%.2f^2} = %.4f$',B(1:3)), 'Interpreter','latex', 'FontSize',16)

figure
[c,h] = contour(Xm, Ym, Zm, [1 1]);
axis('equal')
grid
Ax = gca;
xlabel('X')
ylabel('Y')
title('Contour With Interpolated Data')
% Ax.XAxisLocation = 'origin';
% Ax.YAxisLocation = 'origin';

elpsfcn = @(b,xy) xy(1,:).^2/b(1)^2 + xy(2,:).^2/b(2)^2 - b(3);
opts = optimoptions('fminunc', 'MaxFunctionEvaluations', 5E+3, 'MaxIterations',1E+4);
[B, fv] = fminunc(@(b) norm(elpsfcn(b,c(:,2:end))), rand(3,1), opts)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

B = 3×1

41.5952 15.6322 0.3215

fv = 459.3697

fprintf('Semimajor Axis = %8.4f\nSemiminor Axis = %8.4f\nConstant       = %8.4f\n', B)

Semimajor Axis =  41.5952
Semiminor Axis =  15.6322
Constant       =   0.3215

text(min(xlim)+diff(xlim)*0.4, min(ylim)+diff(ylim)-0.05, sprintf('$\\frac{x^2}{%.2f^2} + \\frac{y^2}{%.2f^2} = %.4f$',B), 'Interpreter','latex', 'FontSize',16)

NOTE — The surface is a reasonable fit to the data, however there are definitely compromises.

.

[Deepshikha Deo]

Thanks a lot for putting so much efforts. I appreciate the time you have given to my problem.

Yes, the curve fitter app is needed to open the sfit file (the second fit named 'fp' is the one). The grid_code in excel file was used to calculate the values for z column in excel sheet, thus it is not needed further but could not be removed too. I have attached the code which contains the equation (function) with the brief description of its unknown parameters. This function was called in curve fitter app as a custom equation since the real equation was too long to type.

function[zVal]=gaussianFit(x,y,a0,ax,ay,phi,x0,y0)

% no constants in the equation

% the meaning of unknowns is:

% x0 and y0 are the center location of the peak value. Values should be within the range of x and y data

% a0 is the peak value or height of Gaussian curve. Its values should be greater than zero and not more than the maximum z value

% phi is the angle of rotation of ellipse formed.

% ax is the major axis of ellipse or the standard deviation in x-direction. It should be greater than zero

% ay is the minor axis of ellipse or the standard deviation in y-direction. It should be greater than zero

% x and y are the x & y coordinates values

% zVal is the UHI values (z values)

zVal=a0*exp(-(((x-x0).*cos(phi)+(y-y0).*sin(phi)).^2/0.5*ax^2)-(((y-y0).*cos(phi)-(x-x0).*sin(phi)).^2/0.5*ay^2));

end

There is also an onramp course with the name "curve fitting onramp" (https://matlabacademy.mathworks.com/details/curve-fitting-onramp/orcf?s_tid=srchtitle_site_search_1_curve%20fitter%20course) for its understanding.

In the curve fitter app, after adding the data, I selected "custom equation" as fit type and typed the function. After manipulating the lower and upper bounds for six parameters, I hit the 'Fit' which generated the surface fit. The following code is generated by export generate code.

Note- The variables' names may wary.

function [fitresult, gof] = createFit2(POINT_X, POINT_Y, UHI)

%CREATEFIT2(POINT_X,POINT_Y,UHI)

% Create a fit.

%

% Data for 'fp' fit:

% X Input: POINT_X from BG_B_Apr2003

% Y Input: POINT_Y from BG_B_Apr2003

% Z Output: UHI from BG_B_Apr2003

% Output:

% fitresult : a fit object representing the fit.

% gof : structure with goodness-of fit info.

%

% See also FIT, CFIT, SFIT.

% Auto-generated by MATLAB on 19-Mar-2024 05:28:55

%% Fit: 'fp'.

[xData, yData, zData] = prepareSurfaceData( POINT_X, POINT_Y, UHI );

% Set up fittype and options.

ft = fittype( 'gaussianFit(x,y,a0,ax,ay,phi,x0,y0)', 'independent', {'x', 'y'}, 'dependent', 'z' );

opts = fitoptions( 'Method', 'NonlinearLeastSquares' );

opts.Display = 'Off';

opts.Lower = [0 0 0 0 77 12];

opts.StartPoint = [0.754686681982361 0.276025076998578 0.679702676853675 0.655098003973841 0.162611735194631 0.118997681558377];

opts.Upper = [Inf Inf Inf 360 78 14];

% Fit model to data.

[fitresult, gof] = fit( [xData, yData], zData, ft, opts );

% Plot fit with data.

figure( 'Name', 'fp' );

h = plot( fitresult, [xData, yData], zData );

legend( h, 'fp', 'UHI vs. POINT_X, POINT_Y', 'Location', 'NorthEast', 'Interpreter', 'none' );

% Label axes

xlabel( 'POINT_X', 'Interpreter', 'none' );

ylabel( 'POINT_Y', 'Interpreter', 'none' );

zlabel( 'UHI', 'Interpreter', 'none' );

grid on

By using export, the output can be exported to the workspace. The following code is used to plot the exported output surface.

plot(fittedmodel) % fittedmodel is the name of the output surface. This will create plot without the point data values

plot(fittedmodel,[x,y],z) % This will plot the surface with point data values. (the name of x, y, and z may vary)

The surface was generated from zero thus giving me the values for the parameters from zero. This is fine for all values except ax and ay. The concept (suggested in research papers), I followed have shown the values for ax and ay from the value 1 of z-axis. I was not able to find the values when z=1 using curve fitter app as it not should exceed the values of ax and ay at 0 of z-axis. After that I came up with the idea I put in the question which I again failed to solve by myself.

I attached this image as an example for better understanding. It shows the 2-D image of the surface (taked from one of the research papers). In this image, the line ax(1K) represents the standard deviation in x-direction whereas the curve is starting from zero.

I know, this can also be done using optimization too but my lack of knowledge became a hinderance to do so.

Once again thanks a lot for the time you have given to my probelm.

Best Regards

[Star Strider]

I am lost.

The values for ‘ax’ and ‘ay’ are:

ax = 8.5711

ay = 9.2511

I guess I didn’t solve your problem. I will delete my answer if you do not accept it.

[Deepshikha Deo]

I got confused with the following values

Semimajor Axis = 107.8973

Semiminor Axis = 41.3402

As it should not be greater than values of ax and ay since these two values are represents axes values at the base of the gaussian surface (when the surface starts from zero).

Technically, you have answered the question but I could not create the surface with my already existing data. Please, tell me how to create that grid to create surface if I already have x and y values as (n x 1) size variable.

Regards

[Star Strider]

The semimajor and semimnor axis values are with respect to the data provided. Since the data values are relatively large in magnitude (77 for ‘X’ and 13 for ‘Y’), the corresponding values are large.

Running it a few more times, and adding a second ellipse parameter estimation with offsets for ‘X’ and ‘Y’, and the objective function now defined as:

elpsfcn = @(b,xy) (xy(1,:)-b(4)).^2/b(1)^2 + (xy(2,:)-b(5)).^2/b(2)^2 - b(3);

with the ‘b(4)’ and ‘b(5)’ offsets themselves (77.5579, 12.9724) appropriately estimated, gives a much better result:

Semimajor Axis = 1.7444

Semiminor Axis = 1.7627

Constant = 0.0026

This ratio is now 0.99, and more representative of the observed ellipse shape. (See the second ‘Contour With Fitted Function’ plot.)

‘Please, tell me how to create that grid to create surface if I already have x and y values as (n x 1) size variable.’

I am not certain what you want, however I probably already provided that in my code, specifically:

xv = linspace(Xu, Xd, 100);

yv = linspace(Yd, Yu, 100);

[Xm,Ym] = ndgrid(xv, yv);

ZF = scatteredInterpolant(T1.X, T1.Y, T1.Z);

Zm = ZF(Xm, Ym);

These are excerpted from the code that created the first surf plot. They interpolate the existing data and create it as the surface.

That is likely the best that I can do with these data.

.

[Deepshikha Deo]

Thank you for your efforts.

Regards

[Star Strider]

As always, my pleasure!