Skip to content
Permalink
Browse files

Merge pull request #348 from aditiiyer/master

Function to compute texture maps from DICOM data and export results to DICOM
  • Loading branch information...
aditiiyer committed Oct 8, 2019
2 parents d7720eb + 6bf571a commit 2c52feddea6287cb55cdad1908698d6a958a0cc8
@@ -1,4 +1,4 @@
function importDICOM(source,destn)
function planC = importDICOM(source,destn)
% function importDICOM(source,destn)
%
% Function to import DICOM from source and write CERR file to destn.
@@ -118,20 +118,24 @@ function importDICOM(source,destn)
%fullOutFileName = fullfile(destn,fileName);

%Check for duplicate name of fullOutFileName
dirOut = dir(destn);
allOutNames = {dirOut.name};
if any(strcmpi([outFileName,'.mat'],allOutNames))
fullOutFileName = [outFileName,'_duplicate_',num2str(rand(1))];
end
if strcmpi(zipFlag,'Yes')
saved_fullFileName = fullfile(destn,[outFileName,'.mat.bz2']);
else
saved_fullFileName = fullfile(destn,[outFileName,'.mat']);
end
if ~exist(fileparts(saved_fullFileName),'dir')
mkdir(fileparts(saved_fullFileName))
if exist('destn','var') && ~isempty(destn)
dirOut = dir(destn);
allOutNames = {dirOut.name};
if any(strcmpi([outFileName,'.mat'],allOutNames))
fullOutFileName = [outFileName,'_duplicate_',num2str(rand(1))];
end
if strcmpi(zipFlag,'Yes')
saved_fullFileName = fullfile(destn,[outFileName,'.mat.bz2']);
else
saved_fullFileName = fullfile(destn,[outFileName,'.mat']);
end


if ~exist(fileparts(saved_fullFileName),'dir')
mkdir(fileparts(saved_fullFileName))
end
save_planC(planC,[], 'passed', saved_fullFileName);
end
save_planC(planC,[], 'passed', saved_fullFileName);

catch

@@ -0,0 +1,93 @@
function generateTextureMapFromDICOM(inputDicomPath,strNameC,configFilePath,outputDicomPath)
% generateTextureMapFromDICOM(inputDicomPath,strNameC,configFilePath,outputDicomPath);
%
% Compute texture maps from DICOM images and export to DICOM.
% -------------------------------------------------------------------------
% INPUTS
% inputDicomPath : Path to DICOM data.
% strNameC : List of structure names.
% A mask of the bounding box around these structures is
% used to compute texture maps.
% configFilePath : Path to config files for texture calculation.
% outputDicomPath : Output directory.
% -------------------------------------------------------------------------
% AI 10/8/19


%% Import DICOM data
planC = importDICOM(inputDicomPath);

%% Get structure mask
indexS = planC{end};
strC = {planC{indexS.structures}.structureName};
strNumV = nan(1,length(strNameC));
for n = 1:length(strNameC)
strNumV(n) = getMatchingIndex(strNameC{n},strC,'EXACT');
end
mask3M = getStrMask(strNumV,planC);

%% Get scan array
scanNum = getStructureAssociatedScan(strNumV(1),planC);
scan3M = getScanArray(scanNum,planC);
CToffset = planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
scan3M = scan3M - CToffset;

%% Read config file
paramS = getRadiomicsParamTemplate(configFilePath);

%% Apply filter
filterType = fieldnames(paramS.imageType);
filterType = filterType{1};
paramS = paramS.imageType.(filterType);
outS = processImage(filterType,scan3M,mask3M,paramS);
fieldName = fieldnames(outS);
fieldName = fieldName{1};
filtScan3M = outS.(fieldName);


%% Create Texture Scans
assocScanUID = planC{indexS.scan}(scanNum).scanUID;
nTexture = length(planC{indexS.texture}) + 1;
planC{indexS.texture}(nTexture).assocScanUID = assocScanUID;
assocStrUID = strjoin({planC{indexS.structures}(strNumV).strUID},',');
planC{indexS.texture}(nTexture).assocStructUID = assocStrUID;
planC{indexS.texture}(nTexture).category = filterType;
planC{indexS.texture}(nTexture).parameters = paramS;
planC{indexS.texture}(nTexture).description = filterType;
planC{indexS.texture}(nTexture).textureUID = createUID('TEXTURE');

[xVals, yVals, zVals] = getScanXYZVals(planC{indexS.scan}(scanNum));
dx = abs(mean(diff(xVals)));
dy = abs(mean(diff(yVals)));
dz = abs(mean(diff(zVals)));
deltaXYZv = [dy dx dz];
zV = zVals;
regParamsS.horizontalGridInterval = deltaXYZv(1);
regParamsS.verticalGridInterval = deltaXYZv(2); %(-)ve for dose
regParamsS.coord1OFFirstPoint = xVals(1);
regParamsS.coord2OFFirstPoint = yVals(end);
regParamsS.zValues = zV;
regParamsS.sliceThickness =[planC{indexS.scan}(scanNum).scanInfo(:).sliceThickness];
assocTextureUID = planC{indexS.texture}(nTexture).textureUID;

planC = scan2CERR(filtScan3M,filterType,'Passed',regParamsS,assocTextureUID,planC);


%% Write filtered image to DICOM
for n = 1:length(planC{indexS.scan})-1
planC = deleteScan(planC, 1);
end

newSeriesInstanceUID = dicomuid;
for n = 1:length(planC{indexS.scan}(1).scanInfo)
planC{indexS.scan}(1).scanInfo(n).seriesInstanceUID = newSeriesInstanceUID;
end

planC = generate_DICOM_UID_Relationships(planC);
if ~exist(outputDicomPath,'dir')
mkdir(outputDicomPath)
end
export_CT_IOD(planC,outputDicomPath,filterType);


end
@@ -10,6 +10,10 @@
%-------------------------------------------------------------------------
%AI 03/16/18

if ~exist('hWait','var')
hWait = [];
end

filterType = strrep(filterType,' ','');

% record the original image size
@@ -53,7 +57,7 @@
paramS.binWidth.val = [];
end

if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
[energy,entropy,sumAvg,corr,...
invDiffMom,contrast,clustShade,...
clustProminence,haralCorr] = textureByPatchCombineCooccur(volToEval,...
@@ -132,7 +136,7 @@

outS.(outname) = out3M;

if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 (n-1)/(length(dirListC)-1) (n-1)/(length(dirListC)-1)]' [0 1 1 0]']);
drawnow;
end
@@ -147,7 +151,7 @@
out3M = out3M(:,:,1:end-1);
end
out3M = flip(out3M,3);
if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 1 1]' [0 1 1 0]']);
drawnow;
end
@@ -161,7 +165,7 @@
scan3M = scan3M(minr:maxr,minc:maxc,mins:maxs);
vol3M = double(mask3M).*double(scan3M);
[outS.SobelMag,outS.SobelDir] = sobelFilt(vol3M);
if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 1 1]' [0 1 1 0]']);
drawnow;
end
@@ -176,7 +180,7 @@
scan3M = scan3M(minr:maxr,minc:maxc,mins:maxs);
vol3M = double(mask3M).*double(scan3M);
outS.LoG_recursive = recursiveLOG(vol3M,paramS.Sigma_mm.val,paramS.VoxelSize_mm.val);
if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 1 1]' [0 1 1 0]']);
drawnow;
end
@@ -187,7 +191,7 @@
vol3M = double(mask3M).*double(scan3M);
outS.Gabor = filtImgGabor(vol3M,paramS.Radius.val,paramS.Sigma.val,...
paramS.AspectRatio.val,paramS.Orientation.val,paramS.Wavlength.val);
if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 1 1]' [0 1 1 0]']);
drawnow;
end
@@ -214,7 +218,7 @@
outV = patchStatM(:,n);
out3M(mask3M) = outV;
outS.(statC{n}) = out3M;
if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 n/length(statC) n/length(statC)]' [0 1 1 0]']);
drawnow;
end
@@ -244,7 +248,7 @@
text3M = convn(paddedVolM,lawsMasksS.(fieldNamesC{i}),'same');
text3M = text3M(6:end-5,6:end-5,6:end-5);
outS.(fieldNamesC{i}) = text3M; % for the entire cubic roi
if exist('hWait','var') && ishandle(hWait)
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 i/numFeatures i/numFeatures]' [0 1 1 0]']);
drawnow;
end
@@ -259,6 +263,10 @@
coLlAGe3M = getCollageFeature(scan3M, mask3M, paramS.Dominant_Dir_Radius.val,...
paramS.Cooccur_Radius.val, paramS.Number_Gray_Levels.val, dir, hWait);
outS.entropy = coLlAGe3M;
if ishandle(hWait)
set(hWait, 'Vertices', [[0 0 1 1]' [0 1 1 0]']);
drawnow;
end

otherwise

@@ -140,10 +140,14 @@
attr.addAll(SOPattr);

clear imgattr imgplaneattr imgpixelattr ctimageattr SOPattr

fileNum = num2str(filenumber + nWritten);
filename = fullfile(destDirPath,['IMG_', repmat('0', [1 5-length(fileNum)]), fileNum]);


if ischar(filenumber)
filename = fullfile(filenameRoot,filenumber);
else
fileNum = num2str(filenumber + nWritten);
filename = fullfile(destDirPath,['IMG_', repmat('0', [1 5-length(fileNum)]), fileNum]);
end

writefile_mldcm(attr, filename);

nWritten = nWritten + 1;

0 comments on commit 2c52fed

Please sign in to comment.
You can’t perform that action at this time.