Radiomics - github code https://gi​thub.com/m​vallieres/​radiomics

Hello,

I have a problem running the following github code used for reading a DICOM image of H&N cancer from The Cancer Imaging Archive(TCIA) in the windows platform. A function named sData is declared, but in the output sData{4} and sData{5} are empty, but I need data in that array to run the next code to compute ROI. I tried with all the image data set but the array is always empty. I am doubtful since I use windows, is there a different condition for line 69 as it is mentioned 'Good enough for Linux, add conditions for MAC and Windows' and I do not know which conditions to use. If any one has previously done this please do help.

function [sData] = readDICOMdir(dicomPath,waitB)
% -------------------------------------------------------------------------
% function [sData] = readDICOMdir(dicomPath,waitB)
% -------------------------------------------------------------------------
% DESCRIPTION: 
% This function reads the DICOM content of a single directory. It then 
% organizes the data it in a cell of structures called 'sData', and 
% computes the region of interest (ROI) defined by a given RTstruct (if 
% present in the directory).
% -------------------------------------------------------------------------
% INPUTS:
% - dicomPath: Full path where the DICOM files to read are located.
% - waitB: Logical boolean. If true, a waiting bar will be displayed.
% -------------------------------------------------------------------------
% OUTPUTS:
% - sData: Cell of structures organizing the content of the volume data, 
%          DICOM headers, DICOM RTstruct* (used to compute the ROI) and 
%          DICOM REGstruct* (used to register a MRI volume to a PET volume)
%          * If present in the directory
%    --> sData{1}: Explanation of cell content
%    --> sData{2}: Imaging data and ROI defintion (if applicable)
%    --> sData{3}: DICOM headers of imaging data
%    --> sData{4}: DICOM RTstruct (if applicable)
%    --> sData{5}: DICOM REGstruct (if applicable)
% -------------------------------------------------------------------------
% AUTHOR(S): 
% - Martin Vallieres <[e-mail]>
% - Sebastien Laberge <[e-mail]>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>, 
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres, Sebastien Laberge
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------


% INITIALIZATION
if waitB
    waitbarHandle = waitbar(0,'Loading DICOM files...','WindowStyle','modal');
end
elements = dir(dicomPath);
nElements = length(elements);
volume = cell(1,1,nElements);
dicomHeaders = [];
RTstruct = [];
REG = [];


% READING DIRECTORY CONTENT
sliceNumber = 0;
for elementNumber = 1:nElements
    elementName = elements(elementNumber).name;
    if ~strcmp(elementName,'.') && ~strcmp(elementName,'..') % Good enough for Linux, add conditions for MAC and Windows.
        elementFullFile = fullfile(dicomPath,elementName);
        if isdicom(elementFullFile)
            tmp = dicominfo(elementFullFile);
            if strcmp(tmp.Modality,'RTSTRUCT')
                RTstruct = tmp;
            elseif strcmp(tmp.Modality,'REG')
                REG = tmp;
            elseif strcmp(tmp.Modality,'MR') || strcmp(tmp.Modality,'PT') || strcmp(tmp.Modality,'CT')
                sliceNumber = sliceNumber + 1;
                volume{sliceNumber} = double(dicomread(elementFullFile));
                dicomHeaders = appendStruct(dicomHeaders,tmp);
            end
        end
    end
    if waitB
        waitbar(elementNumber/nElements,waitbarHandle); 
    end
end
nSlices = sliceNumber; % Total number of slices
volume = volume(1:nSlices); % Suppress empty cells in images


% DETERMINE THE SCAN ORIENTATION
dist = [abs(dicomHeaders(2).ImagePositionPatient(1) - dicomHeaders(1).ImagePositionPatient(1)), ...
        abs(dicomHeaders(2).ImagePositionPatient(2) - dicomHeaders(1).ImagePositionPatient(2)), ...
        abs(dicomHeaders(2).ImagePositionPatient(3) - dicomHeaders(1).ImagePositionPatient(3))];
[~,index] = max(dist);
if index == 1
    orientation = 'Sagittal';
elseif index == 2
    orientation = 'Coronal';
else
    orientation = 'Axial';
end


% SORT THE IMAGES AND DICOM HEADERS
slicePositions = zeros(1,nSlices);
for sliceNumber = 1:nSlices
    slicePositions(sliceNumber) = dicomHeaders(sliceNumber).ImagePositionPatient(index);
end
[~,indices] = sort(slicePositions);
volume = cell2mat(volume(indices));
dicomHeaders = dicomHeaders(indices);


% FILL sData
sData = cell(1,5);
type = dicomHeaders(1).Modality;
if strcmp(type,'PT') || strcmp(type,'CT')
    if strcmp(type,'PT')
        type = 'PET';
    end
    for i=1:size(volume,3)
        volume(:,:,i)=volume(:,:,i)*dicomHeaders(i).RescaleSlope + dicomHeaders(i).RescaleIntercept;
    end
end
type = [type,'scan'];
sData{1} = struct('Cell_1','Explanation of cell content', ...
                  'Cell_2','Imaging data and ROI defintion (if applicable)', ...
                  'Cell_3','DICOM headers of imaging data', ...
                  'Cell_4','DICOM RTstruct (if applicable)', ...
                  'Cell_5','DICOM REGstruct (if applicable)');

sData{2}.scan.volume = volume;
sData{2}.scan.orientation = orientation;
try sData{2}.scan.pixelW = dicomHeaders(1).PixelSpacing(1); catch, sData{2}.scan.pixelW = []; end % Pixel Width
try sData{2}.scan.sliceT = dicomHeaders(1).SliceThickness; catch, sData{2}.scan.sliceT = []; end % Slice Thickness
s1 = round(0.5*nSlices); s2 = round(0.5*nSlices) + 1; % Slices selected to calculate slice spacing
sData{2}.scan.sliceS = sqrt(sum((dicomHeaders(s1).ImagePositionPatient - dicomHeaders(s2).ImagePositionPatient).^2)); % Slice Spacing
sData{2}.type = type;
sData{3} = dicomHeaders;
sData{4} = RTstruct;
sData{5} = REG;


% COMPUTE TUMOR DELINEATION USING RTstruct
if ~isempty(sData{4})
    [sData] = computeROI(sData);
end


if waitB
    close(waitbarHandle)
end

end


% UTILITY FUNCTION
function [structureArray] = appendStruct(structureArray,newStructure)

if isempty(structureArray)
    structureArray = newStructure;
    return
end

structLength = length(structureArray);
fields = fieldnames(structureArray(1));
nFields = length(fields);

for i = 1:nFields
    try
        structureArray(structLength + 1).(fields{i}) = newStructure.(fields{i});
    catch
        structureArray(structLength + 1).(fields{i}) = 'FIELD NOT PRESENT';
    end
end

end

Many thanks

[Walter Roberson]

The line you mentioned is fine for Mac and Windows as well.

Answer 1

[Mario Malic]

Hi Philip,

The code should work, does it?

Edited as per Walter's note

When you use dir specifying folder as an input argument, entries with name '.', '..' will appear somewhere in your dir structure.

This if statement is used just to avoid those two rows, because isdicom will produce an error when you supply to it the path that doesn't lead to an actual file.

 if ~strcmp(elementName,'.') && ~strcmp(elementName,'..')

[M Philip]

Hi,

The code is working, but in the output sData{4} is empty. I want content in sData{4} to run the next script 'computeROI'.

Thanks for the explanation for the code,

if ~strcmp(elementName,'.') && ~strcmp(elementName,'..')

Please help with it. I was stuck with this for many weeks and I am not a proficient user.

[Walter Roberson]

It is not correct to say that the first two entries will be those.

First of all, MS Windows does not define the order of entries returned by the internals used by dir(), saying that the order is up to the filesystem.

Secondly as far as I have been able to research, the NTFS file system does not define the order of entries returned, and I see reports that the actual order has changed twice over the history of NTFS. Because Microsoft does not define the order, it does not have to file any official announcement that the order is changing, since any program that relies on a particular order is using undefined behavior at its own risk.

Thirdly, the actual order depends upon the code points representing the names. I have not determined whether the order is by utf8 bytes or by UTF16 bytes or according to Unicode code points, or if it depends on the current user's locale. If I recall correctly, it depends upon the locale of the user who created the NTFS file system.

What this is leading up to is that there are valid printable filenames in NTFS which are known to sort before period. period is not the first printable character in any standard order. Not all of the printable characters before period are permitted in filenames but there are several that are and files starting with those will be before period and double period in practice (and you should not rely on the order anyhow.)

[Walter Roberson]

at the moment, as we do not know which file you have used, we have no reason to expect that those kinds of data are present in your dicom file.

[M Philip]

Hello @Walter Roberson,

This is the link to the website from where I got the PET/CT images.

https://wiki.cancerimagingarchive.net/display/Public/Head-Neck-PET-CT

This was the input argument:

readDICOMdir('E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290','True')

In the command window output appears as

1×5 cell array

{1×1 struct} {1×1 struct} {1×209 struct} {0×0 double} {0×0 double}

These are the screenshot of the output I got.

Thanks.

[Mario Malic]

The file is a bit too big to download.

This part of the code identifies what file contains, and it looks like, there might not be 'RTSTRUCT' file in here. Is there supposed to be one?

if strcmp(tmp.Modality,'RTSTRUCT')
    RTstruct = tmp;
elseif strcmp(tmp.Modality,'REG')
    REG = tmp;
elseif strcmp(tmp.Modality,'MR') || strcmp(tmp.Modality,'PT') || strcmp(tmp.Modality,'CT')
    sliceNumber = sliceNumber + 1;
    volume{sliceNumber} = double(dicomread(elementFullFile));
    dicomHeaders = appendStruct(dicomHeaders,tmp);
end

Paste the second line in your code like this and search in the command window if you can find RTSTRUCT. Maybe these files have the different name like RT_STRUCT or similar.

tmp = dicominfo(elementFullFile);
tmp.Modality
if strcmp(tmp.Modality,'RTSTRUCT')

[M Philip]

Hello @Mario Malic,

The file is 72GB in size, a bit too big.

As you said, I can find a folder named 'RT Struct' outside the folder which contains these images. But I am not sure how I will set both these folders as the dicom path.

When I include both those folders as input argument, it shows 'too many arguments'.

If I just give the RTStruct folder, it shows error in the line 'isempty'

Is there a way to include both these folders as input?

[Mario Malic]

If you get an error which explains itself, first thing that you should take a look is at function documentation and see how many input arguments does it take (answer is one).

doc dir

Let's take a look at directories first, what is the folder that contains the whole data relevant for that case, including the folder with RTSTRUCT data?

E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290

% Folders
| HN-HMR-039
    | 08-27-1885-TEP PANC. SPHERE ORL-70517
        |08-27-1885-NASOPHARYNX IMRT C-75645
            | 2.000000-O-MAR-51290
                |% now we're within 2.000000-O-MAR-51290
                |dicomfile1.dicom
                |dicomfile2.dicom
                |RTSTRUCT % folder

If the files for the relevant case including RTSTRUCT ones are within "2.000000-O-MAR-51290" folder, you can use this path to call the readDICOMdir

dicomPath = "readDICOMdir('E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290\**','True')"

This is almost like the regular path, but with wildcard '*'. When you add a single asterisk, it will look for files and folders within specified folder. The double asterisk will look within specified folder and in all the folders within it.

---

Regarding the mail: it can be hard for everyone to understand the code very in short time, especially if it comes from another person. It's always best to look for information in the documentation (I haven't read the documentation for a function, which I'll do now), here's the relevant part.

% DESCRIPTION:

% This function reads the DICOM content of a single directory. It then

% organizes the data it in a cell of structures called 'sData', and

% computes the region of interest (ROI) defined by a given RTstruct (if

% present in the directory).

% -------------------------------------------------------------------------

% INPUTS:

% - dicomPath: Full path where the DICOM files to read are located.

% - waitB: Logical boolean. If true, a waiting bar will be displayed.

I don't know anything about DICOM files, or how are they placed within the folders, so, I hope it works, if it doesn't, send a screenshot or maybe briefly describe which folder has the relevant files.

Note: The issue with the double asterisk is, if the folder contains thousands of files, but only 20 are DICOM files to be analysed, it'll be a huge time loss.

[M Philip]

dicomPath = "readDICOMdir('E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290\**','True')"

I tried the above path, but it shows error as below.

Also the 'RTstruct' folder is outside '2.000000-O-MAR-51290' i.e, only '2.000000-O-MAR-51290' is within'08-27-1885-NASOPHARYNX IMRT C-75645'. But 'RTStruct' is outside '08-27-1885-NASOPHARYNX IMRT C-75645' and within '08-27-1885-TEP PANC. SPHERE ORL-70517'

The RTSTRUCT folder contain only one dicom image.

The screenshot of path is given below.

[Mario Malic]

I tested out yesterday isdicom, it throws an error if the path to folder is given to it, that's most likely the error there. I forgot that those are included when you use dir with wildcard.

Alright, then copy the file from RTSTRUCT folder into the folder 2.000000-O-MAR-51290 and call the readDICOMdir with folder your from earlier comment.

readDICOMdir('E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290','True')

[Mario Malic]

No no, copy only the file 1-1.dcm into the folder 2.000000-O-MAR-51290

[M Philip]

Hello,

Thanks for that.

I copied RT struct file into the folder ' 2.000000-O-MAR-51290 '.

But I am getting the following error.

Warning: Attribute (3006,0050) at byte 20918450 might have been incorrectly interpreted as explicit VR. If reading this
file fails, set the 'UseVRHeuristic' parameter to false. 
> In dicominfo (line 63)
In readDICOMdir (line 72)
 
Warning: Not enough data imported.  Attempted to read 761017654 bytes at position 20918466.  Only read 51232612. 
> In dicominfo (line 63)
In readDICOMdir (line 72)
 
Warning: File and dictionary VR values do not match for attribute "ContourData" (3006,0050). Expected: DS. Found: FL.
Consider setting the 'UseDictionaryVR' parameter to true. 
> In dicominfo>convertRawAttr (line 152)
In dicominfo>processMetadata (line 114)
In dicominfo>parseSequence (line 488)
In dicominfo>convertRawAttr (line 209)
In dicominfo>processMetadata (line 114)
In dicominfo>parseSequence (line 488)
In dicominfo>convertRawAttr (line 209)
In dicominfo>processMetadata (line 114)
In dicominfo (line 72)
In readDICOMdir (line 72) 
Unrecognized function or variable 'sliceOK'.

Error in computeROI (line 85)
            pts(:,1) = pts(:,1) - sData{3}(sliceOK).ImagePositionPatient(1);

Error in readDICOMdir (line 148)
    [sData] = computeROI(sData);

[Mario Malic]

There's some progress, it looks like MATLAB is throwing errors on the RTSTRUCT '1-1.dcm' file. You can try to test Name-Value parameters suggested in the warning messages.

% rtsFile is path to rtstruct file
rtsFile = 'E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290\1-1.dcm';

One of the next three lines will probably throw no warnings and the file will be succesfully read. Once you find out, which one is it, you'll have to modify the readDICOMdir function to include it.

tmp = dicominfo(rtsFile, 'UseVRHeuristic', false);
tmp = dicominfo(rtsFile, 'UseDictionaryVR', true);
tmp = dicominfo(rtsFile, 'UseVRHeuristic', false, 'UseDictionaryVR', true);

% readDICOMdir snippet
if isdicom(elementFullFile)
    [~, fileName, ~] = fileparts(elementFullFile);
    if strcmpi(fileName, '1-1') % This modifies how you the RTSTRUCT file is read
        tmp = dicominfo(elementFullFile, 'Name', Value); % The line that doesn't cause the warnings/errors
    else
        tmp = dicominfo(elementFullFile);
    end
    if strcmp(tmp.Modality,'RTSTRUCT')
        RTstruct = tmp;
    elseif strcmp(tmp.Modality,'REG')
        REG = tmp;
    elseif strcmp(tmp.Modality,'MR') || strcmp(tmp.Modality,'PT') || strcmp(tmp.Modality,'CT')
        sliceNumber = sliceNumber + 1;
        volume{sliceNumber} = double(dicomread(elementFullFile));
        dicomHeaders = appendStruct(dicomHeaders,tmp);
    end
end

[M Philip]

Yes, I am happy that there is some progress. In the meanwhile I copied the file within REGCTSIM to '2.000000-O-MAR-51290' but not RTSTRUCT and that gave me value in the sData{5} array which was previously empty.

I just want to know ,Where to type the rtsFILE path?While declaring the function have to include rtsFILE also along with dicompatha and waitB? and also the code starting tmp=?(These questions may sound silly, but I am not experinced at using matlab)

The code snippet you gave is after making the above changes or the change includes that too?

[Mario Malic]

This is for you to try in the command prompt, it's not a part of the function. You need to test which one of these dicominfo lines work without warnings/errors.

% rtsFile is path to rtstruct file
elementFullFile = 'E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290\1-1.dcm';
tmp = dicominfo(elementFullFile, 'UseVRHeuristic', false);
tmp = dicominfo(elementFullFile, 'UseDictionaryVR', true);
tmp = dicominfo(elementFullFile, 'UseVRHeuristic', false, 'UseDictionaryVR', true);

If the third line on the above snippet works well, this is how the part of code in readDICOMdir should look like. If one of the other lines works well, then replace it

if isdicom(elementFullFile)
    [~, fileName, ~] = fileparts(elementFullFile);
    if strcmpi(fileName, '1-1') % This modifies only how you read the RTSTRUCT file
        tmp = dicominfo(elementFullFile, 'UseVRHeuristic', false, 'UseDictionaryVR', true); % -To replace with the line that doesn't cause the warnings/errors
    else
        tmp = dicominfo(elementFullFile);
    end
    if strcmp(tmp.Modality,'RTSTRUCT')
        RTstruct = tmp;
    elseif strcmp(tmp.Modality,'REG')
        REG = tmp;
    elseif strcmp(tmp.Modality,'MR') || strcmp(tmp.Modality,'PT') || strcmp(tmp.Modality,'CT')
        sliceNumber = sliceNumber + 1;
        volume{sliceNumber} = double(dicomread(elementFullFile));
        dicomHeaders = appendStruct(dicomHeaders,tmp);
    end
end

[M Philip]

Thank you.

Line one and three of the below code worked well.

% rtsFile is path to rtstruct file
elementFullFile = 'E:\PhD\h&n data\Head-Neck-PET-CT\HN-HMR-039\08-27-1885-TEP PANC. SPHERE ORL-70517\08-27-1885-NASOPHARYNX IMRT C-75645\2.000000-O-MAR-51290\1-1.dcm';
tmp = dicominfo(elementFullFile, 'UseVRHeuristic', false);
tmp = dicominfo(elementFullFile, 'UseDictionaryVR', true);
tmp = dicominfo(elementFullFile, 'UseVRHeuristic', false, 'UseDictionaryVR', true);

But when I made the corresponding changes in readDICOMdir.m it shows error as below.

Unrecognized function or variable 'sliceOK'.

Error in computeROI (line 85)
            pts(:,1) = pts(:,1) - sData{3}(sliceOK).ImagePositionPatient(1);

Error in readDICOMdir (line 153)
    [sData] = computeROI(sData)

[Mario Malic]

Error is in ComputeROI.m. Since UIDRT wasn't equal to UIDSlice, sliceOK doesn't exist.

Well, this will require some digging into contents of your sData to find out why. Maybe different name for itemSlice or missing data?

    for sliceNumber = 1:nSlices
        
        % Find slice correspondence between volume and RTstruct
        itemSlice = ['Item_',num2str(sliceNumber)];
        UIDrt = sData{4}.ROIContourSequence.(itemContour).ContourSequence.(itemSlice).ContourImageSequence.Item_1.ReferencedSOPInstanceUID;
        for i = 1:sizeScan(3)
            UIDslice = sData{3}(i).SOPInstanceUID;
            if strcmp(UIDrt,UIDslice)
                sliceOK = i; 
                break
            end
        end
        
        % code...
    end

[M Philip]

I will try this code and where to plug this in? (I haven't fully understood the code yet, what each line does)

[Mario Malic]

If you see the error that you posted, it says on which line error occured and in which file.

Unrecognized function or variable 'sliceOK'. 
Error in computeROI (line 85) % file and line
            pts(:,1) = pts(:,1) - sData{3}(sliceOK).ImagePositionPatient(1);
Error in readDICOMdir (line 153) 
    [sData] = computeROI(sData)

Now you have to take a look at computeROI.m and figure out why is this if statement never true

if strcmp(UIDrt,UIDslice)
    sliceOK = i; 
    break
end

[M Philip]

When I run the readDICOMdir.m script does it run the computeROI.m script also? There are three more scripts like getROIonly.m, getROIbox.m and computeSUVmap.m under the same folder 'Utilities' in github. There is no need to run all these separately?

How can I figure out why it is never true? is it from the output sData ?

[Mario Malic]

Of course it does, when you call readDICOMdir.m . Somewhere within the code computeROI is called

% COMPUTE TUMOR DELINEATION USING RTstruct
if ~isempty(sData{4})
    [sData] = computeROI(sData);
end

It means it sends sData into it and starts executing the lines in computeROI.m

• There is no need to run all these separately? I don't think so, but I can't be sure as I am not familiar with the code.
• How can I figure out why it is never true? is it from the output sData ? strcmp is the function that compares strings and if they are not equal, false is returned as result. I have just noticed that, there are two for loops before this snippet below, and I guess you'll have to start to search if the code in the loop works correctly at all, does contourNumber or sliceNumber takes any value other than 1? You might have to start debugging.

for contourNumber = 1:nContours
    % code
    for sliceNumber = 1:nSlices
        % code
        % Find slice correspondence between volume and RTstruct
        for i = 1:sizeScan(3)
            UIDslice = sData{3}(i).SOPInstanceUID;
            if strcmp(UIDrt,UIDslice) % strings are compared on this line
                sliceOK = i;
                break
            end
        end

I think that if you want to solve this issue, there are few things that you need to find out:

• structure of the RTSTRUCT file, what does it contain? These two variables should match, why do they don't? Does the data output from imaging devices differ between different manufacturers?

UIDrt = sData{4}.ROIContourSequence.(itemContour).ContourSequence.(itemSlice).ContourImageSequence.Item_1.ReferencedSOPInstanceUID;
UIDslice = sData{3}(i).SOPInstanceUID;

If there are plenty of tasks similar to this one are ahead of you, I think you have to brush up on your programming skills, it should make things easier in the long term.

[M Philip]

Hello @Mario Malic

Learning matlab is something new for me. I am learning a new skill now. If you have any suggestions on some websites from where I can learn this and practice, it would be helpful. I would be using more of this in the future too. How can I rectify find and rectify errors easily, does it com only with practice.?

Regardsing the RTSTRUCT, since I got it from a website with Public Access I am not sure what it contains. I have contacted the person behind the code and data but have not got any response. What else I can do to solve it.

I tried the same code with another path.

I got the error as

Undefined function 'computeBoundingBox' for input arguments of type 'double'.

Error in computeROI (line 106)
    [boxBound] = computeBoundingBox(mask);

Error in readDICOMdir (line 151)
    [sData] = computeROI(sData);

There are almost 300 folders and many subfolders within each. I am not able to find out what each folder contains by just seeing the name, that is another problem which I face. So I am not sure whether the RTStruct file is the only one which I have to copy to the image older or some other files have to be copied for the script to run.

[M Philip]

Also to debug the computeROI.m script, what will be the input argument? Is it sData? How to do that? I could do for readDICOM.dir but for computeROI.m I am not sure.

[Mario Malic]

Is MATLAB your first programming language? You should start with MATLAB Onramp. There's Cody on Mathworks Central which is a some sort of game/competition with different problem groups across all levels, from beginner to expert. It's great way to start with learning MATLAB.

How can you rectify errors quickly? What does this error tells you?

Undefined function 'computeBoundingBox' for input arguments of type 'double'.
[boxBound] = computeBoundingBox(mask);

Now I went over to this function to read the description and it says

% DESCRIPTION: 
% This function computes the smallest box containing the whole region of 
% interest (ROI). It is adapted from the function compute_boundingbox.m
% of CERR <http://www.cerr.info/>.
% -------------------------------------------------------------------------
% INPUTS:
% - mask: 3D array, with 1's inside the ROI, and 0's outside the ROI.
% -------------------------------------------------------------------------
% OUTPUTS:
% - boxBound: Bounds of the smallest box containing the ROI. 
%             Format: [minRow, maxRow;
%                      minColumn, maxColumns;
%                      minSlice, maxSlice]
% -------------------------------------------------------------------------

I check INPUTS and it mentions a 3d array with 1's and 0's. Maybe the variable mask already has 1s and 0s, but it's in the wrong variable type - double. It doesn't say what variable type should be input, if it's 1s and 0s, it probably accepts logical array. Why is this case not covered? Probably because the code was written on the older version of MATLAB and newer version has some changes that happened in their functions, so the code has to be adapted to cover it.

Then I take a look at mask, where is it defined, where its value is changed, lines 60, 99, 101? If you want try to fix without checking anything, you could simply just add following line before the computeBoundingBox line, but I wouldn't encourage you to do so.

if ~islogical(mask)
    mask = logical(mask);
end

If you take a look at line 99, it has function or(some code inside), and it returns a logical array, but you send a double(numeric) array into function, it means that line if condition from line 79 is false, which means that pts_temp is empty and line 101 is executed

if ~isempty(pts_temp) % line 79

If you take a look at it, it calls zeros which creates a double(array), and now you see the cause of your error (maybe).

mask(:,:,sliceOK) = zeros(sizeScan(1),sizeScan(2)); % line 101

How to fix it? Actually, question is can you or should you?

Line below could work, if the computeBoundingBox function actually accepts logical array.

mask(:,:,sliceOK) = false(sizeScan(1),sizeScan(2));

Let's stop a bit with the code, and not to try to fix the errors when they are thrown. There might be so many of them until you get some result.

But still, it is essential to understand the what's going on in the code first, and then the what each line in the code does.

In the git repository readme there are other files that you can download and see how they are structured in the folder, but I doubt it's the reason. Maybe the author wrote code for one standard of dicom files, whereas in your files standard is different.

Maybe you could ask the author about his opinion on your situation and maybe he can shed a suggestion or two about code, files, books, or anything that can help you.

----

If the variable type doesn't sound familiar to you, you should find an introductory course on programming in MATLAB, or maybe even programming in general. I had two college introductory courses in C, then MATLAB afterwards, but I have learned a lot on MATLAB Answers and Cody by asking and looking for answers. There's also a question that is worthy to read - link.

[M Philip]

Matlab is not my first programming language. I have learned Python as part of my curriculum but never used otherwise.

The author of the code has used the same images from the Cancer Imaging Archive for their study. I am also using the same set of images and same code written by them.

I tried using the line

mask(:,:,sliceOK) = false(sizeScan(1),sizeScan(2));

It shows error as:

Unrecognized function or variable 'computeBoundingBox'.

Error in computeROI (line 106)
    [boxBound] = computeBoundingBox(mask);

Error in readDICOMdir (line 151)
    [sData] = computeROI(sData);

Thank you for sharing the last link and for pointing out to some resources. I will try to use those to better understand Matlab. I want to understand each line of the code and how it is executed and how to correct errors.

Thanks for being patient with me and my questions.

[Mario Malic]

That error is a bit weird now, because it worked before, now MATLAB can't find the function computeBoundingBox. This function is in another folder so make sure MATLAB can find it by adding it to path.

https://www.mathworks.com/help/matlab/ref/pathtool.html, click the button 'add with subfolders' and select the folder that contains that repository.

It can be overwhelming with the code at the moment, but I think that someone with decent knowledge of MATLAB would be a bit stuck with this as well, because it requires knowledge about medical imaging, related software and files in addition to MATLAB.

You're welcome, unfortunately it didn't end up with a solution, but direction is not so bad either, just stay with it.

Good luck!

Edit:

Apparently, this function reads RTSTRUCT file, take a look. There are more resources if you Google.

https://www.mathworks.com/help/images/ref/dicomcontours.html

https://github.com/ulrikls/dicomrt2matlab

[M Philip]

Hey Mario,

It worked now!!!! Thanks very much.You were very helpful. I wonder how you figure out the errors so easily and fix it!

Now I have to understand what each line of the code means and whether the output I get is what I want.

I added the other folder to the path as you said. The first link is broken. But I found it from Matlab documentation. Thank you for explaining each step of it.

This is how the output looks like.

Can I ask similar questions in the future if I am stuck?

Thanks again.

[Mario Malic]

Yay! I am glad it worked!

This code worked? Or, did you use one of those linked functions?

mask(:,:,sliceOK) = false(sizeScan(1),sizeScan(2));

Well, I do not know if I will be able to help, but here's the thing with this forum. If you can describe your question and problem precisely, provide files needed to run the code, it's very very likely that someone will help you. There are few people here who know a lot about medical imaging, including DICOM files.

Asking questions the right way

[M Philip]

I used the code which you gave.

mask(:,:,sliceOK) = false(sizeScan(1),sizeScan(2));

Thanks for sharing the useful links. I will learn how to ask questions which are very precise.

I really appreciate your help!

Thanks again.