Radiomics - github code https://github.com/mvallieres/radiomics

Question

M Philip 2021-2-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/751909-radiomics-github-code-https-github-com-mvallieres-radiomics

评论： M Philip 2021-2-23

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 <mart.vallieres@gmail.com>
% - Sebastien Laberge <sebastien.laberge.3000@gmail.com>
% -------------------------------------------------------------------------
% 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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Walter Roberson 2021-2-21

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mario Malic 2021-2-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/751909-radiomics-github-code-https-github-com-mvallieres-radiomics#answer_629364

编辑：Mario Malic 2021-2-21

在 MATLAB Online 中打开

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,'..')

28 个评论
显示 26更早的评论隐藏 26更早的评论

Walter Roberson 2021-2-21

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.)

Mario Malic 2021-2-22

在 MATLAB Online 中打开

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 2021-2-22

在 MATLAB Online 中打开

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 2021-2-22

在 MATLAB Online 中打开

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

Mario Malic 2021-2-22

编辑：Mario Malic 2021-2-22

在 MATLAB Online 中打开

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

Mario Malic 2021-2-22

在 MATLAB Online 中打开

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 2021-2-23

在 MATLAB Online 中打开

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.

Mario Malic 2021-2-23

在 MATLAB Online 中打开

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.

Mario Malic 2021-2-23

编辑：Mario Malic 2021-2-23

在 MATLAB Online 中打开

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 2021-2-23

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Radiomics - github code https://github.com/mvallieres/radiomics

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

28 个评论
显示 26更早的评论隐藏 26更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

28 个评论 显示 26更早的评论隐藏 26更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Radiomics - github code https://github.com/mvallieres/radiomics

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

28 个评论
显示 26更早的评论隐藏 26更早的评论