How to extract row, column, slice for image spect

Hi all,

i have one image SPECT. When i read the information in dicominfo, it state that there is 130 rows, 130 colums, and 42 slices. But when i used below formula, it cannot give the rows,columns, slice location. Anybody can help me?

My image name is I10.dcm

Below is my code that i used

a = dicomread('I10.dcm' );
[r,c,slice] =  findND (a(:)>=800 & a(:)<=850 )  % min pixel value is 0, max pixel value is 850

% where code for findND is below

function varargout=findND(X,varargin)
%Find non-zero elements in ND-arrays. Replicates all behavior from find.
% The syntax is equivalent to the built-in find, but extended to
% multi-dimensional input.
%
% [...] = findND(X,K) returns at most the first K indices. K must be a
% positive  scalar of any type.
%
% [...] = findND(X,K,side) returns either the first K or the last K
% inidices. The input side  must be a char, either 'first' or 'last'. The
% default behavior is 'first'.
%
% [I1,I2,I3,...,In] = findND(X,...) returns indices along all the
% dimensions of X.
%
% [I1,I2,I3,...,In,V] = findND(X,...) returns indices along all the
% dimensions of X, and additionally returns a vector containg the values.
%
% Note for Matlab 6.5:
% The syntax with more than one input is present in the online doc for R14
% (Matlab 7.0), so this might be the latest release without support for
% this syntax.
%
% Compatibility:
% Matlab: should work on all releases (tested on R2017b, R2012b and 6.5)
% Octave: tested on 4.2.1
% OS:     should work cross-platform
%
% Version: 1.1
% Date:    2017-12-29
% Author:  H.J. Wisselink
% Email=  'h_j_wisselink*alumnus_utwente_nl';
% Real_email = regexprep(Email,{'*','_'},{'@','.'})
%
% Changes from v1.0 to v1.1:
% - Added support for Matlab 6.5 (R13).
% - Fixed a minor bug where the orientation of the output vector did not
%   match the orientation of the built-in find function. It also always
%   returned 1 as value, instead of the true value.
% - Added support upper case third input.

%Parse inputs
%validateattributes(X,{'numeric','logical'},{'nonempty'},1)
if ~(isnumeric(X) || islogical(X)) || numel(X)==0
    error('Expected first input (X) to be a non-empty numeric or logical array.')
end
switch nargin
    case 1
        %[...] = findND(X);
        side='first';
        K=inf;
    case 2
        %[...] = findND(X,K);
        side='first';
        K=varargin{1};
        %validateattributes(K,{'numeric','logical'},{'scalar','positive'},2)
        if ~(isnumeric(K) || islogical(K)) || numel(K)~=1 || any(K<0)
            error('Expected second input (K) to be a positive numeric or logical scalar.')
        end
    case 3
        %[...] = FIND(X,K,'first');
        K=varargin{1};
        %validateattributes(K,{'numeric','logical'},{'scalar','positive'},2)
        if ~(isnumeric(K) || islogical(K)) || numel(K)~=1 || any(K<0)
            error('Expected second input (K) to be a positive numeric or logical scalar.')
        end
        side=varargin{2};
        if ~isa(side,'char') ||...
                ~( strcmpi(side,'first') || strcmpi(side,'last'))
            error('Third input must be either ''first'' or ''last''.')
        end
        side=lower(side);
    otherwise
        error('Incorrent number of inputs.')
end

%parse outputs
nDims=length(size(X));
%allowed outputs: 0, 1, nDims, nDims+1
if nargout>1 && nargout<nDims
    error('Incorrect number of output arguments.')
end

varargout=cell(nargout,1);
v=version;v=str2double(v(1:3));
if v<7
    %The find(X,k,side) syntax was introduced between 6.5 and 7
    if nargout>nDims
        [ind,col_index_equal_to_one,val]=find(X(:));%#ok no tilde pre-R2009b
        %X(:) converts X to a column vector. Treating X(:) as a matrix
        %forces val to be the actual value, instead of the column index.
        if length(ind)>K
            if strcmp(side,'first')
                %select first K outputs
                ind=ind(1:K);
                val=val(1:K);
            else
                %select last K outputs
                ind=ind((end-K):end);
                val=val((end-K):end);
            end
        end
        [varargout{1:(end-1)}] = ind2sub(size(X),ind);
        varargout{end}=val;
    else
        ind=find(X);
        if length(ind)>K
            if strcmp(side,'first')
                %select first K outputs
                ind=ind(1:K);
            else
                %select last K outputs
                ind=ind((end-K):end);
            end
        end
        [varargout{:}] = ind2sub(size(X),ind);
    end
else
    if nargout>nDims
        %Tilde (~) to ignore outputs was introduced in R2009b. It is
        %probably faster to ignore the extra output than to use an if.
        [ind,col_index_equal_to_one,val]=find(X(:),K,side);%#ok
        %X(:) converts X to a column vector. Treating X(:) as a matrix
        %forces val to be the actual value, instead of the column index.
        [varargout{1:(end-1)}] = ind2sub(size(X),ind);
        varargout{end}=val;
    else
        ind=find(X,K,side);
        [varargout{:}] = ind2sub(size(X),ind);
    end
end

Answer 1

[Image Analyst]

With this line of code

[r,c,slice] =  findND (a(:)>=800 & a(:)<=850 )

when you do (:) that turns (the poorly named) a into a column vector, and so a(:)>=800 & a(:)<=850 is also a column vector. Thus all information about how manu rows, columns, and slices a has is lost, and findND() will not have it. So findND can only deliver the non-zero rows, NOT the columns and slices. Why not use regular find() in a loop?

for sliceIndex = 1 : size(a, 3)
    mask = a(:,:,sliceIndex) >= 800 & a(:,:,sliceIndex) <= 850;
    [rows, columns] = find(mask); % rows and columns for this slice only, not all of them.
    % Etc.  Now add these rows and columns to a master list of them over all slices.  (You do this)...
end

[mohd akmal masud]

Thanks Image Analyst,

I've tried it. I wrote the code like below.

 a = dicomread('I10.dcm' );
for sliceIndex = 1 : size(a, 3)
    mask = a(:,:,18) >= 600 & a(:,:,28) <= 850;
    [rows, columns] = find(mask) % rows and columns for this slice only, not all of them.
    % Etc.  Now add these rows and columns to a master list of them over all slices.  (You do this)...
end            

But it just appeared like this,

rows =

0×1 empty double column vector

columns =

0×1 empty double column vector

I want to extract the rows,column and the slice number that contain pixel value from 800 till 850. And hope the answer like below (I give the example)

rows =

122

103

104

columns =

121

99

100

slice=

12

18

21

[Image Analyst]

Then there is no row and column that is non-zero in BOTH slice 18 and slice 28.

Why did you inspect only those two slices in a loop over all slices?

You can make a logical map of where those values are by doing

mask = (a >= 800) & (a <= 850);

You should be able to do everything you want with that mask. You might not, and probably don't, need the individual row, column, and slice indexes.

[mohd akmal masud]

Hi Image Analyst,

it work i used this below formula

a = dicomread('I10.dcm' );
for sliceIndex = 1 : size(a, 3)
  mask = (a >= 800) & (a <= 850);
    [rows, columns, slice] = find(mask) % rows and columns for this slice only, not all of them.
    % Etc.  Now add these rows and columns to a master list of them over all slices.  (You do this)...
end   

but, the no. of rows is correct, but something wrong with no. of columns and no. of slice. the answer i got is like below

rows =

73

73

74

75

columns =

3705

3705

3835

3836

slice =

4 x 1 logical array

1

1

1

1

i think the slice should be from 1 to 42 (because the image have 42 slice),

then the columns should be 1 to 130 (because the dimension image is130 x 130)

Have you any ideas to solve it?

[Image Analyst]

No, find() does not work with 3-D arrays. When I said to make the mask over all slices I meant BEFORE the loop, like

a = dicomread('I10.dcm' );
mask = (a >= 800) & (a <= 850);
for sliceIndex = 1 : size(a, 3)
    % Get 2-D mask at this slice level.
    thisMask = mask(:,:,sliceIndex);
    % Extract rows and columns for this slice only, not all of them.
    [rows, columns] = find(thisMask); 
    % Etc.  Now add these rows and columns to a master list of them over all slices.  (You do this)...
end

but again, I don't think you need the rows and columns. Like I said before, it's possible/probably that you can do everything you need to do with the binary mask array rather than having long lists of x,y,z coordinates of the mask array.

[mohd akmal masud]

oh i see.. then do you have suggestion/idea to plot the 3D or plot voxel that have value from 500 to 850 for example. please help me

because my purpose to have list of x,y,z is to plot the voxel to make 3D as well. can you help me. Thanks in advance.

[Image Analyst]

If you want to do 3-D visualization, you're best off using a program made for that like Avizo.

MATLAB is generally limited to cut-away slices, iso surfaces, or 3-D scatterplots. Look up 3-D visualization in the help.

[mohd akmal masud]

ok Thank you so much image analyst for your explanation from the top till the end. im very appreciate it.

but if you dont mind, can you share with me how to build up iso surface just using the x and y points?

[Image Analyst]

The isosurface() function takes x, y, and z values that are output from meshgrid.

[rows, columns, slices] = size(mask);
[x, y, z] = meshgrid(1:columns, 1:rows, 1:slices);

Then you can set the value for each location. Try initializing v to all zeros, then setting mask. Untested code:

for k = 1 : numel(x)
    thisx = x(k)
    thisy = y(k)
    thisz = z(k)
    v(k) = mask(thisy, thisx, thisz);
end

then call isosurface():

isosurface(x, y, z, v, 1);

I haven't done or tested this. If that doesn't work, you can probably figure it out as well as I can.

[mohd akmal masud]

Thanks Image Analyst, i will try first

[mohd akmal masud]

Hi Image Analyst,

i was try what you give like below, it work,

The code i used like below

 a = dicomread('I10.dcm' );
mask = (a >= 800) & (a <= 850);
[rows, columns, slices] = size(mask);
[x, y, z] = meshgrid(1:columns, 1:rows, 1:slices);
for k = 1 : numel(x)
    thisx = x(k)
    thisy = y(k)
    thisz = z(k)
    v(k) = mask(thisy, thisx, thisz);
end

isosurface(x, y, z, v, 1)

But some error below

Error using isosurface (line 75)

V must be a 3D array.

Error in cubacubacuba (line 12)

isosurface(x, y, z, v, 1)

[Image Analyst]

Then try passing in mask instead of v.

[mohd akmal masud]

you mean like this?

 a = dicomread('I10.dcm' );
mask = (a >= 800) & (a <= 850);
[rows, columns, slices] = size(mask);
[x, y, z] = meshgrid(1:columns, 1:rows, 1:slices);
for k = 1 : numel(x)
    thisx = x(k)
    thisy = y(k)
    thisz = z(k)
    mask(k) = mask(thisy, thisx, thisz);
end

isosurface(x, y, z, mask, 1)

[Image Analyst]

No, like this (I think, untested)

 a = dicomread('I10.dcm' );
mask = (a >= 800) & (a <= 850);
[rows, columns, slices] = size(mask);
[x, y, z] = meshgrid(1:columns, 1:rows, 1:slices);
isosurface(x, y, z, mask, 1)

[mohd akmal masud]

its same

[mohd akmal masud]

still error.

[Image Analyst]

Well then I don't know. Sorry, but I've never used isosurface. You'll have to work it out yourself or call tech support.

[mohd akmal masud]

ok thanks Image Analyst,

so to double confirm with you, let say i want to extract the rows and columns for some slice, is it right i used the code below? let say for slice no. 30

a = dicomread('I10.dcm' );
mask = (a >= 800) & (a <= 850);
for sliceIndex = 1 : size(a, 3)
    % Get 2-D mask at this slice level.
    thisMask = mask(:,:,30);
    % Extract rows and columns for this slice only, not all of them.
    [rows, columns] = find(thisMask); 
    % Etc.  Now add these rows and columns to a master list of them over all slices.  (You do this)...
end

Answer 2

[mohd akmal masud]

Sorry all, another question i have but i wrote at this space. please help me

Dear all,

this is my code to view CT image by slice

P = zeros(256, 256, 72);
for K = 1 : 72
    petname = sprintf('I4%03d.dcm', K);
    P(:,:,K) = dicomread(petname);
end
imshow3D(P)

then, this is my code for view SPECT image by slice,

Noted: all my 512 slice SPECT image stored in one file.

[spect map]=dicomread('128x128');
info = dicominfo('128x128');
gp=info.SliceThickness;
spect=(squeeze(spect));%smooth3
aa=size(spect);aa=aa(3);
imshow3D(spect);

Anybody can help me to fuse both SPECT and CT images for all slice?