Fitting Curve with scattered dots

Hi,

I have to plot the beam propagation of a Laser. I can plot the scattered dots, which we have measured in our experiments. Now I wanna Fit those dots. But it doesn't work. It says "Error using fit>iFit Insufficient data. You need at least 3 data points to fit this model." But I already have 22 data points. Can anybody help me?

close all;
clear variables;
%Files and Pics
data_list = dir('*.bmp');
%Zugreifen auf Ordner und BIlder und Namen zuweisen
name = {data_list.name};
%Laden der schwarzen Filterdatei
%black = uigetfile("E:\Master_Technologiemanagement\Spezialisierung\Technik\Laser\Praktikum\SB_V2\Pictures\dunkel\460_dunkel.bmp");
%Konvertieren in Matrix
imBlack = imread("E:\Master_Technologiemanagement\Spezialisierung\Technik\Laser\Praktikum\SB_V2\Pictures\dunkel\460_dunkel.bmp");
lambda = 0.633; %Wellenlänge Laser in um
x_array = [];
y_array = [];
z_array = [];
for i = 1:length(name)
    picture = imread(name{i});  %Einlesen des Bildes in Matrix
    current_name = name{i};
    zaxe = str2double(name{i}(1:3)); %Auslesen des z Wertes aus Dateinamen
    
    imBlack = double(imBlack);
    %picture(:,:,i)=picture;    %speichert das aktuelle (i-te) 
    %Bild als 2D Matrix in der 3D Matrix (:,:,i)
    picture = double(picture);
    
    filtered = minus(picture,imBlack);
    %Ermittlung ROI
    ROI = im2double(filtered);
    ROIdiff = [1,1];
    ROIoffset =[0, 0];
    %1.Moment
    [row,column] = size(filtered);
    Summe = 0;
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*x;
        end
    end
    first_Moment_X = Summe;
    Summe = 0;
    Integral = sum(sum(filtered));
    Schwerp_X = first_Moment_X./Integral;
    Schwerp_X = double(Schwerp_X);
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*y;
        end
    end
    first_Moment_Y=Summe;
    Schwerp_Y = first_Moment_Y/Integral;
    %2.Moment
    [row,column] = size(filtered);
    Summe = 0.;
    Integral = sum(sum(filtered));
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(x-Schwerp_X).^2;
        end
    end
    sec_Moment_X = Summe;
    Summe = 0;
    Sigma_X = sec_Moment_X./Integral;
    Sigma_X = double(Sigma_X);
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(y-Schwerp_Y).^2;
        end
    end
    sec_Moment_Y=Summe;
    Summe = 0;
    Sigma_Y = sec_Moment_Y./Integral;
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(y-Schwerp_Y).*(x-Schwerp_X);
        end
    end
    sec_Moment_XY=Summe;
    Summe = 0;
    Sigma_XY = sec_Moment_XY./Integral;
    %Berechnung des Strahldurchmessers nach Formel 18 und 19 in Norm 
    dxPixel = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y+(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4*(Sigma_XY.^2))));
    dyPixel = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y-(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4.*(Sigma_XY.^2))));
    
    %Gesamtbild mit ROI, Schwerpunkt und Durchmesser
    % imshow(picture);
    axis on
    hold on;
    %  rectangle('position',[ROIoffset(2), ROIoffset(1),2*dxPixel,2*dyPixel],'EdgeColor','b');
    %  plot(ROIoffset(2)+SchwerpX, ROIoffset(1)+Schwerp_Y,'r+', 'MarkerSize', 5, 'LineWidth', 1);
    t=-pi:0.01:pi;
    xOff=ROIoffset(2)+Schwerp_X+dxPixel.*0.5.*cos(t);
    yOff=ROIoffset(1)+Schwerp_Y+dyPixel.*0.5.*sin(t);
    % plot(xOff,yOff,'g')
    
    hold on;
    scatter(zaxe,dxPixel,50,'x','MarkerFacecolor','blue');
    
end
%Plotten der einzelnen Bilder
hold on;
[fitx,gofx]= fit(zaxe,dxPixel,"poly2");
fithandle = plot(fitx);
set(fithandle,'color','[0 .7 .7]');
[fitxw,gofxw] = fit(zaxe,dxPixel,'poly2','weights',1./dxPixel);

[Cris LaPierre]

Share your data. You can attach it to your post using the paperclip icon.

[Bersin Tekmen]

@Cris LaPierre done

[Matt J]

@Bersin Tekmen I suggest you attach a .mat file containing zaxe and dxPixel.

[Bersin Tekmen]

@Matt J all i have about zaxe and dxPixel are those 3 lines of code here attached above..

[Matt J]

@Berkin Bilgic which 3 lines and why can't you attach the results that those lines produce?

Answer 1

[Cris LaPierre]

There are a couple issues to fix.

The error you are currently getting is because zaxe and dxPixel are scalars. In order to use them as inputs to the fit function, they need to be column vectors. That is fixed by using indexing to capture the result of each loop and assign it to a specific row of the corresponding vector.

zaxe(i,1) = str2double(name{i}(1:3)); %Auslesen des z Wertes aus Dateinamen

...

dxPixel(i,1) = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y+(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4*(Sigma_XY.^2))));

dyPixel(i,1) = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y-(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4.*(Sigma_XY.^2))));

You will then get a new error that one of your values (in dxPixel) is NaN, which is not allowed. You can fix that with a call to prepareCurveData.

[zaxe,dxPixel] = prepareCurveData(zaxe,dxPixel);

The final code might look like this.

% Load images
unzip('bmps.zip')

%Files and Pics
data_list = dir('*.bmp');
%Zugreifen auf Ordner und BIlder und Namen zuweisen
name = {data_list.name};

%Konvertieren in Matrix
imBlack = imread("460_dunkel.bmp");
lambda = 0.633; %Wellenlänge Laser in um
x_array = [];
y_array = [];
z_array = [];
for i = 1:length(name)
    picture = imread(name{i});  %Einlesen des Bildes in Matrix
    current_name = name{i};
    zaxe(i,1) = str2double(name{i}(1:3)); %Auslesen des z Wertes aus Dateinamen

    imBlack = double(imBlack);
    %picture(:,:,i)=picture;    %speichert das aktuelle (i-te)
    %Bild als 2D Matrix in der 3D Matrix (:,:,i)
    picture = double(picture);

    filtered = minus(picture,imBlack);
    %Ermittlung ROI
    ROI = im2double(filtered);
    ROIdiff = [1,1];
    ROIoffset =[0, 0];
    %1.Moment
    [row,column] = size(filtered);
    Summe = 0;
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*x;
        end
    end
    first_Moment_X = Summe;
    Summe = 0;
    Integral = sum(sum(filtered));
    Schwerp_X = first_Moment_X./Integral;
    Schwerp_X = double(Schwerp_X);
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*y;
        end
    end
    first_Moment_Y=Summe;
    Schwerp_Y = first_Moment_Y/Integral;
    %2.Moment
    [row,column] = size(filtered);
    Summe = 0.;
    Integral = sum(sum(filtered));
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(x-Schwerp_X).^2;
        end
    end
    sec_Moment_X = Summe;
    Summe = 0;
    Sigma_X = sec_Moment_X./Integral;
    Sigma_X = double(Sigma_X);
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(y-Schwerp_Y).^2;
        end
    end
    sec_Moment_Y=Summe;
    Summe = 0;
    Sigma_Y = sec_Moment_Y./Integral;
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(y-Schwerp_Y).*(x-Schwerp_X);
        end
    end
    sec_Moment_XY=Summe;
    Summe = 0;
    Sigma_XY = sec_Moment_XY./Integral;
    %Berechnung des Strahldurchmessers nach Formel 18 und 19 in Norm
    dxPixel(i,1) = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y+(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4*(Sigma_XY.^2))));
    dyPixel(i,1) = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y-(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4.*(Sigma_XY.^2))));

    %Gesamtbild mit ROI, Schwerpunkt und Durchmesser
    t=-pi:0.01:pi;
    xOff=ROIoffset(2)+Schwerp_X+dxPixel(i).*0.5.*cos(t);
    yOff=ROIoffset(1)+Schwerp_Y+dyPixel(i).*0.5.*sin(t);

end

scatter(zaxe,dxPixel,50,'x','MarkerFacecolor','blue');
[zaxe,dxPixel] = prepareCurveData(zaxe,dxPixel);

Warning: Removing NaN and Inf from data

%Plotten der einzelnen Bilder
hold on;
[fitx,gofx]= fit(zaxe,dxPixel,"poly2");
fithandle = plot(fitx);
set(fithandle,'color','[0 .7 .7]');
[fitxw,gofxw] = fit(zaxe,dxPixel,'poly2','weights',1./dxPixel);
hold off

[Bersin Tekmen]

Damn I am so happy right now. Thank you so much. This one worked for me. Thank you.

Answer 2

[Voss]

The problem is that zaxe and dxPixel are overwritten on each iteration of the for loop, so after the loop finishes, they have the value they had on the last iteration. That's a problem for fit because they are both scalars.

If you want to use all zaxe and dxPixel values used in each iteration of the for loop you can make them into vectors with one element per for-loop iteration, and use those vectors in fit after the loop. Here I've done that and the vectors are called zaxe_all and dxPixel_all.

(Also, note that the loop iterates over all .bmp files, including 460_dunkel.bmp and 560_fokus.bmp. I got a NaN somewhere that fit didn't like, corresponding to 460_dunkel.bmp, so I removed that from name. It is unclear whether 560_fokus.bmp should also be removed.)

(Another thing I want to point out is you can use zaxe_all and dxPixel_all to create your scatter plot after the loop, rather than doing a scatter plot for each iteration inside the loop. However, you're specifying a constant MarkerFaceColor but using a Marker that's not affected by MarkerFaceColor, so I wasn't sure what the intent was, so I didn't change that.)

close all;
clear variables;
%Files and Pics
data_list = dir('*.bmp');
%Zugreifen auf Ordner und BIlder und Namen zuweisen
name = {data_list.name}

name = 1×23 cell array
    {'460.bmp'}    {'460_dunkel.bmp'}    {'470.bmp'}    {'480.bmp'}    {'490.bmp'}    {'500.bmp'}    {'540.bmp'}    {'544.bmp'}    {'548.bmp'}    {'552.bmp'}    {'556.bmp'}    {'560_fokus.bmp'}    {'564.bmp'}    {'568.bmp'}    {'572.bmp'}    {'576.bmp'}    {'580.bmp'}    {'582.bmp'}    {'620.bmp'}    {'630.bmp'}    {'640.bmp'}    {'650.bmp'}    {'660.bmp'}

%Laden der schwarzen Filterdatei
%black = uigetfile("E:\Master_Technologiemanagement\Spezialisierung\Technik\Laser\Praktikum\SB_V2\Pictures\dunkel\460_dunkel.bmp");
%Konvertieren in Matrix
% imBlack = imread("E:\Master_Technologiemanagement\Spezialisierung\Technik\Laser\Praktikum\SB_V2\Pictures\dunkel\460_dunkel.bmp");
imBlack = imread("460_dunkel.bmp");

name(strcmp(name,'460_dunkel.bmp')) = []; % remove 460_dunkel.bmp

zaxe_all = zeros(size(name));
dxPixel_all = zeros(size(name));

lambda = 0.633; %Wellenlänge Laser in um
x_array = [];
y_array = [];
z_array = [];
for i = 1:length(name)
    picture = imread(name{i});  %Einlesen des Bildes in Matrix
    current_name = name{i};
    zaxe = str2double(name{i}(1:3)); %Auslesen des z Wertes aus Dateinamen
    
    imBlack = double(imBlack);
    %picture(:,:,i)=picture;    %speichert das aktuelle (i-te) 
    %Bild als 2D Matrix in der 3D Matrix (:,:,i)
    picture = double(picture);
    
    filtered = minus(picture,imBlack);
    %Ermittlung ROI
    ROI = im2double(filtered);
    ROIdiff = [1,1];
    ROIoffset =[0, 0];
    %1.Moment
    [row,column] = size(filtered);
    Summe = 0;
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*x;
        end
    end
    first_Moment_X = Summe;
    Summe = 0;
    Integral = sum(sum(filtered));
    Schwerp_X = first_Moment_X./Integral;
    Schwerp_X = double(Schwerp_X);
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*y;
        end
    end
    first_Moment_Y=Summe;
    Schwerp_Y = first_Moment_Y/Integral;
    %2.Moment
    [row,column] = size(filtered);
    Summe = 0.;
    Integral = sum(sum(filtered));
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(x-Schwerp_X).^2;
        end
    end
    sec_Moment_X = Summe;
    Summe = 0;
    Sigma_X = sec_Moment_X./Integral;
    Sigma_X = double(Sigma_X);
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(y-Schwerp_Y).^2;
        end
    end
    sec_Moment_Y=Summe;
    Summe = 0;
    Sigma_Y = sec_Moment_Y./Integral;
    for y=1:row
        for x=1:column
            Summe = Summe + double(filtered(y,x))*(y-Schwerp_Y).*(x-Schwerp_X);
        end
    end
    sec_Moment_XY=Summe;
    Summe = 0;
    Sigma_XY = sec_Moment_XY./Integral;
    %Berechnung des Strahldurchmessers nach Formel 18 und 19 in Norm 
    dxPixel = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y+(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4*(Sigma_XY.^2))));
    dyPixel = 2*sqrt(2)*sqrt(Sigma_X + Sigma_Y-(sign(Sigma_X-Sigma_Y).*sqrt((Sigma_X-Sigma_Y).^2+4.*(Sigma_XY.^2))));
    
    %Gesamtbild mit ROI, Schwerpunkt und Durchmesser
    % imshow(picture);
    axis on
    hold on;
    %  rectangle('position',[ROIoffset(2), ROIoffset(1),2*dxPixel,2*dyPixel],'EdgeColor','b');
    %  plot(ROIoffset(2)+SchwerpX, ROIoffset(1)+Schwerp_Y,'r+', 'MarkerSize', 5, 'LineWidth', 1);
    t=-pi:0.01:pi;
    xOff=ROIoffset(2)+Schwerp_X+dxPixel.*0.5.*cos(t);
    yOff=ROIoffset(1)+Schwerp_Y+dyPixel.*0.5.*sin(t);
    % plot(xOff,yOff,'g')
    
    hold on;
    scatter(zaxe,dxPixel,50,'x','MarkerFacecolor','blue');
    
    zaxe_all(i) = zaxe;
    dxPixel_all(i) = dxPixel;
    
end
%Plotten der einzelnen Bilder
hold on;
[fitx,gofx]= fit(zaxe_all(:),dxPixel_all(:),"poly2");
fithandle = plot(fitx);
set(fithandle,'color',[0 .7 .7]);

[fitxw,gofxw] = fit(zaxe_all(:),dxPixel_all(:),'poly2','weights',1./dxPixel_all);