Average of multiple matrices

Hello dear forum members

I have a question, I have 36 images and I want to calculate this formula (M=Bi - mean value of all images)

that is, I have to do a matrix as an average of the 36 images and then subtract the image from individual images Bi and at the end I sum the results.

does someone have an idea, how I realize this

Your Latifa

Answer 1

[Karim]

This depends how your images are stored in memory, if they are stored in a 3D array you can use the mean function

% generate 36 matrices, stored in a 3D array
IM = rand(60,60,36);

% here '3' indicates the 3d dimension (i.e. the slices)
MyMean = mean( IM, 3)

MyMean = 60×60

0.5423 0.5302 0.4499 0.6049 0.5529 0.4937 0.4482 0.5046 0.4249 0.5116 0.5654 0.5206 0.5490 0.4849 0.4883 0.4870 0.4785 0.5159 0.4727 0.4874 0.4605 0.5692 0.5124 0.5299 0.3905 0.5232 0.5547 0.4935 0.5231 0.4855 0.5468 0.4794 0.4764 0.5231 0.5032 0.5078 0.4658 0.4594 0.4497 0.5204 0.5166 0.4808 0.5565 0.4713 0.4848 0.4805 0.4335 0.5317 0.5321 0.4965 0.4896 0.5081 0.4681 0.3973 0.4686 0.4634 0.5425 0.5439 0.5559 0.4928 0.5300 0.5006 0.5265 0.4788 0.5069 0.5307 0.4282 0.5432 0.4681 0.4531 0.5471 0.5047 0.5714 0.4917 0.5110 0.5373 0.5117 0.4247 0.4278 0.4275 0.5413 0.5046 0.5123 0.4798 0.4597 0.5371 0.5881 0.4942 0.5025 0.5276 0.4752 0.4677 0.4788 0.4446 0.5567 0.4630 0.4647 0.5636 0.5216 0.5469 0.5187 0.4974 0.5457 0.5155 0.4299 0.5183 0.3784 0.4394 0.4890 0.4367 0.5318 0.4033 0.4858 0.5446 0.4996 0.4741 0.5789 0.4643 0.4695 0.5074 0.4532 0.3975 0.6091 0.4424 0.5648 0.4693 0.4671 0.4325 0.5529 0.4172 0.5023 0.4736 0.5148 0.5017 0.5227 0.5115 0.5303 0.4667 0.5141 0.5075 0.4537 0.4380 0.4438 0.5337 0.4947 0.5352 0.4409 0.4874 0.5363 0.5250 0.4288 0.5157 0.4929 0.4885 0.4616 0.5684 0.5134 0.5041 0.4447 0.4677 0.4846 0.4331 0.4827 0.5668 0.4760 0.4696 0.4597 0.5153 0.4943 0.5200 0.5452 0.5485 0.5282 0.5239 0.5404 0.4633 0.4966 0.4686 0.5248 0.4939 0.4564 0.4663 0.5359 0.5494 0.4919 0.5289 0.6235 0.4576 0.4335 0.5506 0.5256 0.4485 0.5047 0.4727 0.5368 0.4898 0.3592 0.4526 0.4088 0.5358 0.4045 0.3875 0.4460 0.4536 0.5408 0.4599 0.4638 0.5069 0.5380 0.5099 0.6138 0.3802 0.5928 0.4648 0.5081 0.4707 0.4212 0.4195 0.5440 0.4571 0.5042 0.5226 0.5296 0.5771 0.3977 0.5137 0.5115 0.4058 0.5647 0.4902 0.5190 0.5141 0.4769 0.5969 0.4781 0.4203 0.5499 0.4942 0.4371 0.4919 0.6080 0.5764 0.4546 0.4690 0.4312 0.5565 0.4930 0.4178 0.4755 0.4580 0.4496 0.4649 0.4300 0.4207 0.5368 0.4766 0.5311 0.4218 0.5017 0.4895 0.5588 0.5115 0.4546 0.5221 0.5602 0.5118 0.4663 0.5374 0.4802 0.4736 0.4812 0.5693 0.4882 0.5565 0.4741 0.5420 0.5190 0.5387 0.4813 0.5916 0.5528 0.5074 0.4795 0.5309 0.4927 0.5546 0.4946 0.5187 0.5807 0.5787 0.5524 0.5082 0.5658 0.5423 0.4987 0.6166 0.4134 0.5001 0.4329 0.4975

% subtract the mean from each slice
for i = 1:size(IM,3)
    IM(:,:,i) = IM(:,:,i) - MyMean;
end

[Latifa Bouguessaa]

hello Karim thanks for your help

I'm trying to extract the mean of individual matrices in a loop but unfortunately it doesn't work. here is my code:

--------------------------------------

for i=1:size(subMat,3)

M = (subMat(:,:,36)- mean(subMat,3));

end

-----------------------------------

Error using -

Integers can only be combined with integers of the same class, or scalar doubles.

Error in test1 (line 150)

M = (subMat(:,:,36)- mean(subMat,3));

--------------------------------------------------------------

I hope you have a solution

kind regards Latifa

[Image Analyst]

Are your 36 images on disk? Or do you already have them as 36 2-D or 36 3-D/color images, or as a 36 slice 3-D image with each slice being a gray scale image?

[Latifa Bouguessaa]

hello the images are in gray values 60X60 consist of several ROIs of a CT image

[Walter Roberson]

@Latifa Bouguessaa

Could you confirm that what you have is a single image, and that the 36 images of interest are portions of the one large image? If so then what information do you have about the position of the smaller images within the larger image?

[Latifa Bouguessaa]

yes I can confirm that. the small pictures are part of a larger picture. angle between radius between center of roi and angle are information about the position

Best wishes latifa

[Karim]

@Latifa Bouguessaa, i updated the answer to demonstrate the indexing for the substration. Notice the difference in indexing with the loop example in your comment. The loop will work for data stored in a 3d matrix. In this example each slice represents a 2d matrix which holds e.g. the grayscale values of an image.

However, it is not clear if this will actually help you. Do you already have the data (i.e. each image) stored as slices of a 3d matrix? Can you share some example code or data?

[Latifa Bouguessaa]

Hello Karim,

Your ROIs are sliced but I can't edit them in the loop.

I'll post my CODE here and the DICOM image that I need to process. This code is the noise power spectrum of the CT images.

Thanks for your help Karim

kind regards Latifa

[Latifa Bouguessaa]

Unfortunately, the Dicom image cannot be uploaded

[Walter Roberson]

zip the image and upload the zip

[Latifa Bouguessaa]

here is the DICOM image

[Karim]

I tried to fix the code, see below for the results. The error was due to the difference in variable types: the mean command results in 'double' type while the data in 'subImage' is of type 'int16'.

Hence I added a bit of code to convert the the result of the mean to type 'int16'. This allows the code to proceed, but I do not know if this is allowed for your computation. An alternative is to convert subImage to doubles. And look if this influences your result.

Some small remark, note that the loop where you evaluate 'spe' is not really usefull, at each iteration you will overwrite the variable. However, i noticed that a few lines later you divide by 36. Was it your intention to add all the matrices?

unzip("https://www.mathworks.com/matlabcentral/answers/uploaded_files/1060430/Testdaten_0300.zip")

Citra = dicomread('Testdaten_0300.dcm');
info = dicominfo('Testdaten_0300.dcm');
Potong=info.RescaleIntercept;
Miring=info.RescaleSlope;
FV=info.ReconstructionDiameter;
Spasi=info.PixelSpacing;
imshow(Citra,'DisplayRange',[]);

% Umwandelung der CT Daten in HU
[a b]=size(Citra);
D1=int16(Citra);
D2=zeros(a,b);
D2=int16(D2);
for k=1:a
    for l=1:b
        D2(k,l)=D1(k,l)*Miring+Potong;
    end
end
figure(1)
imshow(D2,'DisplayRange',[]);

% Automatische Segmentierung
B=zeros(a,b);
for i=1:a
    for j=1:b
        if D2(i,j) < (800-1024)
            B(i,j)=0;
        else
            B(i,j)=1;
        end
    end
end
Obyek=sum(sum(B));
if Obyek > 1
    BR = bwmorph(B,'remove');
    BL=bwlabel(BR);
    NM=max(max(BL));
    for put=1:NM
        CB=zeros(a,b);
        for i=1:a
            for j=1:b
                if BL(i,j) == put
                    CB(i,j)=1;
                end
            end
        end
        DB = imfill(CB,'holes');
        LA(put)=bwarea(DB);
    end
    terbesar=max(LA);
    for iter=1:NM
        if (LA(iter)==terbesar)
            urutan =iter;
        end
    end
    CB=zeros(a,b);
    for i=1:a
        for j=1:b
            if BL(i,j) == urutan
                CB(i,j)=1;
            end
        end
    end
    DB = imfill(CB,'holes');
else
    DB=zeros(a,b);
end
figure(2)
imshow(D2,'DisplayRange',[]);
DBP = bwmorph(DB,'remove');
[n m]=size(DBP);
p=0;
for i=1:n
    for j=1:m
        if DBP(i,j)==1
            p=p+1;
            x(p)=j;
            y(p)=i;
        end
    end
end
hold on
scatter(x,y,'w', '.')
% Berechnung des Zentrums des Phantoms
x0=round(a/2); y0=round(b/2);
N=0;
xpos=0; ypos=0;
for i=1:a
    for j=1:b
        if DB(i,j) > 0
            N=N+1;
            xpos=xpos+i;
            ypos=ypos+j;
        end
    end
end
yb=round(xpos/N);
xb=round(ypos/N);
plot(xb,yb,'w.', 'MarkerSize', 5)
%Festlegung des Schrittwinkels und des Radius
[rows, columns, numberOfColorChannels] = size(D2);
xCenter = xb;
yCenter = yb;
allTheta = 0 : 10 : 359;
radius = 100;
% Kreis erstellen
viscircles([xCenter, yCenter], radius, 'Color','y', 'LineWidth', 1);
boxWidth = 60;
co = 1;
for k = 1 : length(allTheta)
    % Winkel
    theta = allTheta(k);
    % Punkte auf dem Kreis
    x = radius * cosd(theta) + xCenter;
    y = radius * sind(theta) + yCenter;
    % ROIs erstellen
    xul = x - boxWidth/2;
    yul = y - boxWidth/2;
    pos = [xul, yul, boxWidth, boxWidth];
    rectangle('Position', pos, 'EdgeColor', 'g', 'LineWidth', 1)
    axis square;
    grid on;
    subImage(:,:,co) = imcrop(D2, pos);
    co=co+1;
end

figure(4)
for i=1:36
    %Darstellung der Rois
    subplot(6,6,i)
    imagesc(subImage(:,:,i))
end

%mittelwert = mean(subImage,3);
%Berechnung des NPS
%Ermittelung der Pixelabstände
Deltax=Spasi(1);
Deltay=Spasi(2);
%Berechnung des Faktors
f=(Deltax*Deltay)/(2*size(subImage,1)*size(subImage,2));

% look at the datatype op 'subImage'
class( subImage )

ans = 'int16'

% take the mean of the images and convert to int16 for matrix subtration
mymean = int16( mean( subImage , 3) );

% this step makes little sense... at each iteration you overwrite the variable 'spe'
for i = 1:size(subImage,3)
    spe = subImage(:,:,i) - mymean;
end
spe=fftshift(abs(fft(spe))).^2;
spel1=spe;
nps=f*spel1./36;
nps2=movmean(nps,40);
l=size(nps);
[a b]=size(subImage(:,:,1));
spaki=Spasi(1,1);
[f g]=size(nps2);
%Frequenz
frequency=(1/((b)*spaki));
sumbu=0:frequency:(((g-1)*frequency));
%kurve=nps(1,round(g/2):g);
figure(5)
%Ploten der NPS Kurve
plot(sumbu, 10*log(nps2), 'LineWidth',0.5)
xlabel('Spatial Frquency  [mm]^-1')
ylabel('NPS [HU]^2 [mm]^2')

set(gcf, 'color', 'w')

[Latifa Bouguessaa]

Hello Karim

First of all thanks for your help.

When calculating the NOISE Power Spectrum of the CT images, several Rois were specified in a homogeneous CT image. A 1D curve of the NPS must then be calculated from these ROIs.

The rois are stored as a slice, which then leads to errors when subtracting the mean image from each image. I'll upload the formula so you can understand what I have to program.

I hope you have an idea to solve this problem

many thanks again

Kind regards Latifa

[Karim]

In the routine above i indicated why the error occurs when you try to substract the mean (type double) from the images (type intigers): matlab can only substract matrices of similair type. After fixing a few lines the code does run and no longer leads to errors. I did run it in the answer, you can see the images in the comment.

While the code does run, I cannot comment on whenter or not it produces somthing correct or something that makes sense. As indicated, I suspect that some things are not yet correct (e.g. the compuation of 'spe') however I'm not familiair with noise power spectrums of CT images.

[Walter Roberson]

M = (subMat(:,:,36)- mean(subMat,3));

That line takes the mean over the third dimension, giving a 2D double result. You then try to subtract that 2D double result from a 2D uint8 matrix. However, when you do computations on uint8 matrices, the two operands for the operation must both be uint8 unless one of the operands is a scalar double.

An alternate way of writing the code would be

M = (subMat(:,:,36)- mean(subMat,3,'native'));

which would cause mean() of uint8 data to return a uint8 result.

Reminder: when you subtract uint8 from uint8, if the subtraction would return negative, zero is returned instead. So all values in slice 36 that were less than the mean value of subMat over the third dimension, will be transformed into 0 when you use uint8 calculations. If you want the possibility of negative values then you need to conver the uint8 to double before you do the subtraction.

[Latifa Bouguessaa]

Hello Walter Roberson, Karim and Image Analyst

thank you all very much for your help