Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (feat_muscle_subdivisions) Browse

R8218/dicominfoSCase.m
7434b144f712feat_muscle_subdivisions

Last Change
Raw File

dicominfoSCase.m
View Options

function outputArg = dicominfoSCase(varargin)
%DICOMINFOSCASE Get info from dicom fields and save in csv, xls, mat, and SQL
% This function should be started from /shoulder/methods/matlab/database
% A logfile is written during execution in
% /shoulder/methods/matlab/database/log/dicominfoSCase.log
% If the function is called without arguments the csv file is overwritten,
% otherwise it is updated. The same applies for xls, mat and SQL.
% The fonction can be run from server (lbovenus) with
% cd /home/shoulder/methods/matlab/database; matlab -nojvm -nodisplay -nosplash -r "dicominfoSCase;quit"
% For 700 cases, it takes less than 60 seconds (when executed from server)
% Progress can be be checked throuhth log file with following shell command
% cd /home/shoulder/methods/matlab/database;tail -f log/dicominfoSCase.log
% This funtion replaces
% /home/shoulder/methods/matlab/database/CTinfoCT.m
% which only works matlab R2015a on a Windows system
% Syntax: outputArg = dicominfoSCase(varargin)
%
% Input:
% varargin: can be nothing (all cases), 'N', 'P', or SCase(s)
%
% Output:
% outputArg: 1 if ok and 0 otherwise
%
% Example: CTdicomInfo, CTdicominfo('P'), CTdicominfo('P315')
% Other M-files required:
% Subfunctions: listSCase
% See also: measureSCase
%
% Author: Alexandre Terrier
% EPFL-LBO
% Creation date: 2018-07-01
% Revision date: 2018-08-29
%
% TO DO:
% Re-check for double slash in paths
% We should implement the ShoulderCase class
% We should save data in SCase/matlab directory
% We might add family name, given name (to check anonymity)
% We should first read csv and then update
% We should update MySQL
% Check 'Index exceeds matrix dimensions' when 'using dicomdisp'
% Structure SCaseDicom should match csv (and xlx aand sql) firld names
% rather than dicom
%% Initialization
% This should be in a function call
dataDir = '../../../data/'; % location of data
xlsDataDir = strcat(dataDir, 'Excel/xlsFromMatlab/'); % Location of the XLS ShoulderDatabase
logDir = 'log'; % log file in current directory
%% Open log file
filename = 'dicominfoSCase.log';
filename = [logDir '/' filename];
logFileId = fopen(filename, 'w');
fprintf(logFileId, 'dicominfoSCase log file\n');
fprintf(logFileId, '\nDate: ');
fprintf(logFileId, datestr(datetime));
fprintf(logFileId, '\n');
fprintf(logFileId, '\nGet dicom info of SCase');
fprintf(logFileId, '\n');
%% Get list of SCases
if nargin == 0 % dicominfoSCase called without argument
SCaseList = listSCase; %listSCase is an external function which provides a list with all the SCases in Excel DB (commented by JSM)
update = 0;
elseif nargin == 1
update = 1;
inputArg = varargin{1};
SCaseList = listSCase(inputArg);
end
%% Read dicom info
tic; % Start stopwatch timer
fprintf(logFileId, '\n\nRead dicom meta data\n'); % log file message
% Struct below correponds to regular dicom field names, except SCase_ID
SCaseDicom = struct(...
'SCase_id' ,[],...
'PatientID' ,[],...
'PatientSex' ,[],...
'PatientAge' ,[],...
'PatientBirthDate' ,[],...
'PatientSize' ,[],...
'PatientWeight' ,[],...
'PatientFamilyName' ,[],...
'PatientGivenName' ,[],...
'AcquisitionDate' ,[],...
'ConvolutionKernel' ,[],...
'PixelSpacing' ,[],...
'SpacingBetweenSlices' ,[],...
'SliceThickness' ,[],...
'KVP' ,[],...
'XrayTubeCurrent' ,[],...
'InstitutionName' ,[],...
'StudyDescription' ,[],...
'SeriesDescription' ,[],... % Might describe shoulder side
'Manufacturer' ,[],...
'ManufacturerModelName' ,[],...
'ProtocolName' ,[]...
);
% Loop over all SCases collected above
for SCase_idx=1:1:numel(SCaseList)
SCase_id = SCaseList(SCase_idx).id;
dicomDir = [SCaseList(SCase_idx).dir '/dicom'];
dicomFiles = dir([dicomDir, '/*.dcm']);
dicomFile1 = [dicomDir, '/', dicomFiles(10).name]; % seems better to check 10th than 1st !
dicomInfoError = false;
fprintf(logFileId, '\n');
fprintf(logFileId, SCaseList(SCase_idx).id);
fprintf(logFileId, ': ');
% Try first with matlab function dicominfo
try
dicomInfo1 = dicominfo(dicomFile1, 'UseVRHeuristic', true);
fprintf(logFileId, 'dicominfo ok');
% If dicominfo didn't worked, try with matlab function dicomdisp
catch ME
dicomInfoError = true;
fprintf(logFileId, '\n');
fprintf(logFileId, 'dicominfo error --> using dicomdisp');
fprintf(logFileId, ['\n' ME.message]);
dicomdispStr = evalc('dicomdisp(dicomFile1)');
end
% If dicominfo didn't worked, continue with output of dicomdisp
if dicomInfoError % get dicom info from dicomdisp()
PatientID = extractAfter(dicomdispStr,'PatientID');
PatientID = extractAfter(PatientID,'[');
PatientID = extractBefore(PatientID,']');
PatientSex = extractAfter(dicomdispStr,'PatientSex');
PatientSex = extractAfter(PatientSex,'[');
PatientSex = extractBefore(PatientSex,']');
PatientAge = extractAfter(dicomdispStr,'PatientAge');
PatientAge = extractAfter(PatientAge,'[');
PatientAge = extractBefore(PatientAge,']');
PatientBirthDate = extractAfter(dicomdispStr,'PatientBirthDate');
PatientBirthDate = extractAfter(PatientBirthDate,'[');
PatientBirthDate = extractBefore(PatientBirthDate,']');
PatientSize = extractAfter(dicomdispStr,'PatientSize');
PatientSize = extractAfter(PatientSize,'[');
PatientSize = extractBefore(PatientSize,']');
PatientSize = str2double(PatientSize);
PatientWeight = extractAfter(dicomdispStr,'PatientWeight');
PatientWeight = extractAfter(PatientWeight,'[');
PatientWeight = extractBefore(PatientWeight,']');
PatientWeight = str2double(PatientWeight);
AcquisitionDate = extractAfter(dicomdispStr,'AcquisitionDate');
AcquisitionDate = extractAfter(AcquisitionDate,'[');
AcquisitionDate = extractBefore(AcquisitionDate,']');
ConvolutionKernel = extractAfter(dicomdispStr,'ConvolutionKernel');
ConvolutionKernel = extractAfter(ConvolutionKernel,'[');
ConvolutionKernel = extractBefore(ConvolutionKernel,']');
PixelSpacing = extractAfter(dicomdispStr,'PixelSpacing');
PixelSpacing = extractAfter(PixelSpacing,'[');
PixelSpacing = extractBefore(PixelSpacing,'\');
SpacingBetweenSlices = extractAfter(dicomdispStr,'SpacingBetweenSlices');
SpacingBetweenSlices = extractAfter(SpacingBetweenSlices,'[');
SpacingBetweenSlices = extractBefore(SpacingBetweenSlices,']');
SpacingBetweenSlices = str2double(SpacingBetweenSlices);
SliceThickness = extractAfter(dicomdispStr,'SliceThickness');
SliceThickness = extractAfter(SliceThickness,'[');
SliceThickness = extractBefore(SliceThickness,']');
SliceThickness = str2double(SliceThickness);
KVP = extractAfter(dicomdispStr,'KVP');
KVP = extractAfter(KVP,'[');
KVP = extractBefore(KVP,']');
KVP = str2double(KVP);
XrayTubeCurrent = extractAfter(dicomdispStr,'XrayTubeCurrent');
XrayTubeCurrent = extractAfter(XrayTubeCurrent,'[');
XrayTubeCurrent = extractBefore(XrayTubeCurrent,']');
XrayTubeCurrent = str2double(XrayTubeCurrent);
InstitutionName = extractAfter(dicomdispStr,'InstitutionName');
InstitutionName = extractAfter(InstitutionName,'[');
InstitutionName = extractBefore(InstitutionName,']');
Manufacturer = extractAfter(dicomdispStr,'Manufacturer');
Manufacturer = extractAfter(Manufacturer,'[');
Manufacturer = extractBefore(Manufacturer,']');
ManufacturerModelName = extractAfter(dicomdispStr,'ManufacturerModelName');
ManufacturerModelName = extractAfter(ManufacturerModelName,'[');
ManufacturerModelName = extractBefore(ManufacturerModelName,']');
ProtocolName = extractAfter(dicomdispStr,'ProtocolName');
ProtocolName = extractAfter(ProtocolName,'[');
ProtocolName = extractBefore(ProtocolName,']');
SeriesDescription = extractAfter(dicomdispStr,'SeriesDescription');
SeriesDescription = extractAfter(SeriesDescription,'[');
SeriesDescription = extractBefore(SeriesDescription,']');
else % dicomInfo exists
if isfield(dicomInfo1,'PatientID')
PatientID = dicomInfo1.PatientID;
else
PatientID = '';
end
if isfield(dicomInfo1,'PatientSex')
PatientSex = dicomInfo1.PatientSex;
else
PatientSex = '';
end
if isfield(dicomInfo1,'PatientBirthDate')
PatientBirthDate = dicomInfo1.PatientBirthDate;
else
PatientBirthDate = '';
end
if isfield(dicomInfo1,'PatientAge')
PatientAge = dicomInfo1.PatientAge;
else
PatientAge = '';
end
if isfield(dicomInfo1,'PatientSize')
PatientSize = dicomInfo1.PatientSize;
else
PatientSize = '';
end
if isfield(dicomInfo1,'PatientWeight')
PatientWeight = dicomInfo1.PatientWeight;
else
PatientWeight='';
end
if isfield(dicomInfo1,'AcquisitionDate')
AcquisitionDate = dicomInfo1.AcquisitionDate;
else
AcquisitionDate = '';
end
if isfield(dicomInfo1,'ConvolutionKernel')
ConvolutionKernel = dicomInfo1.ConvolutionKernel;
else
ConvolutionKernel = '';
end
if isfield(dicomInfo1,'PixelSpacing')
PixelSpacing = dicomInfo1.PixelSpacing(1);
PixelSpacing = num2str(PixelSpacing);
% PixelSpacing = [,PixelSpacing,]; % Commented because I (AT) don't
% understand the effect
else
PixelSpacing='';
end
if isfield(dicomInfo1,'SpacingBetweenSlices')
SpacingBetweenSlices = dicomInfo1.SpacingBetweenSlices;
else % SpacingBetweenSlices estimated from 2 slices
% fprintf(logFileID, '\n');
fprintf(logFileId, ', SpacingBetweenSlices estimated from 2 slices');
dicomFile2 = [dicomDir, '/', dicomFiles(11).name];
dicomInfo2 = dicominfo(dicomFile2, 'UseVRHeuristic', true);
ImagePositionPatient1 = dicomInfo1.ImagePositionPatient(3);
ImagePositionPatient2 = dicomInfo2.ImagePositionPatient(3);
posdDiff = abs(ImagePositionPatient2 - ImagePositionPatient1);
InstanceNumber1 = dicomInfo1.InstanceNumber;
InstanceNumber2 = dicomInfo2.InstanceNumber;
InstanceNumberDiff = abs(InstanceNumber2 - InstanceNumber1);
if InstanceNumberDiff == 0
InstanceNumberDiff=1;
end
SpacingBetweenSlices = posdDiff/InstanceNumberDiff;
end
if isfield(dicomInfo1,'SliceThickness')
SliceThickness = dicomInfo1.SliceThickness;
else
SliceThickness='';
end
if isfield(dicomInfo1,'KVP')
KVP = dicomInfo1.KVP;
else
KVP= '';
end
if isfield(dicomInfo1,'XrayTubeCurrent')
XrayTubeCurrent = dicomInfo1.XrayTubeCurrent;
else
XrayTubeCurrent = '';
end
if isfield(dicomInfo1,'InstitutionName')
InstitutionName = dicomInfo1.InstitutionName;
else
InstitutionName = '';
end
if isfield(dicomInfo1,'Manufacturer')
Manufacturer = dicomInfo1.Manufacturer;
else
Manufacturer ='';
end
if isfield(dicomInfo1,'ManufacturerModelName')
ManufacturerModelName = dicomInfo1.ManufacturerModelName;
else
ManufacturerModelName = '';
end
if isfield(dicomInfo1,'ProtocolName')
ProtocolName = dicomInfo1.ProtocolName;
else
ProtocolName = '';
end
if isfield(dicomInfo1,'SeriesDescription')
SeriesDescription = dicomInfo1.SeriesDescription;
else
SeriesDescription = '';
end
% End of reading dimcom info
end
%% Format data
PatientBirthDate = datetime(PatientBirthDate,'InputFormat','yyyyMMdd');
PatientBirthDate = datestr(PatientBirthDate, 'yyyy-mm-dd');
PatientAge = str2double(strtok(PatientAge,'Y'));
PatientSize = round(100 * PatientSize);
if PatientSize == 0
PatientSize = '';
end
PatientWeight = round(PatientWeight);
if PatientWeight == 0
PatientWeight = '';
end
AcquisitionDate = datetime(AcquisitionDate,'InputFormat','yyyyMMdd');
AcquisitionDate = datestr(AcquisitionDate, 'yyyy-mm-dd');
PixelSpacing = str2double(PixelSpacing);
PixelSpacing = round(PixelSpacing, 3);
SpacingBetweenSlices = round(SpacingBetweenSlices, 3);
SliceThickness = round(SliceThickness, 3);
KVP = round(KVP);
XrayTubeCurrent = round(XrayTubeCurrent);
% Replace any ',' by ' ' to avoid problem in csv
InstitutionName = strrep(InstitutionName,',',' ');
Manufacturer = strrep(Manufacturer,',',' ');
ManufacturerModelName = strrep(ManufacturerModelName,',',' ');
ProtocolName = strrep(ProtocolName,',',' ');
SeriesDescription = strrep(SeriesDescription,',',' ');
%% Add data in SCaseDicom structure
% This could be directly in SCase.CT object
SCaseDicom(SCase_idx).id = SCase_id;
SCaseDicom(SCase_idx).PatientID = PatientID;
SCaseDicom(SCase_idx).PatientSex = PatientSex;
SCaseDicom(SCase_idx).PatientBirthDate = PatientBirthDate;
SCaseDicom(SCase_idx).PatientAge = PatientAge;
SCaseDicom(SCase_idx).PatientSize = PatientSize;
SCaseDicom(SCase_idx).PatientWeight = PatientWeight;
SCaseDicom(SCase_idx).AcquisitionDate = AcquisitionDate;
SCaseDicom(SCase_idx).ConvolutionKernel = ConvolutionKernel;
SCaseDicom(SCase_idx).PixelSpacing = PixelSpacing;
SCaseDicom(SCase_idx).SpacingBetweenSlices = SpacingBetweenSlices;
SCaseDicom(SCase_idx).SliceThickness = SliceThickness;
SCaseDicom(SCase_idx).KVP = KVP;
SCaseDicom(SCase_idx).XrayTubeCurrent = XrayTubeCurrent;
SCaseDicom(SCase_idx).InstitutionName = InstitutionName;
SCaseDicom(SCase_idx).Manufacturer = Manufacturer;
SCaseDicom(SCase_idx).ManufacturerModelName = ManufacturerModelName;
SCaseDicom(SCase_idx).ProtocolName = ProtocolName;
SCaseDicom(SCase_idx).SeriesDescription = SeriesDescription;
end
fprintf(logFileId, '\n\nElapsed time (s): %s', num2str(toc)); % stopwatch timer
%% Write csv file
% This might be a function. The input would be a filename and a structure
% data
dataHeader = [...
'SCase_ID,' ...
'patient_IPP,' ...
'patient_gender,' ...
'patient_birthDate,' ...
'patient_age,' ...
'patient_height,' ...
'patient_weight,' ...
'CT_date,' ...
'CT_kernel,' ...
'CT_pixelSize,' ...
'CT_sliceSpacing,' ...
'CT_sliceThickness,' ...
'CT_tension,' ...
'CT_current,' ...
'CT_institution,' ...
'CT_manufacturer,' ...
'CT_model,' ...
'CT_protocol,' ...
'CT_SeriesDescription' ...
];
csvFilename = strcat(xlsDataDir, 'CTdicomInfo.csv');
% If dicominfoSCase is called with argument, we have to update the csv.
if update
% Read the csv and put data in a structure like SCaseDicom
% We should chech that file exists
csvTable = readtable(csvFilename,'DatetimeType', 'text');
% If 'DatetimeType' is specified as 'text', then the type for imported
% date and time data depends on the value specified in the 'TextType'
% parameter: If 'TextType' is 'char', then readtable returns dates as a
% cell array of character vectors. If 'TextType' is 'string', then
% readtable returns dates as an array of strings.
% Set csv data in temporary SCaseDicomTmp structure
SCaseDicom_csv = table2struct(csvTable);
% Update the SCaseDicomTmp structure with new data
for iSCaseDicom = 1:length(SCaseDicom) % loop on new SCase(s)
iSCaseDicomTmp = strcmp({SCaseDicom_csv.SCase_id}, SCaseDicom(iSCaseDicom).id); % 0/1 array
iSCaseDicomTmp = find(iSCaseDicomTmp); % index with non-zero, correponding to SCaseID
% We should test iSCaseDicomTmp
SCaseDicom_csv(iSCaseDicomTmp).patient_IPP = SCaseDicom(iSCaseDicom).PatientID;
SCaseDicom_csv(iSCaseDicomTmp).patient_gender = SCaseDicom(iSCaseDicom).PatientSex;
SCaseDicom_csv(iSCaseDicomTmp).patient_birthDate = SCaseDicom(iSCaseDicom).PatientBirthDate;
SCaseDicom_csv(iSCaseDicomTmp).patient_age = SCaseDicom(iSCaseDicom).PatientAge;
SCaseDicom_csv(iSCaseDicomTmp).patient_height = SCaseDicom(iSCaseDicom).PatientSize;
SCaseDicom_csv(iSCaseDicomTmp).patient_weight = SCaseDicom(iSCaseDicom).PatientWeight;
SCaseDicom_csv(iSCaseDicomTmp).CT_date = SCaseDicom(iSCaseDicom).AcquisitionDate;
SCaseDicom_csv(iSCaseDicomTmp).CT_kernel = SCaseDicom(iSCaseDicom).ConvolutionKernel;
SCaseDicom_csv(iSCaseDicomTmp).CT_pixelSize = SCaseDicom(iSCaseDicom).PixelSpacing;
SCaseDicom_csv(iSCaseDicomTmp).CT_sliceSpacing = SCaseDicom(iSCaseDicom).SpacingBetweenSlices;
SCaseDicom_csv(iSCaseDicomTmp).CT_sliceThickness = SCaseDicom(iSCaseDicom).SliceThickness;
SCaseDicom_csv(iSCaseDicomTmp).CT_tension = SCaseDicom(iSCaseDicom).KVP;
SCaseDicom_csv(iSCaseDicomTmp).CT_current = SCaseDicom(iSCaseDicom).XrayTubeCurrent;
SCaseDicom_csv(iSCaseDicomTmp).CT_institution = SCaseDicom(iSCaseDicom).InstitutionName;
SCaseDicom_csv(iSCaseDicomTmp).CT_manufacturer = SCaseDicom(iSCaseDicom).Manufacturer;
SCaseDicom_csv(iSCaseDicomTmp).CT_model = SCaseDicom(iSCaseDicom).ManufacturerModelName;
SCaseDicom_csv(iSCaseDicomTmp).CT_protocol = SCaseDicom(iSCaseDicom).ProtocolName;
SCaseDicom_csv(iSCaseDicomTmp).CT_SeriesDescription = SCaseDicom(iSCaseDicom).SeriesDescription;
end
% SCaseDicom must now be updated
for iSCaseDicom = 1:length(SCaseDicom)
SCaseDicom(iSCaseDicom).id = SCaseDicom_csv(iSCaseDicom).SCase_id;
SCaseDicom(iSCaseDicom).PatientID = SCaseDicom_csv(iSCaseDicom).patient_IPP;
SCaseDicom(iSCaseDicom).PatientSex = SCaseDicom_csv(iSCaseDicom).patient_gender;
SCaseDicom(iSCaseDicom).PatientBirthDate = SCaseDicom_csv(iSCaseDicom).patient_birthDate;
SCaseDicom(iSCaseDicom).PatientAge = SCaseDicom_csv(iSCaseDicom).patient_age;
SCaseDicom(iSCaseDicom).PatientSize = SCaseDicom_csv(iSCaseDicom).patient_height;
SCaseDicom(iSCaseDicom).PatientWeight = SCaseDicom_csv(iSCaseDicom).patient_weight;
SCaseDicom(iSCaseDicom).AcquisitionDate = SCaseDicom_csv(iSCaseDicom).CT_date;
SCaseDicom(iSCaseDicom).ConvolutionKernel = SCaseDicom_csv(iSCaseDicom).CT_kernel;
SCaseDicom(iSCaseDicom).PixelSpacing = SCaseDicom_csv(iSCaseDicom).CT_pixelSize;
SCaseDicom(iSCaseDicom).SpacingBetweenSlices = SCaseDicom_csv(iSCaseDicom).CT_sliceSpacing;
SCaseDicom(iSCaseDicom).SliceThickness = SCaseDicom_csv(iSCaseDicom).CT_sliceThickness;
SCaseDicom(iSCaseDicom).KVP = SCaseDicom_csv(iSCaseDicom).CT_tension;
SCaseDicom(iSCaseDicom).XrayTubeCurrent = SCaseDicom_csv(iSCaseDicom).CT_current;
SCaseDicom(iSCaseDicom).InstitutionName = SCaseDicom_csv(iSCaseDicom).CT_institution;
SCaseDicom(iSCaseDicom).Manufacturer = SCaseDicom_csv(iSCaseDicom).CT_manufacturer;
SCaseDicom(iSCaseDicom).ManufacturerModelName = SCaseDicom_csv(iSCaseDicom).CT_model;
SCaseDicom(iSCaseDicom).ProtocolName = SCaseDicom_csv(iSCaseDicom).CT_protocol;
SCaseDicom(iSCaseDicom).SeriesDescription = SCaseDicom_csv(iSCaseDicom).CT_SeriesDescription;
end
end
% Write csv file
csvFileId = fopen(csvFilename, 'w');
fprintf(csvFileId,dataHeader);
for SCase_idx=1:1:numel(SCaseDicom) % Loop on all SCase provided by listSCase()
fprintf(csvFileId,'\n%s,%s,%s,%s,%f,%f,%f,%s,%s,%f,%f,%f,%f,%f,%s,%s,%s,%s,%s',...
SCaseDicom(SCase_idx).id,...
SCaseDicom(SCase_idx).PatientID,...
SCaseDicom(SCase_idx).PatientSex,...
SCaseDicom(SCase_idx).PatientBirthDate,...
SCaseDicom(SCase_idx).PatientAge,...
SCaseDicom(SCase_idx).PatientSize,...
SCaseDicom(SCase_idx).PatientWeight,...
SCaseDicom(SCase_idx).AcquisitionDate,...
SCaseDicom(SCase_idx).ConvolutionKernel,...
SCaseDicom(SCase_idx).PixelSpacing,...
SCaseDicom(SCase_idx).SpacingBetweenSlices,...
SCaseDicom(SCase_idx).SliceThickness,...
SCaseDicom(SCase_idx).KVP,...
SCaseDicom(SCase_idx).XrayTubeCurrent,...
SCaseDicom(SCase_idx).InstitutionName,...
SCaseDicom(SCase_idx).Manufacturer,...
SCaseDicom(SCase_idx).ManufacturerModelName,...
SCaseDicom(SCase_idx).ProtocolName,...
SCaseDicom(SCase_idx).SeriesDescription...
);
end
fclose(csvFileId);
fprintf(logFileId, ['\n\nSaved csv file: ' csvFilename]);
%% Write xls file (Windows only)
[~,message] = csv2xlsSCase(csvFilename);
fprintf(logFileId, message);
%% Update MySQL
% To be done
%{
addpath('XLS_MySQL_DB');
MainExcel2SQL;
%}
%% Final log
fprintf(logFileId, '\n\nSCases treated: %s', num2str(numel(SCaseDicom))); % Number of SCases treated
fprintf(logFileId, '\nElapsed time (s): %s', num2str(toc)); % stopwatch timer
fclose(logFileId); % Close log file
outputArg = 1; % Might be extended
end
c4science · Help