numerical integration using trap()

Hi.

I am trying to solve the functions below numerically using MATLAB

this is how I did it using trapz()

%1.Moment
[row,column] = size(filtered);
Summe = 0;
for y=1:row
    for x=1:column
        trapz(trapz(x.*double(filtered(y,x))));
    end
end
first_Moment_X = Summe;
Summe = 0;
Schwerp_X = first_Moment_X./en_density;
Schwerp_X = double(Schwerp_X);
for y=1:row
    for x=1:column
        trapz(trapz(y.*double(filtered(y,x))));
    end
end
first_Moment_Y=Summe;
Schwerp_Y = first_Moment_Y./en_density;

But I get zero as result. Does anyone know what I do wrong?

Answer 1

[Star Strider]

Since ‘filtered’ is apparently a matrix, perhaps —

filtered = randn(10,10);

%1.Moment
[row,column] = size(filtered);
x=1:column;
Summe = trapz(trapz(x.*double(filtered)));
en_density = trapz(trapz(filtered));
first_Moment_X = Summe;
Summe = 0;
Schwerp_X = first_Moment_X./en_density;
Schwerp_X = double(Schwerp_X)

Schwerp_X = 4.1243

y=1:row;
Summe = trapz(trapz(y.*double(filtered)));
first_Moment_Y=Summe;
Schwerp_Y = first_Moment_Y./en_density

Schwerp_Y = 4.1243

This will produce a non-zero result (unless ‘filtered’ is such that the double integral is zero) and unless it is a square matrix, the multiplication of ‘x’ or ‘y’ may need to be transposed so that the multiplication is conformable to the matrix dimensions. (I have no idea what ‘en_density’ is, so I created it by doing a couble integral of ‘filtered’.)

Make appropriate changes to get the desired result.

The zero result appears to be the integration of individual points rather than of the entire matrix, for example —

Q = trapz(trapz(rand))

Q = 0

.

[Bersin Tekmen]

Thanks for your answer @Star Strider. I see you don't use a for-loop anymore? You're right ‘filtered’ is a image. So I read the pixels and get a matrix. But dont I need a for-loop to iterate through the matrix?

%1.Moment
[row,column] = size(filtered);
Summe = 0;
for y=1:row
    trapz(trapz(y.*double(filtered)));
    for x=1:column
        trapz(trapz(x.*double(filtered)));
    end
end
first_Moment_X = Summe;
Summe = 0;
Schwerp_X = first_Moment_X./en_density;
Schwerp_X = double(Schwerp_X);
first_Moment_Y=Summe;
Schwerp_Y = first_Moment_Y./en_density;
Schwerp_Y = double(Schwerp_Y);

This is now how I am trying it.

But it takes long time to solve rn...

[Star Strider]

‘But dont I need a for-loop to iterate through the matrix?’

I did not realise that ‘filtered’ is an image. The loop depends on what you are integrating. If there are three channels, then a separate loop for each channel could be necessary, however the trapz calls would be the same otherwise in my code. Integrating individual points will always produce in a zero result.

It will also be necessary to save the results of the trapz calls to a variable (either a scalar or a vector), with or without a loop. If they are not specifically saved, they will get overwritten with the next result that is not saved to a specific variable.

Experiment with your image with my code to see if the result is appropriate.

[Bersin Tekmen]

like that? @Star Strider

%1.Moment
[row,column] = size(filtered);

for y=1:row
    first_Moment_Y = trapz(trapz(y(i,1).*double(filtered)));
    for x=1:column
        first_Moment_X = trapz(trapz(x(i,1).*double(filtered)));
    end
end
Schwerp_X = first_Moment_X./en_density;
Schwerp_X = double(Schwerp_X);
Schwerp_Y = first_Moment_Y./en_density;
Schwerp_Y = double(Schwerp_Y);

[Torsten]

1 first_Moment_X and first_Moment_Y are always overwritten in the loop

2 the loop variables are the x and y coordinates ? I doubt it.

[Star Strider]

@Bersin Tekmen — I would not use the loops. Do something similar to what I wrote.

[Bersin Tekmen]

@Star Strider with your code I get: Arrays have incompatible sizes for this operation. And yes you were right "en_density" is calculated like you did. It is the energy density of the laser.

[Bersin Tekmen]

@Torsten i tried the code below. But the code is running like 2h meanwhile... I just stopped it...

%1.Moment
[row,column] = size(filtered);

for y=1:row
    first_Moment_Y(i,1) = trapz(trapz(y.*double(filtered)));
    for x=1:column
        first_Moment_X(i,1) = trapz(trapz(x.*double(filtered)));
    end
end
Schwerp_X = first_Moment_X./en_density;
Schwerp_X = double(Schwerp_X);
Schwerp_Y = first_Moment_Y./en_density;
Schwerp_Y = double(Schwerp_Y);

'2 the loop variables are the x and y coordinates ? I doubt it.'

What can it be? Can you make any suggestion? x and y are the rows and columns of the image I have to integrate though the whole pic.. and there are 22 pics with 1024x1280 pixels...

[Star Strider]

@Bersin Tekmen — The array incompatibility may relate to the sizes of ‘x’ and ‘y’ with respect to the matrix. I have no idea what they are, so I cannot suggest a specific way to do the multiplications. (In my code, the matrix is square, so everything works.)

[Bersin Tekmen]

@Star Strider thats the whole code. Does it help you to understand it?

Here you can see all declerations and variables:

close all;
clear variables;
%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");
name(strcmp(name,"460_dunkel.bmp")) = []; % remove 460_dunkel.bmp
lambda = 0.633; %Wellenlänge Laser in um
Pixelgroesse = 4.65; %Pixelgroesse in um
for i = 1:22
    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];
    
    en_density = trapz(trapz(filtered));

This part now is for the formula i showed above

    %1.Moment
    [row,column] = size(filtered);
    
    for y=1:row
        first_Moment_Y(i,1) = trapz(trapz(y.*double(filtered)));
        for x=1:column
            first_Moment_X(i,1) = trapz(trapz(x.*double(filtered)));
        end
    end
    Schwerp_X = first_Moment_X./en_density;
    Schwerp_X = double(Schwerp_X);
    Schwerp_Y = first_Moment_Y./en_density;
    Schwerp_Y = double(Schwerp_Y);

Then there is the part for a very similar formula:

    %2.Moment
    [row,column] = size(filtered);
    Summe = 0.;
    for y=1:row
        sec_Moment_Y(i,1) = trapz(trapz(((y-Schwerp_Y).^2).*double(filtered)));
        for x=1:column
            sec_Moment_X(i,1) = trapz(trapz(((x-Schwerp_X).^2).*double(filtered)));
        end
    end
    
    Sigma_X = sec_Moment_X./en_density;
    Sigma_X = double(Sigma_X);
    Sigma_Y = sec_Moment_Y./en_density;
    Sigma_Y = double(Sigma_Y);

And then there is a plot and and a print and so on.. but I think that doesn't matter for the problem I have here.

As you can see I have an image with 1024x1280 pixels... so x = 1024 and y is 1280.

So. How can I just solve that equation with the help of matlab now?

[Star Strider]

I am lost.

I stil believe that the multiplication of the vector ‘x’ or ‘y’ by the ‘filtered’ matrix should be done using either ordinary matrix-vector multiplication or element-wise (array) matrix-vector multiplication, however I have no idea what the code actually does, so I cannot determine which would be correct. The vector-matrix multiplications must be conformable (the dimensions must match) regardless.

[Torsten]

@Bersin Tekmen

Similar to the integration for the first moments, I wonder how you want to do numerical integration if the results depend on a variable, namely z. Don't you need some kind of symbolic integration here ? Or do you have to perform integration for each element in z-direction separately ?