How to do Fminsearch- Curve fitting

I am trying to perfrom Fminsearch to fit the image intensity decaying data points to the wave propogation model , to determine constants

I am not sure how to go about the fminsearch. Kindly requesting help.

Below is my code:

% Objective: 
 %{
   1. Import a CT scan image( CT scan of an Cu wire)
   2. Extract the intensity from the CT scan image(CT scan of an Cu wire) :
   which will be the observed data
   3. Define the exponential decay model describing the Wave propagation through a medium
   4. Perform fminsearch to determine the constants ( imaginary constants)
   in the decay model.
%}

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;


%% Importing the CT scan/test image

% fullname = get_full_filename('20220422_5anglestep000.tif')
test_image = imread('test_image.tif');
figure(1)
imshow(test_image)

% Getting the dimensions of the image.
[rows, columns] = size(test_image);



%% % Computung and displaying the histogram.
figure(5)
[pixelCount, grayLevels] = imhist(test_image);
 bar(pixelCount);
grid on;
title('Histogram of test image');
xlim([0 grayLevels(end)]); % Scale x axis manually.

%% Extracting the intensities from the medium

% Image represents the wave propgation through a medium in terms of
% intensity vs thickness(displacement) of the medium

% Looking to extract the intensity as a function of displacment from image

% k: displacement ( Z axis locations)
% wave_obs = intensity (Y intenisity)


figure(3)
middleRow = round(size(test_image, 1) / 2) % Find middle row.
wave_obs = test_image(middleRow, :) % Get all columns in the middle row.
% Now plot the intensity.
plot(wave_obs, 'b-', 'LineWidth', 2);
grid on;
xlabel('X Column');
ylabel('Gray Level');


%% Fitting wave_obs data points to the below defined exponential decay model  

%% Defining the equation wave = exp(i(1-del)k*z * exp(-beta)*z

%{
i -  signifying imaginary part
z -  array representing the length of the medium : # of wave_obs data points
k -  refractive index
del - constant (to be determined)
beta -constant (to be determined)
wave_obs- data points to fit

%}

k=1.5;

%% fminsearch : HELP NEEDED!

wave_pred =  fminsearch(wave_prediction);

predicted_wave = wave_pred(z,k)



function wave_prediction = wave_pred(z,k)
    del =5;
    beta =7;
    k=1.5;
    z= 0:0.2:20;
    wave1 = exp(-1i*(1-del)*k*z);
    wave2 = exp(-beta)*z;
    wave_prediction = @(del,beta)wave1.*wave2;
end

[Anand Ra]

@Walter Roberson kindly requesting your help for the above

[Walter Roberson]

wave_pred =  fminsearch(wave_prediction);

At that point in the code, you do not have any function or variable named wave_prediction so it is not obvious what you are trying to do there.

Remember also that fminsearch() needs an initial search point.

predicted_wave = wave_pred(z,k)

I wonder if you should be calling that first and assigning the result to wave_prediction, since the output of the function is a function handle ?

    wave_prediction = @(del,beta)wave1.*wave2;

That function handle accepts two parameters, del and beta, and ignores them both and always returns the product of wave1 and wave2. But if you are planning to make that the function handle used by fminsearch() then you have the difficulty that fminsearch() will only pass in a single parameter (that is a vector)

wave1 and wave2 are both defined in terms of the vector z, so the result of executing the function handle would be a vector. But fminsearch() needs a scalar to be returned.

You are trying to do a curve fitting, which requires that you have existing known data points and known corresponding value. It is not immediately clear to me which variables in your code hold the existing values ?

[Anand Ra]

@Walter Roberson, Thanks, trying to understand your response. Btw the variable wave_obs is holding the data points to be fitted.

Answer 1

[Matt J]

Currently, you are minimizing your prediction function. You need to be minimizing, not the prediction, but a function that compares your prediction to measured data. You need to provide the measured data as well as an initial gues [del0,beta0] if the unknown parameters

fminsearch(@(p,meas) wave_prediction_Error(p, meas)  ,  [del0,beta0]);

function Error = wave_prediction_Error(params, measurements)

    del =params(1);
    beta =(2);
    
    k=1.5;
    z= 0:0.2:20;

    waveprediction = exp(-1i*(1-del)*k*z) .* exp(-beta)*z;
 
    
    Error = norm( waveprediction - measurements);
end

[Anand Ra]

@Matt J Sorry, I am trying to follow your response and struggling. The data I am trying to fit is stored in wave_obs ( attached). I can consider the initial values of del and beta as 0.01.

below is the changes I made based on your feedback:

del0=0.1;
beta0=0.1;

fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs)  ,  [del0,beta0]);


figure(9)
 z= 0:0.2:20;
plot(z, wave_obs, 'ro', z, waveprediction, 'b-');

function Error = wave_prediction_Error(params, wave_obs)
    del =params(1);
    beta =(2);
    
    k=1.5;
    z= 0:0.2:20;
    waveprediction = exp(-1i*(1-del)*k*z) .* exp(-beta)*z;
 
    
    Error = norm( waveprediction - wave_obs);
end

[Matt J]

And? What was the result?

[Anand Ra]

@Matt J oops! Sorry forget to attach the error. Trying to troubleshoot.

Error using *

Incorrect dimensions for matrix

multiplication. Check that the number

of columns in the first matrix

matches the number of rows in the

second matrix. To perform elementwise

multiplication, use '.*'.

Error in

absorption_model>@(p,meas)wave_prediction_Error(p,wave_obs)

(line 120)

fminsearch(@(p,meas)

wave_prediction_Error(p, wave_obs) ,

[del0,beta0]);

Error in fminsearch (line 201)

fv(:,1) = funfcn(x,varargin{:});

Error in absorption_model (line 120)

fminsearch(@(p,meas)

wave_prediction_Error(p, wave_obs) ,

[del0,beta0]);

Related documentation

[Matt J]

I made a few modifications to your code below. However, the fundamental problem is that your wave_prediction only has 101 points, but your wave_obs has much more. Therefore, we cannot do the vector subtraction in this line

Error = norm( waveprediction - wave_obs);

load('wave_obs.mat', 'wave_obs')

wave_obs=double(wave_obs);

del0=0.1;
beta0=0.1;

fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs)  ,  [del0,beta0]);

Arrays have incompatible sizes for this operation.

Error in solution>wave_prediction_Error (line 24)
Error = norm( waveprediction - wave_obs);

Error in solution (line 8)
fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs) , [del0,beta0]);

Error in fminsearch (line 201)
fv(:,1) = funfcn(x,varargin{:});

figure(9)
 z= 0:0.2:20;
plot(z, wave_obs, 'ro', z, waveprediction, 'b-');

function Error = wave_prediction_Error(params, wave_obs)
    del =params(1);
    beta =(2);
    
    k=1.5;
    z= 0:0.2:20;
    wave_prediction = exp(-1i*(1-del)*k*z) .* exp(-beta).*z;
 
    
    Error = norm( wave_prediction - wave_obs);
end

[Matt J]

It is also suspicious that your wave_obs are purely real-valued, but your wave_prediction is complex-valued. Therefore, it is hard to see how your model could ever agree with wave_obs.

[Anand Ra]

@Matt J

Actually modified the input data. Attached the new data set.

As far as fitting with your code is concerned,

1. Looking to plot observed vs predicted data for comparison.
2. And how to obtain values for the constants del and beta.

Below is the code

load('wave_obs.txt', 'wave_obs')

wave_obs=double(wave_obs);

del0=0;

beta0=0;

fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs) , [del0,beta0]);

% z= 0:1:1449;

% plot(z,wave_obs)

figure(1)

function Error = wave_prediction_Error(params, wave_obs)

del =params(1);

beta =(2);

z= 0:1:1450;

k=1.5;

wave_prediction = exp(-1i*(1-del)*k*z) .* exp(-beta).*z;

Error = norm( wave_prediction - wave_obs);

% plot(z,wave_obs)

plot(z,wave_obs,'r-',z,wave_prediction,'b-');

end

[Matt J]

The optimal del and beta will be the output of fminsearch:

params=fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs)  ,  [del0,beta0]);

delOpt=params(1);
betaOpt=params(2)

wave_prediction=@(z)exp(-1i*(1-delOpt)*k*z) .* exp(-betaOpt).*z;

However, there is still the problem that your prediction function produces complex values, so what would it mean to "plot" it? You could certainly plot the absolute values:

z= 0:1:1450;

plot(z,abs(wave_obs),'rx',z,abs( wave_prediction(z)) 'b-')

[Anand Ra]

@Matt J Thank you. You have helped me to correct a couple of things by forcing to think about fitting. -The del and beta values are not zero. And yes, you are right, I am looking to plot the absolute values.

However, not sure why the code keeps running and the plot displayed is disappeared.

The final values of del and beta are not displayed.

load('wave_obs.txt', 'wave_obs')

wave_obs=double(wave_obs);

del0=1.6419807E-06; % updated

beta0=7.30026262E-10;%updated

% fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs) , [del0,beta0]);

params=fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs) , [del0,beta0]);

delOpt=params(1);

betaOpt=params(2);

% wave_prediction=@(z)exp(-1i*(1-delOpt)*k*z) .* exp(-betaOpt).*z;

% z= 0:1:1449;

% plot(z,wave_obs)

figure(1)

function Error = wave_prediction_Error(params, wave_obs)

delOpt=params(1);

betaOpt=params(2);

z= 1:1:1451;

k=(2*pi)/0.00000008;

wave_prediction = exp(-1i*(1-del)*k*z) .* exp(-beta).*z;

Error = norm( wave_prediction - wave_obs);

figure(2);

plot(z,abs(wave_obs),'rx',z,abs(wave_prediction(z)),'b-')

% plot(z,wave_obs,'r-',z,wave_prediction,'b-');

grid on

grid minor

legend('Observed','Predicted');

xlabel('Displacement (Thickness'); ylabel('Intensity');

title('exponential decay model describing the Wave propagation through a medium');

end

[Matt J]

For some reason, you moved the plotting commands inside your objective function. That's not where it belongs:

load('wave_obs.txt', 'wave_obs')
wave_obs=double(wave_obs);



 del0=1.6419807E-06; % updated
 beta0=7.30026262E-10;%updated
% fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs)  ,  [del0,beta0]);

params=fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs)  ,  [del0,beta0]);

delOpt=params(1);
betaOpt=params(2);



final_wave_prediction=@(z)exp(-1i*(1-delOpt)*k*z) .* exp(-betaOpt).*z;

figure(2);
plot(z,abs(wave_obs),'rx',z,abs(final_wave_prediction(z)),'b-')
grid on
grid minor
legend('Observed','Predicted');
xlabel('Displacement (Thickness'); ylabel('Intensity');
title('exponential decay model describing the Wave propagation through a medium');



function Error = wave_prediction_Error(params, wave_obs)
      
    delOpt=params(1);
    betaOpt=params(2);

    z= 1:1:1451;
    k=(2*pi)/0.00000008;
   
    wave_prediction = exp(-1i*(1-del)*k*z) .* exp(-beta).*z;
    Error = norm( wave_prediction - wave_obs);

 

end
    

[Anand Ra]

@Matt J For some reason running the above code, throws below errors.. I am unable to determine the fixes.

'del' requires ROS Toolbox.

Error in absorption_fitting>wave_prediction_Error

(line 38)

wave_prediction = exp(-1i*(1-del)*k*z) .*

exp(-beta).*z;

Error in

absorption_fitting>@(p,meas)wave_prediction_Error(p,wave_obs)

(line 11)

params=fminsearch(@(p,meas) wave_prediction_Error(p,

wave_obs) , [del0,beta0]);

Error in fminsearch (line 201)

fv(:,1) = funfcn(x,varargin{:});

Error in absorption_fitting (line 11)

params=fminsearch(@(p,meas) wave_prediction_Error(p,

wave_obs) , [del0,beta0]);

[Matt J]

Change this back to what we had before (plus a few tweaks that I've inserted),

function Error = wave_prediction_Error(params, wave_obs)
      
    del=params(1);
    beta=params(2);
    z= 1:numel(wave_obs);
    k=(2*pi)/0.00000008;
   
    wave_prediction = exp(-1i*(1-del)*k*z) .* exp(-beta).*z;
    Error = norm( wave_prediction(:) - wave_obs(:));
 
end

[Anand Ra]

@Matt J Ok and did that. The predicted values on the plot is zero. And I am unable to verify the final estimates of del and beta.

I hope I incorporated the changes you mentioned correctly. Below is the incorporated code.

%% Loading the intensity data

load('wave_obs.txt', 'wave_obs')

wave_obs=double(wave_obs);

%% Defining the model of wave decay = exp(i(1-del)k*z * exp(-beta)*z

%{

i - signifying imaginary part

z - array representing the length of the medium

k - refractive index

del - constant ( to be determined)

beta -constant ( to be determined)

%}

%% Initial values of del and beta

del0=1.6419807E-06; % updated

beta0=7.30026262E-10;%updated

% fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs) , [del0,beta0]);

params=fminsearch(@(p,meas) wave_prediction_Error(p, wave_obs) , [del0,beta0]);

delOpt=params(1);

betaOpt=params(2);

% Objective function

final_wave_prediction=@(z)exp(-1i*(1-delOpt)*k*z) .* (exp(-betaOpt).*k*z);

% Plotting

figure(2);

plot(z,abs(wave_obs),'rx',z,abs(final_wave_prediction(z)),'b-')

grid on

grid minor

legend('Observed','Predicted');

xlabel('Displacement (Thickness)');

ylabel('Intensity');

title('exponential decay model describing the Wave propagation through a medium')

function Error = wave_prediction_Error(params, wave_obs)

del=params(1);

beta=params(2);

z= 1:numel(wave_obs);

k=(2*pi)/0.00000008;

wave_prediction = exp(-1i*(1-del)*k*z) .* exp(-beta).*z;

Error = norm( wave_prediction(:) - wave_obs(:));

end

[Matt J]

Your model is still wrong. It is still complex-valued whereas wave_obs is real.

[Anand Ra]

Attached is the snapshot of the model.

[Anand Ra]

@Matt J Attached (above) the model. At the top you can see the model equation.

[Matt J]

It's complex. Your measurements are real.

[Matt J]

Also, the absolute value of the model function in your photo cannot exceed 1 for beta*k>0. However, wave_obs has a minimum value of,

>> min(wave_obs)

ans =

1500000

[Matt J]

Thirdly, it seems unlikely that in the graph in your photo z is a sequence of integers 1:N like it is in your code

Surely z is supposed to have units as well?

[Anand Ra]

@Matt J

I found this mathworks reference for fitting data to a complex function. But I am struggling to decode it.

https://www.mathworks.com/matlabcentral/answers/99986-how-can-i-fit-a-complex-valued-function-to-my-data

[Anand Ra]

Agreed about z.

Quick background summary.

The data set I am trying to fit is the intensity data produced by the CT scan of a copper wire.

So to answer your question about "z" , it would be the thickness of the copper wire which is about 1.8-2mm.

[Anand Ra]

@Matt J .

Continuing previous comment. : As you can see the data set I have is decaying when plotted. So, I am trying to fit this data set to this decaying model.

[Matt J]

No, the figure shows that the decaying curves follow equations and which are not the model equation curves you have implemented.

Also, the curves decay starting from 1 at z=0. Your wave_obs data are all well above 1. Presumably, there was some sort of data normalization that you were supposed to do.

Also, the data you have attached are not monotonically decaying. They have oscillations:

[Anand Ra]

@Matt J Thanks a lot for the support. Very much appreciated. I will look at the math again to see if I can simplify to work on the fitting.