Average of multiple matrices
Show older comments
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
Accepted Answer
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)
% subtract the mean from each slice
for i = 1:size(IM,3)
IM(:,:,i) = IM(:,:,i) - MyMean;
end
15 Comments
Latifa Bouguessaa
on 9 Jul 2022
Image Analyst
on 9 Jul 2022
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
on 9 Jul 2022
Walter Roberson
on 10 Jul 2022
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
on 10 Jul 2022
Karim
on 10 Jul 2022
@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
on 10 Jul 2022
Latifa Bouguessaa
on 10 Jul 2022
Walter Roberson
on 10 Jul 2022
zip the image and upload the zip
Latifa Bouguessaa
on 10 Jul 2022
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 )
% 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
on 11 Jul 2022
Karim
on 11 Jul 2022
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
on 11 Jul 2022
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
on 12 Jul 2022
More Answers (0)
Categories
Find more on Images in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
