Fit a custom made function to a certain trend within data from a matrix

Dear all,

I have been struggling for a while with the two issues that I show you below.

I have a matrix which represents the image of a pixelated neutron detector. The units of such a matrix is Count/minute. If I plot the matrix with the command surface(X, Y, Z), being Z the matrix, I get the following figure:

In this figure we can see two different trends: the x=y line (plus its width), which is what we call the "Specular Ridge", and a tiny contribution given by the white dashed line, which is what we call the "Zemann Splitting". The Zemann splitting is given by the equation y = sqrt(x^2 - (1.47e-7*H*L^2)), where H is the magnetic field and L is the wavelength of the nuetrons (5.183 AA in my case).

I would like to find the optimum H which fits best to the data inside the red circle by using the above mentioned equation (H should be something like 6-7). And later I would like to make a ROI (Region Of Interest) to count how many counts lay within such area.

I have attached the matrix (file 00607_UD_subplot.txt) and X and Y coordinates (files 00607_UD_xcoordinates.txt and 00607_UD_ycoordinates.txt).

The code to visualize the figure is the following:

A = load('00607_UD_subplot.txt');

X = load('00607_UD_xcoordinates.txt');

Y = load('00607_UD_ycoordinates.txt');

surface(X, Y, A)

axis square

colormap(map);

caxis([5e-4 100])

set(gca,'colorscale','log')

hold on

plot(X, sqrt((X.^2) - (1.47e-7*6.5*5.183^2)), 'Color','white','LineStyle','--', 'LineWidth', 2)

xlabel('\alpha_{i} [rad]', 'fontsize', 18);

ylabel('\alpha_{f} [rad]', 'fontsize', 18);

cb1 = colorbar;

cb1.Label.String = 'Neutron Intensity (CTS/MIN)';

cb1.FontSize = 16;

cb1.Label.FontSize = 20;

cb1.LineWidth = 1.5;

Any help / suggestion would be helpful and very welcomedl!!!

Cheers,

Jose.

Answer 1

[Star Strider]

TThe problem here is how to isolate the values you want to regress. I worked on this for a while, and then hit on the idea of using a contour to define the central area bounds (and then fudged those values a bit in the end) to isolate only the values you want tto regress. I plotted them separately at the end, not on the original plot.

Try this —

A = load('00607_UD_subplot.txt');
X = load('00607_UD_xcoordinates.txt');
Y = load('00607_UD_ycoordinates.txt');

[Xm,Ym] = meshgrid(X, Y);
[x1,x2] = bounds(X)

x1 = 0

x2 = 0.0147

[y1,y2] = bounds(Y)

y1 = 2.5368e-05

y2 = 0.0149

[a1,a2] = bounds(A,'all')

a1 = 0

a2 = 0.3331

BL = [x1 1; x2 1] \ [y1; y2]

BL = 2×1

1.0173 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

yfcn = @(b,x) sqrt(x.^2 - (1.47e-7*b(1)*b(2).^2));

Lx = (Xm >= 0.005) & (Xm <= 0.0075);

figure
surface(X, Y, A)
axis square
colormap(turbo);
caxis([5e-4 100])
set(gca,'colorscale','log')
hold on
plot(X, sqrt((X.^2) - (1.47e-7*6.5*5.183^2)), 'Color','white','LineStyle','--', 'LineWidth', 2)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

xlabel('\alpha_{i} [rad]', 'fontsize', 18);
ylabel('\alpha_{f} [rad]', 'fontsize', 18);
cb1 = colorbar;
cb1.Label.String = 'Neutron Intensity (CTS/MIN)';
cb1.FontSize = 16;
cb1.Label.FontSize = 20;
cb1.LineWidth = 1.5;

Lvl = 0.020;

figure
surf(Xm, Ym, A, 'EdgeColor','interp')
hold on
[cm,ch] = contour3(Xm, Ym, A, [1 1]*Lvl, '-r', 'LineWidth',2);
hold  off
xlabel('\alpha_{i} [rad]', 'fontsize', 18);
ylabel('\alpha_{f} [rad]', 'fontsize', 18);
colormap(turbo)

find(cm(1,:) == Lvl)

ans = 1×2

1 8

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

S = [numel(cm(1,9:end)) numel(unique(cm(1,9:end)))]

S = 1×2

159 120

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(cm(1,2:7), cm(2,2:7), '.-')
hold on
plot(cm(1,9:end), cm(2,9:end), '.-')
hold off
grid
axis('equal')

[gx,gy] = gradient([cm(1,9:end); cm(2,9:end)].');

figure
quiver(cm(1,9:end).', cm(2,9:end).', -gx(:,1), -gx(:,2), 0.25)
grid
axis('equal')

[xm,xmix] = min(cm(1,9:end))

xm = 7.2756e-04

xmix = 80

xc = cm(1,xmix:end);
yc = cm(2,xmix:end);

Lc = (xc >= 0.005) & (xc <= 0.0075);
xcv = xc(Lc);
ycv = yc(Lc);

ps = polyshape([xcv flip(xcv)], [zeros(size(ycv)) flip(ycv)])

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

ps = 
  polyshape with properties:

      Vertices: [16x2 double]
    NumRegions: 1
      NumHoles: 0

AL = A > 0;
Xmz = Xm(AL);
Ymz = Ym(AL);
Az = A(AL);
inp = inpolygon(Xmz, Ymz, [xcv flip(xcv)], [zeros(size(ycv)) flip(ycv)]-0.0005);
nip = nnz(inp)

nip = 11

Xmip = Xmz(inp);
Ymip = Ymz(inp);
Aip = A(inp);

yfcn = @(b,x) sqrt(x.^2 - (1.47e-7*b(1)*b(2).^2));
[B,fv] = fminsearch(@(b) norm(Ymip - yfcn(b,Xmip)), rand(2,1)*10);
resnorm = fv

resnorm = 0.0038

fprintf('\nH = %.5f\nL = %.5f\n',B)

H = 7.10612
L = 4.92885

xp = linspace(min(Xmip)-1.5E-4, max(Xmip));

yp = yfcn(B, xp);

figure
scatter(Xmip, Ymip, 's', 'filled', 'DisplayName','Data')
hold on
plot(xp, yp, '--r', 'DisplayName','Regression Fit')
hold off
grid
axis('equal')
legend('Location','best')

This creates a polygon from the lower contour line and the x-axis, then isolates the points using the inpolygon function.

The regression is straightforward, however other optimization toolbox functions may give more accurate results for the parameters.

NOTE — The parameters are multiplied by each other, so their estimated values are are not entirely independent. (This is not a good model, acttually, and ideally should have a second equation that also allows the parameters to be defined to be more independent.)

.

[Jose Manuel Gomez Guzman]

Works perfectly fine. Thanks!!!

[Star Strider]

As always, my pleasure!

Yours is probably the most challenging problem I have seen here, at least recently.