How to rotate image 3D

Dear All,

I was create the coding below. The input data (K.dcm) as attached data. But the final is below

What I want is like below:

%% Spect 
clc
clear all
close all
[spect map]=dicomread('K.dcm');
info = dicominfo('K.dcm');
%gp=info.SliceThickness;
spect=(squeeze(spect));%smooth3
aa=size(spect);aa=aa(3);
imshow3D(spect)
% 1. volume1 
BCKG1 =262;
MAX1 =1182; 
bckgratio1=MAX1/BCKG1;
seedR1 = 70; seedC1 = 81; seedP1 = 3;
W1 = graydiffweight(spect, seedC1, seedR1, seedP1 , 'GrayDifferenceCutoff', 1000000);
thresh1 =(5*0.00001*(bckgratio1^2))+(0.0008*bckgratio1)-0.0012;
[BW1, D1] = imsegfmm(W1, seedC1, seedR1, seedP1, thresh1);
% figure, imshow3D(BW1)
T1 = regionprops('table', BW1,'Area','Centroid')
figure
imshow3D(BW1);
% 2. volume2 
BCKG2 =181 ;
MAX2 =566; 
bckgratio2=MAX2/BCKG2;
seedR2 = 74; seedC2 = 73; seedP2 = 3;
W2 = graydiffweight(spect, seedC2, seedR2, seedP2 , 'GrayDifferenceCutoff', 1000000);
thresh2 =(5*0.00001*(bckgratio2^2))+(0.0008*bckgratio2)-0.0012;
[BW2, D2] = imsegfmm(W2, seedC2, seedR2, seedP2, thresh2);
% figure, imshow3D(BW1)
T2 = regionprops('table', BW2,'Area','Centroid')
figure
imshow3D(BW2);
% 3. volume3 
BCKG3 =181 ;
MAX3 =557; 
bckgratio3=MAX3/BCKG3;
seedR3 = 68; seedC3 = 69; seedP3 = 3;
W3 = graydiffweight(spect, seedC3, seedR3, seedP3 , 'GrayDifferenceCutoff', 1000000);
thresh3 =(5*0.00001*(bckgratio3^2))+(0.0008*bckgratio3)-0.0012;
[BW3, D3] = imsegfmm(W3, seedC3, seedR3, seedP3, thresh3);
% figure, imshow3D(BW1)
T3 = regionprops('table', BW3,'Area','Centroid')
figure
imshow3D(BW3);
% 4. volume4 
BCKG4 =181 ;
MAX4 =536; 
bckgratio4=MAX4/BCKG4;
seedR4 = 71; seedC4 = 59; seedP4 = 3;
W4 = graydiffweight(spect, seedC4, seedR4, seedP4 , 'GrayDifferenceCutoff', 1000000);
thresh4 =(5*0.00001*(bckgratio4^2))+(0.0008*bckgratio4)-0.0012;
[BW4, D4] = imsegfmm(W4, seedC4, seedR4, seedP4, thresh4);
% figure, imshow3D(BW1)
T4 = regionprops('table', BW4,'Area','Centroid')
figure
imshow3D(BW4);
% 5. volume5 
BCKG5 =181 ;
MAX5 =504; 
bckgratio5=MAX5/BCKG5;
seedR5 = 81; seedC5 = 54; seedP5 = 3;
W5 = graydiffweight(spect, seedC5, seedR5, seedP5 , 'GrayDifferenceCutoff', 1000000);
thresh5 =(5*0.00001*(bckgratio5^2))+(0.0008*bckgratio5)-0.0012;
[BW5, D5] = imsegfmm(W5, seedC5, seedR5, seedP5, thresh5);
% figure, imshow3D(BW1)
T5 = regionprops('table', BW5,'Area','Centroid')
figure
imshow3D(BW5);
% 6. volume6 
BCKG6 =294 ;
MAX6 =520; 
bckgratio6=MAX6/BCKG6;
seedR6 = 77; seedC6 = 44; seedP6 = 3;
W6 = graydiffweight(spect, seedC6, seedR6, seedP6 , 'GrayDifferenceCutoff', 1000000);
thresh6 =(5*0.00001*(bckgratio6^2))+(0.0008*bckgratio6)-0.0012;
[BW6, D6] = imsegfmm(W6, seedC6, seedR6, seedP6, thresh6);
% figure, imshow3D(BW1)
T6 = regionprops('table', BW6,'Area','Centroid')
figure
imshow3D(BW6);
% 7. volume7 (Salivary Gland) 
BCKG7 =181 ;
MAX7 =1380; 
bckgratio7=MAX7/BCKG7;
seedR7 = 50; seedC7 = 67; seedP7 = 65;
W7 = graydiffweight(spect, seedC7, seedR7, seedP7 , 'GrayDifferenceCutoff', 1000000);
thresh7 =(5*0.00001*(bckgratio7^2))+(0.0008*bckgratio7)-0.0012;
[BW7, D7] = imsegfmm(W7, seedC7, seedR7, seedP7, thresh7);
% figure, imshow3D(BW1)
T7 = regionprops('table', BW7,'Area','Centroid')
figure
imshow3D(BW7);
% 8. volume8 (Salivary Gland) 
BCKG8 =181 ;
MAX8 =1240; 
bckgratio8=MAX8/BCKG8;
seedR8 = 47; seedC8 = 73; seedP8 = 64;
W8 = graydiffweight(spect, seedC8, seedR8, seedP8 , 'GrayDifferenceCutoff', 1000000);
thresh8 =(5*0.00001*(bckgratio8^2))+(0.0008*bckgratio8)-0.0012;
[BW8, D8] = imsegfmm(W8, seedC8, seedR8, seedP8, thresh8);
% figure, imshow3D(BW1)
T8 = regionprops('table', BW8,'Area','Centroid')
figure
imshow3D(BW8);
allBW = BW1 | BW2 | BW3 | BW4 | BW5 | BW6 | BW7 | BW8;
imshow3D(allBW);


totalcountseachblob = sum(spect(allBW==1))
% Get a list of all files in the folder with the desired file name pattern.
myFolder = ('C:\Users\User\Desktop\New folder\K');
filePattern = fullfile(myFolder, '*.dcm'); % Change to whatever pattern you need.
theFiles = dir(filePattern);
for L = 1 : length(theFiles)
    baseFileName = theFiles(L).name;
    fullFileName = fullfile(theFiles(L).folder, baseFileName);
    fprintf(1, 'Now reading %s\n', fullFileName);
    % Now do whatever you want with this file name,
    % such as reading it in as an image array with imread()
    RZ(:,:,L) = dicomread(fullFileName);
end
RZ=double(RZ);
[N M L]=size(RZ);
rats=size(RZ)./size(spect);
[x,y,z]=meshgrid(rats(1):rats(1):M,rats(2):rats(2):N,rats(3):rats(3):L);
RZ=interp3(RZ,x,y,z);

e=10000;%isovalue, air at 0 (Faceskin800, air1250)
figure
whitebg('black')
colordef none
axis off
pause
fprintf('\nFull Reconstruction. Please wait...\n');

q = patch(isosurface(RZ,e));
axis equal
set(q,'FaceColor','w','EdgeColor','none');
alpha(q,0.2)
material shiny 
hold

camlight('headlight')
%               title([num2str(e)]);
view([-22 -4])

p = patch(isosurface(allBW));
axis equal
set(p,'FaceColor','c','EdgeColor','none');
alpha(p,1)
hold

[Cris LaPierre]

If you zip your dcm file, you can attach it to your post using the paperclip icon. We could then provide more specific instructions.

In the meantime, you may consider using medicalvolume to load your dicom instead of dicomread. It will take into consideration the orientation of the acquired data.

The MathWorks released the Medical Imaging Toolbox last year and then this year launched a Medical Image Processing course on Coursera. You may find that helpful, and you can enroll for free.

[mohd akmal masud]

Dear @Cris LaPierre,

here is my data as attached data

[mohd akmal masud]

Dear @Cris LaPierre, I solved. thank you.

[Adam Danz]

I noticed your use of colordef and whitebg.

These functions are removed in an upcoming release. Instead, set your figure color to black and your axes color to "none". Alternatively, set the figure's theme to "dark" if you're using the beta desktop of a future release.

Answer 1

[Cris LaPierre]

I noticed you must be using R2024a instead of R2024b. Some of the functionality this code uses is no longer available in R2024b, as Adam pointed out.

The two volumes you are looking at do not have the same coordinate systems. In particular, the coordinate information collected in your spect image is incorrect (green). This is what is making it hard to align the two images. You are manually trying to align them, but you shouldn't have to, as the meta data is supposed to give you the real-world transformation.

The volumes are collected using different plane mappings.

• CT mapping is sagittal, coronal, transverse
• SPECT mapping is coronal, sagittal, transverse

However, as best I can tell, the XY plane mapping for your SPECT data is swapped, causing a rotation between the two images. In addition, the origins in the two data sets are very different. When viewed using volshow, the scaling is handled automatically for you, but because the voxel spacing is different, plotting the raw data yourself will result in volumes with different sizes.

You were able to get a solution, but figured I'd share how you could do this if the meta data were correct.

% Correct the meta data
% load SPECT and CT volumes
volNM = medicalVolume('K.dcm')
volNM.Voxels = permute(volNM.Voxels,[2 1 3]); % this organizes the data with the plane mapping
volCT = medicalVolume('./K')

% Register the two volumes
[optimizer,metric] = imregconfig("multimodal");
optimizer.InitialRadius = 0.001;
RCT = imref3d(volCT.VolumeGeometry.VolumeSize,volCT.VoxelSpacing(1),volCT.VoxelSpacing(2),volCT.VoxelSpacing(3));
RNM = imref3d(volNM.VolumeGeometry.VolumeSize,volNM.VoxelSpacing(1),volNM.VoxelSpacing(2),volNM.VoxelSpacing(3));
[regVoxels,newT] = imregister(volCT.Voxels, RCT, volNM.Voxels, RNM,"translation", optimizer, metric);

% create a new volume for CT data in SPECT spacing
newCT = volNM;
newCT.Voxels = permute(regVoxels,[2,1,3]); % this changes CT data to match SPECT plane mappipng

Now run the original code (I made some minor changes)

spect = volNM.Voxels;
aa=size(spect,3);
imshow3D(spect)

% 1. volume1 
BCKG1 =262;
MAX1 =1182; 
bckgratio1=MAX1/BCKG1;
seedR1 = 70; seedC1 = 81; seedP1 = 3;
W1 = graydiffweight(spect, seedC1, seedR1, seedP1 , 'GrayDifferenceCutoff', 1000000);
thresh1 =(5*0.00001*(bckgratio1^2))+(0.0008*bckgratio1)-0.0012;
[BW1, D1] = imsegfmm(W1, seedC1, seedR1, seedP1, thresh1);
% figure, imshow3D(BW1)
T1 = regionprops('table', BW1,'Area','Centroid')
figure
imshow3D(BW1);

% 2. volume2 
BCKG2 =181 ;
MAX2 =566; 
bckgratio2=MAX2/BCKG2;
seedR2 = 74; seedC2 = 73; seedP2 = 3;
W2 = graydiffweight(spect, seedC2, seedR2, seedP2 , 'GrayDifferenceCutoff', 1000000);
thresh2 =(5*0.00001*(bckgratio2^2))+(0.0008*bckgratio2)-0.0012;
[BW2, D2] = imsegfmm(W2, seedC2, seedR2, seedP2, thresh2);
% figure, imshow3D(BW1)
T2 = regionprops('table', BW2,'Area','Centroid')
figure
imshow3D(BW2);

% 3. volume3 
BCKG3 =181 ;
MAX3 =557; 
bckgratio3=MAX3/BCKG3;
seedR3 = 68; seedC3 = 69; seedP3 = 3;
W3 = graydiffweight(spect, seedC3, seedR3, seedP3 , 'GrayDifferenceCutoff', 1000000);
thresh3 =(5*0.00001*(bckgratio3^2))+(0.0008*bckgratio3)-0.0012;
[BW3, D3] = imsegfmm(W3, seedC3, seedR3, seedP3, thresh3);
% figure, imshow3D(BW1)
T3 = regionprops('table', BW3,'Area','Centroid')
figure
imshow3D(BW3);

% 4. volume4 
BCKG4 =181 ;
MAX4 =536; 
bckgratio4=MAX4/BCKG4;
seedR4 = 71; seedC4 = 59; seedP4 = 3;
W4 = graydiffweight(spect, seedC4, seedR4, seedP4 , 'GrayDifferenceCutoff', 1000000);
thresh4 =(5*0.00001*(bckgratio4^2))+(0.0008*bckgratio4)-0.0012;
[BW4, D4] = imsegfmm(W4, seedC4, seedR4, seedP4, thresh4);
% figure, imshow3D(BW1)
T4 = regionprops('table', BW4,'Area','Centroid')
figure
imshow3D(BW4);

% 5. volume5 
BCKG5 =181 ;
MAX5 =504; 
bckgratio5=MAX5/BCKG5;
seedR5 = 81; seedC5 = 54; seedP5 = 3;
W5 = graydiffweight(spect, seedC5, seedR5, seedP5 , 'GrayDifferenceCutoff', 1000000);
thresh5 =(5*0.00001*(bckgratio5^2))+(0.0008*bckgratio5)-0.0012;
[BW5, D5] = imsegfmm(W5, seedC5, seedR5, seedP5, thresh5);
% figure, imshow3D(BW1)
T5 = regionprops('table', BW5,'Area','Centroid')
figure
imshow3D(BW5);

% 6. volume6 
BCKG6 =294 ;
MAX6 =520; 
bckgratio6=MAX6/BCKG6;
seedR6 = 77; seedC6 = 44; seedP6 = 3;
W6 = graydiffweight(spect, seedC6, seedR6, seedP6 , 'GrayDifferenceCutoff', 1000000);
thresh6 =(5*0.00001*(bckgratio6^2))+(0.0008*bckgratio6)-0.0012;
[BW6, D6] = imsegfmm(W6, seedC6, seedR6, seedP6, thresh6);
% figure, imshow3D(BW1)
T6 = regionprops('table', BW6,'Area','Centroid')
figure
imshow3D(BW6);

% 7. volume7 (Salivary Gland) 
BCKG7 =181 ;
MAX7 =1380; 
bckgratio7=MAX7/BCKG7;
seedR7 = 50; seedC7 = 67; seedP7 = 65;
W7 = graydiffweight(spect, seedC7, seedR7, seedP7 , 'GrayDifferenceCutoff', 1000000);
thresh7 =(5*0.00001*(bckgratio7^2))+(0.0008*bckgratio7)-0.0012;
[BW7, D7] = imsegfmm(W7, seedC7, seedR7, seedP7, thresh7);
% figure, imshow3D(BW1)
T7 = regionprops('table', BW7,'Area','Centroid')
figure
imshow3D(BW7);

% 8. volume8 (Salivary Gland) 
BCKG8 =181 ;
MAX8 =1240; 
bckgratio8=MAX8/BCKG8;
seedR8 = 47; seedC8 = 73; seedP8 = 64;
W8 = graydiffweight(spect, seedC8, seedR8, seedP8 , 'GrayDifferenceCutoff', 1000000);
thresh8 =(5*0.00001*(bckgratio8^2))+(0.0008*bckgratio8)-0.0012;
[BW8, D8] = imsegfmm(W8, seedC8, seedR8, seedP8, thresh8);
% figure, imshow3D(BW1)
T8 = regionprops('table', BW8,'Area','Centroid')
figure
imshow3D(BW8);
allBW = BW1 | BW2 | BW3 | BW4 | BW5 | BW6 | BW7 | BW8;
imshow3D(allBW);

totalcountseachblob = sum(spect(allBW==1))

RZ=double(newCT.Voxels);

% e=10000;%isovalue, air at 0 (Faceskin800, air1250)
e = 700
figure
whitebg('black')
colordef none
axis off
% pause
fprintf('\nFull Reconstruction. Please wait...\n');

q = patch(isosurface(RZ,e));
axis equal
set(q,'FaceColor','w','EdgeColor','none');
alpha(q,0.2)
material shiny 
hold on

camlight('headlight')
%               title([num2str(e)]);
view([-22 -4])

p = patch(isosurface(allBW));
axis equal
set(p,'FaceColor','c','EdgeColor','none');
alpha(p,1)
hold off

The registration wasn't perfect, but for a first pass, the results are decent.

[Cris LaPierre]

I made another attempt this morning using a different approach. Once the orientation is corrected, aligning the images just involves translations. I made an attempt to do this manually while still keeping the data as medical volumes. Here's what that code looked like.

websave('./NewFolder.zip','https://drive.google.com/uc?export=download&id=1WFrqURk0mEaroiqLCzgoecItJWLZ_ljn');
unzip('NewFolder.zip')

volNM = medicalVolume('New folder/K.dcm')

Warning: ImagePositionPatient attribute not detected in file metadata. Consider using medicalImage instead of medicalVolume.

Warning: ImageOrientationPatient attribute not detected in file metadata. Consider using medicalImage instead of medicalVolume.

volNM = 
  medicalVolume with properties:

                 Voxels: [130x130x90 int16]
         VolumeGeometry: [1x1 medicalref3d]
           SpatialUnits: "mm"
            Orientation: "transverse"
           VoxelSpacing: [4.6640 4.6640 4.6640]
           NormalVector: [0 0 1]
       NumCoronalSlices: 130
      NumSagittalSlices: 130
    NumTransverseSlices: 90
           PlaneMapping: ["coronal"    "sagittal"    "transverse"]
               Modality: "NM"
          WindowCenters: 772
           WindowWidths: 1545

volNM.Voxels = permute(volNM.Voxels,[2 1 3]); % So the data matches the plane mapping
volCT = medicalVolume('./New folder/K')

volCT = 
  medicalVolume with properties:

                 Voxels: [256x256x203 single]
         VolumeGeometry: [1x1 medicalref3d]
           SpatialUnits: "mm"
            Orientation: "transverse"
           VoxelSpacing: [2 2 2]
           NormalVector: [0 0 1]
       NumCoronalSlices: 256
      NumSagittalSlices: 256
    NumTransverseSlices: 203
           PlaneMapping: ["sagittal"    "coronal"    "transverse"]
               Modality: "CT"
          WindowCenters: [203x1 double]
           WindowWidths: [203x1 double]

% create a translation vector (dx,dy,dz)
t = [-300,-277,-685];

% Create spatial referencing information
R = medicalref3d(volCT.VolumeGeometry.VolumeSize,volCT.VolumeGeometry.Position-t,volCT.VolumeGeometry.VoxelDistances);
R = orient(R,volCT.VolumeGeometry.PatientCoordinateSystem);

% Create a new CT volume using the new spatial reference info
newCT = medicalVolume(volCT.Voxels,R);

% resample the newCT volume into the SPECT coordinate system
newCT = resample(newCT,volNM.VolumeGeometry)

newCT = 
  medicalVolume with properties:

                 Voxels: [130x130x90 single]
         VolumeGeometry: [1x1 medicalref3d]
           SpatialUnits: "unknown"
            Orientation: "transverse"
           VoxelSpacing: [4.6640 4.6640 4.6640]
           NormalVector: [0 0 -1]
       NumCoronalSlices: 130
      NumSagittalSlices: 130
    NumTransverseSlices: 90
           PlaneMapping: ["coronal"    "sagittal"    "transverse"]
               Modality: "unknown"
          WindowCenters: []
           WindowWidths: []

At this point, I used the following code to inspect the alignment of the two volumes. Because the modalities generate very different images, some judgement is needed to decide what is 'good enough'.

v = viewer3d();
volshow(newCT,"Colormap",[1 0 1],"RenderingStyle","Isosurface","IsosurfaceValue",0.3,Parent=v)
volshow(volNM,"Colormap",[0 1 0],Parent=v)

I then ran the same code I posted above for processing the images (your code with a few modifications).

spect = volNM.Voxels;
aa=size(spect,3);
% imshow3D(spect)

% 1. volume1 
BCKG1 =262;
MAX1 =1182; 
bckgratio1=MAX1/BCKG1;
seedR1 = 70; seedC1 = 81; seedP1 = 3;
W1 = graydiffweight(spect, seedC1, seedR1, seedP1 , 'GrayDifferenceCutoff', 1000000);
thresh1 =(5*0.00001*(bckgratio1^2))+(0.0008*bckgratio1)-0.0012;
[BW1, D1] = imsegfmm(W1, seedC1, seedR1, seedP1, thresh1);
% figure, imshow3D(BW1)
T1 = regionprops('table', BW1,'Area','Centroid')

T1 = 1x2 table
    Area             Centroid         
    ____    __________________________

     64     79.984    70.328    3.0469

% figure
% imshow3D(BW1);

% 2. volume2 
BCKG2 =181 ;
MAX2 =566; 
bckgratio2=MAX2/BCKG2;
seedR2 = 74; seedC2 = 73; seedP2 = 3;
W2 = graydiffweight(spect, seedC2, seedR2, seedP2 , 'GrayDifferenceCutoff', 1000000);
thresh2 =(5*0.00001*(bckgratio2^2))+(0.0008*bckgratio2)-0.0012;
[BW2, D2] = imsegfmm(W2, seedC2, seedR2, seedP2, thresh2);
% figure, imshow3D(BW1)
T2 = regionprops('table', BW2,'Area','Centroid')

T2 = 1x2 table
    Area             Centroid         
    ____    __________________________

     40     73.175      73.5       3.2

% figure
% imshow3D(BW2);

% 3. volume3 
BCKG3 =181 ;
MAX3 =557; 
bckgratio3=MAX3/BCKG3;
seedR3 = 68; seedC3 = 69; seedP3 = 3;
W3 = graydiffweight(spect, seedC3, seedR3, seedP3 , 'GrayDifferenceCutoff', 1000000);
thresh3 =(5*0.00001*(bckgratio3^2))+(0.0008*bckgratio3)-0.0012;
[BW3, D3] = imsegfmm(W3, seedC3, seedR3, seedP3, thresh3);
% figure, imshow3D(BW1)
T3 = regionprops('table', BW3,'Area','Centroid')

T3 = 1x2 table
    Area             Centroid         
    ____    __________________________

     31     69.323    68.516    3.3548

% figure
% imshow3D(BW3);

% 4. volume4 
BCKG4 =181 ;
MAX4 =536; 
bckgratio4=MAX4/BCKG4;
seedR4 = 71; seedC4 = 59; seedP4 = 3;
W4 = graydiffweight(spect, seedC4, seedR4, seedP4 , 'GrayDifferenceCutoff', 1000000);
thresh4 =(5*0.00001*(bckgratio4^2))+(0.0008*bckgratio4)-0.0012;
[BW4, D4] = imsegfmm(W4, seedC4, seedR4, seedP4, thresh4);
% figure, imshow3D(BW1)
T4 = regionprops('table', BW4,'Area','Centroid')

T4 = 1x2 table
    Area             Centroid         
    ____    __________________________

     79     59.139    71.253    4.3671

% figure
% imshow3D(BW4);

% 5. volume5 
BCKG5 =181 ;
MAX5 =504; 
bckgratio5=MAX5/BCKG5;
seedR5 = 81; seedC5 = 54; seedP5 = 3;
W5 = graydiffweight(spect, seedC5, seedR5, seedP5 , 'GrayDifferenceCutoff', 1000000);
thresh5 =(5*0.00001*(bckgratio5^2))+(0.0008*bckgratio5)-0.0012;
[BW5, D5] = imsegfmm(W5, seedC5, seedR5, seedP5, thresh5);
% figure, imshow3D(BW1)
T5 = regionprops('table', BW5,'Area','Centroid')

T5 = 1x2 table
    Area             Centroid         
    ____    __________________________

     27     53.481    80.556     2.963

% figure
% imshow3D(BW5);

% 6. volume6 
BCKG6 =294 ;
MAX6 =520; 
bckgratio6=MAX6/BCKG6;
seedR6 = 77; seedC6 = 44; seedP6 = 3;
W6 = graydiffweight(spect, seedC6, seedR6, seedP6 , 'GrayDifferenceCutoff', 1000000);
thresh6 =(5*0.00001*(bckgratio6^2))+(0.0008*bckgratio6)-0.0012;
[BW6, D6] = imsegfmm(W6, seedC6, seedR6, seedP6, thresh6);
% figure, imshow3D(BW1)
T6 = regionprops('table', BW6,'Area','Centroid')

T6 = 1x2 table
    Area          Centroid      
    ____    ____________________

     5      43.6    77.2       3

% figure
% imshow3D(BW6);

% 7. volume7 (Salivary Gland) 
BCKG7 =181 ;
MAX7 =1380; 
bckgratio7=MAX7/BCKG7;
seedR7 = 50; seedC7 = 67; seedP7 = 65;
W7 = graydiffweight(spect, seedC7, seedR7, seedP7 , 'GrayDifferenceCutoff', 1000000);
thresh7 =(5*0.00001*(bckgratio7^2))+(0.0008*bckgratio7)-0.0012;
[BW7, D7] = imsegfmm(W7, seedC7, seedR7, seedP7, thresh7);
% figure, imshow3D(BW1)
T7 = regionprops('table', BW7,'Area','Centroid')

T7 = 1x2 table
    Area             Centroid         
    ____    __________________________

    211     67.531    49.995    65.223

% figure
% imshow3D(BW7);

% 8. volume8 (Salivary Gland) 
BCKG8 =181 ;
MAX8 =1240; 
bckgratio8=MAX8/BCKG8;
seedR8 = 47; seedC8 = 73; seedP8 = 64;
W8 = graydiffweight(spect, seedC8, seedR8, seedP8 , 'GrayDifferenceCutoff', 1000000);
thresh8 =(5*0.00001*(bckgratio8^2))+(0.0008*bckgratio8)-0.0012;
[BW8, D8] = imsegfmm(W8, seedC8, seedR8, seedP8, thresh8);
% figure, imshow3D(BW1)
T8 = regionprops('table', BW8,'Area','Centroid')

T8 = 1x2 table
    Area             Centroid         
    ____    __________________________

    213     72.418     48.31    64.221

% figure
% imshow3D(BW8);
allBW = BW1 | BW2 | BW3 | BW4 | BW5 | BW6 | BW7 | BW8;
% imshow3D(allBW);

totalcountseachblob = sum(spect(allBW==1))

totalcountseachblob = 372964

RZ=double(newCT.Voxels);

% e=10000;%isovalue, air at 0 (Faceskin800, air1250)
e = 700;
figure
whitebg('black')

Warning: whitebg will be removed in a future release.

colordef none

Warning: colordef will be removed in a future release.

axis off
% pause
fprintf('\nFull Reconstruction. Please wait...\n');

Full Reconstruction. Please wait...

q = patch(isosurface(RZ,e));
axis equal
set(q,'FaceColor','w','EdgeColor','none');
alpha(q,0.2)
material shiny 
hold on

camlight('headlight')
%               title([num2str(e)]);
view([-22 -4])

p = patch(isosurface(allBW));
axis equal
set(p,'FaceColor','c','EdgeColor','none');
alpha(p,1)
hold off

Here's some code for creating a 3D visulaization of the results using volshow.

% CT Lung visualization parameters

alphamap = (1:256>=60).*(1:256<=80)*0.15;

v = viewer3d('BackgroundColor','w');

volshow(RZ,Alphamap=alphamap,Parent=v);

volshow(allBW,"Colormap",[0 1 0],Parent=v)