Skip to content
Permalink
Browse files

Using pre-calculated pyradiomics values to compare and test

  • Loading branch information...
rkpandya committed May 6, 2019
1 parent 50c79fc commit d2a2383b36dac5fb57d29974cb4ecccee0ac08c2
@@ -3,7 +3,7 @@
% RKP, 03/22/2018



%% Load image
firstOrderParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
'Unit_Testing','tests_for_cerr','test_first_order_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
@@ -34,36 +34,41 @@
firstOrderS.robustMeanAbsDev, firstOrderS.rms, firstOrderS.skewness, ...
firstOrderS.std, firstOrderS.var, firstOrderS.entropy];

%% Calculate features using pyradiomics

testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
mask3M = zeros(size(testM),'logical');
[rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
[maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
mask3M(:,:,uniqueSlices) = maskBoundBox3M;
dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
pixelSize = [dx dy dz]*10;
teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);

pyradFirstorderNamC = {'Energy', 'TotalEnergy','InterquartileRange','Kurtosis',...
'Maximum', 'Mean','MeanAbsoluteDeviation','Median','medianAbsDev',...
'Minimum','10Percentile','90Percentile','InterquartileRange',...
'RobustMeanAbsoluteDeviation','RootMeanSquared','Skewness',...
'StandardDeviation','Variance','Entropy'};

pyradFirstorderNamC = strcat(['wavelet', '_', dirString, '_firstorder_'],pyradFirstorderNamC);
%pyradFirstorderNamC = strcat(['original', '_firstorder_'],pyradFirstorderNamC);

pyRadFirstOrderV = [];
for i = 1:length(pyradFirstorderNamC)
if isfield(teststruct,pyradFirstorderNamC{i})
pyRadFirstOrderV(i) = teststruct.(pyradFirstorderNamC{i});
else
pyRadFirstOrderV(i) = NaN;
end
end
% %% Calculate features using pyradiomics
%
% testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
% single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
% mask3M = zeros(size(testM),'logical');
% [rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
% [maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
% mask3M(:,:,uniqueSlices) = maskBoundBox3M;
% dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
% pixelSize = [dx dy dz]*10;
% teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);
%
% pyradFirstorderNamC = {'Energy', 'TotalEnergy','InterquartileRange','Kurtosis',...
% 'Maximum', 'Mean','MeanAbsoluteDeviation','Median','medianAbsDev',...
% 'Minimum','10Percentile','90Percentile','InterquartileRange',...
% 'RobustMeanAbsoluteDeviation','RootMeanSquared','Skewness',...
% 'StandardDeviation','Variance','Entropy'};
%
% pyradFirstorderNamC = strcat(['wavelet', '_', dirString, '_firstorder_'],pyradFirstorderNamC);
% %pyradFirstorderNamC = strcat(['original', '_firstorder_'],pyradFirstorderNamC);
%
% pyRadFirstOrderV = [];
% for i = 1:length(pyradFirstorderNamC)
% if isfield(teststruct,pyradFirstorderNamC{i})
% pyRadFirstOrderV(i) = teststruct.(pyradFirstorderNamC{i});
% else
% pyRadFirstOrderV(i) = NaN;
% end
% end
%
% %% Compare
% diffFirstOrderV = (cerrFirstOrderV - pyRadFirstOrderV) ./ cerrFirstOrderV * 100

diffFirstOrderV = (cerrFirstOrderV - pyRadFirstOrderV) ./ cerrFirstOrderV * 100
%% Compare with previously calculated values of pyradiomics first order features
saved_pyRadFirstOrderV = [87498183.4756395,1049978201.70767,30.4772882461548,16.3069672735214,836.854187011719,-0.191955631352121,48.4171827693622,0.0908048748970032,NaN,-718.366882324219,-73.1680297851562,64.9023208618165,30.4772882461548,16.4187532424825,96.2495122333083,0.546840923205712,NaN,9263.93175818536,3.17432864851443];
diffFirstOrderV = (cerrFirstOrderV - saved_pyRadFirstOrderV) ./ cerrFirstOrderV * 100
@@ -2,7 +2,7 @@
%
% RKP, 03/22/2018


%% Load image
glcmParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
'Unit_Testing','tests_for_cerr','test_glcm_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
@@ -61,3 +61,7 @@

%% Compare
glcmDiffV = (cerrGlcmV - pyRadGlcmV) ./ cerrGlcmV * 100

%% Compare using previously calculated values of pyradiomics glcm
saved_pyRadGlcmV = [1751.49603382875,40.8299981390930,1699291.05899833,-11302.3472463506,399.183404315141,61.6492931961095,0.731298422981024,3.58225481213654,3.08032041137040,48.5207221130057,NaN,0.0325243174278277,6.85603719995355,0.522724816541355,0.480259668311807,-0.208253976228619,0.884543363584260,0.995164157047163,0.971313972188856,0.373754218660873,NaN,4.82370970733959,NaN];
glcmDiffV = (cerrGlcmV - saved_pyRadGlcmV) ./ cerrGlcmV * 100
@@ -2,7 +2,7 @@
%
% RKP, 03/22/2018


%% Load image
glcmParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
'Unit_Testing','tests_for_cerr','test_glcm_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
@@ -18,6 +18,7 @@

scanType = 'wavelet';
dirString = 'HHH';

%% Calculate features using CERR
harFeat3DdirS = calcGlobalRadiomicsFeatures...
(scanNum, strNum, paramS, planC);
@@ -28,37 +29,41 @@
harlCombS.secondInfCorr, harlCombS.invDiffMomNorm, harlCombS.invDiffNorm, harlCombS.invVar, ...
harlCombS.sumAvg, harlCombS.sumEntropy, harlCombS.sumVar];

%% Calculate features using pyradiomics
% image and mask for a structure
testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
mask3M = zeros(size(testM),'logical');
[rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
[maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
mask3M(:,:,uniqueSlices) = maskBoundBox3M;

dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
pixelSize = [dx dy dz]*10;

teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);

%teststruct = PyradWrapper(testM, mask3M, scanType, dirString);
pyradGlcmNamC = {'Autocorrelation', 'JointAverage', 'ClusterProminence', 'ClusterShade', 'ClusterTendency', ...
'Contrast', 'Correlation', 'DifferenceAverage', 'DifferenceEntropy', 'DifferenceVariance', 'Dissimilarity', ...
'JointEnergy', 'JointEntropy','Id','Idm', 'Imc1' , ...
'Imc2', 'Idmn','Idn','InverseVariance', 'sumAverage', 'SumEntropy', 'sumVariance'};

pyradGlcmNamC = strcat(['wavelet','_', dirString,'_glcm_'],pyradGlcmNamC);
pyRadGlcmV = [];
for i = 1:length(pyradGlcmNamC)
if isfield(teststruct,pyradGlcmNamC{i})
pyRadGlcmV(i) = teststruct.(pyradGlcmNamC{i});
else
pyRadGlcmV(i) = NaN;
end
end
% %% Calculate features using pyradiomics
% % image and mask for a structure
% testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
% single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
% mask3M = zeros(size(testM),'logical');
% [rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
% [maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
% mask3M(:,:,uniqueSlices) = maskBoundBox3M;
%
% dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
% pixelSize = [dx dy dz]*10;
%
% teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);
%
% %teststruct = PyradWrapper(testM, mask3M, scanType, dirString);
% pyradGlcmNamC = {'Autocorrelation', 'JointAverage', 'ClusterProminence', 'ClusterShade', 'ClusterTendency', ...
% 'Contrast', 'Correlation', 'DifferenceAverage', 'DifferenceEntropy', 'DifferenceVariance', 'Dissimilarity', ...
% 'JointEnergy', 'JointEntropy','Id','Idm', 'Imc1' , ...
% 'Imc2', 'Idmn','Idn','InverseVariance', 'sumAverage', 'SumEntropy', 'sumVariance'};
%
% pyradGlcmNamC = strcat(['wavelet','_', dirString,'_glcm_'],pyradGlcmNamC);
% pyRadGlcmV = [];
% for i = 1:length(pyradGlcmNamC)
% if isfield(teststruct,pyradGlcmNamC{i})
% pyRadGlcmV(i) = teststruct.(pyradGlcmNamC{i});
% else
% pyRadGlcmV(i) = NaN;
% end
% end
%
% %% Compare
% glcmDiffV = (cerrGlcmV - pyRadGlcmV) ./ cerrGlcmV * 100

%% Compare
glcmDiffV = (cerrGlcmV - pyRadGlcmV) ./ cerrGlcmV * 100
%% Compare using previously calculated values of pyradiomics glcm
saved_pyRadGlcmV = [929.516284841440,30.4950632056943,5316.83199378276,-0.303806680516285,20.9971321089971,22.7275310589439,-0.0380659828547803,2.78099962843313,3.04970358686672,14.7880425011987,NaN,0.0446949242970363,6.28336538646016,0.479852104826365,0.424210035556516,-0.0921069324860780,0.664903434331990,0.993049866792659,0.955511453133062,0.367836063144979,NaN,3.75942560575858,NaN];
glcmDiffV = (cerrGlcmV - saved_pyRadGlcmV) ./ cerrGlcmV * 100
@@ -41,64 +41,3 @@
subplot(1,3,2), imagesc(cerrLog3M(:,:,slc)), title('CERR recursive LOG')
subplot(1,3,3), imagesc(log3M(:,:,slc)), title('ITK recursive LOG')


%
% %% CERR LoG
% %scanArray3M = flip(onesM,3);
% kernelSize = 5;
% sigmaVoxelV = [4 4 1]; %[0.79 0.79 6.5];
% filtScan3M = LoGFilt(testM,kernelSize,sigmaVoxelV);
%
% %% Analysis and comparison between Pyradiomics and CERR
% maxMat = max(filtScan3M(:));
% maxPy = max(pyFiltM(:));
% diff3M = (filtScan3M/maxMat - pyFiltM/maxPy);
% quantile(abs(diff3M(:)),0.95)
%
% figure('Name', 'diff3M'), imagesc(diff3M(:,:,13))
%
% slcPyM = pyFiltM(:,:,5);
% slcMatM = filtScan3M(:,:,21);
% figure, imagesc(slcPyM), title('PyRad Gauss')
% figure, imagesc(slcMatM), title('CERR Gauss')
%
% lap3M = -ones(3,3,3)/26;
% lap3M(2,2,2) = 1;
% lofFiltScan3M = convn(filtScan3M,lap3M);
% figure, imagesc(lofFiltScan3M(:,:,21)), title('CERR LOG')
%
% lofPyrad3M = convn(pyFiltM,-lap3M);
% figure, imagesc(lofPyrad3M(:,:,5)), title('Pyrad LOG')
%
%
% %figure, imagesc(slcMatM/64.2026)
% figure('Name', 'slcPyM 13'), imagesc(slcPyM/1.6886e+03)
% figure('Name', 'slcMatM 13'), imagesc(flipud(slcMatM)/64.2026)
% diffM = slcPyM/1.6886e+03 - flipud(slcMatM)/64.2026;
% figure('Name', 'slc diff'), imagesc(diffM)
%
% % figure, imagesc(slcMatM/64.2026)
% % figure, imagesc(slcPyM/1.6886e+03)
% % figure, imagesc(flipud(slcMatM)/64.2026)
% % figure, imagesc(slcPyM/1.6886e+03)
% % diffM = slcPyM/1.6886e+03 - flipud(slcMatM)/64.2026;
% % figure, imagesc(diffM)
%
%
% %
% % %Matlab Gaussian fileter result:
% % matlabGaussian = imgaussfilt3(testM,4);
% % validM2 = matlabGaussian(6:end-6,6:end-6,6:end-6);
% % validM2./validPyM
% %
% % %matlabGaussian./pydata
% % gaussDiffV2 = (matlabGaussian - pydata) ./ matlabGaussian * 100
% %
% %
% % filename = 'C:\Users\pandyar1\AppData\Local\Temp\log-sigma-2-0-mm-3D.nrrd';
% % data3M = nrrdread(filename);
% % data3M = flipdim(data3M,3);
% % figure; imagesc(data3M(:,:,20))
% % figure; imagesc(filtScan3M(:,:,2))
%
%
@@ -2,7 +2,7 @@
%
% RKP, 03/22/2018


%% load image
gldmParamFileName = fullfile(fileparts(fileparts(getCERRPath)),...
'Unit_Testing','tests_for_cerr','test_ngldm_radiomics_extraction_settings.json');
cerrFileName = fullfile(fileparts(fileparts(getCERRPath)),...
@@ -26,46 +26,49 @@
ngldmS.gln, ngldmS.glnNorm, ngldmS.dcn, ngldmS.dcnNorm,...
ngldmS.dcp, ngldmS.glv, ngldmS.dcv, ngldmS.entropy, ngldmS.energy];

%% Calculate features using pyradiomics

testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
mask3M = zeros(size(testM),'logical');
[rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
[maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
mask3M(:,:,uniqueSlices) = maskBoundBox3M;

scanType = 'original';

dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
pixelSize = [dx dy dz]*10;

teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);

%teststruct = PyradWrapper(testM, mask3M, scanType);

pyradNgldmNamC = {'SmallDependenceEmphasis', 'LargeDependenceEmphasis',...
'LowGrayLevelCountEmphasis', 'HighGrayLevelCountEmphasis', 'SmallDependenceLowGrayLevelEmphasis', ...
'SmallDependenceHighGrayLevelEmphasis', 'LargeDependenceLowGrayLevelEmphasis', ...
'LargeDependenceHighGrayLevelEmphasis', 'GrayLevelNonUniformity', 'GrayLevelNonUniformityNorm', ...
'DependenceNonUniformity', 'DependenceNonUniformityNormalized', ...
'DependencePercentage', 'GrayLevelVariance', 'DependenceVariance', ...
'DependenceEntropy', 'DependenceEnergy'};


pyradNgldmNamC = strcat(['original', '_gldm_'],pyradNgldmNamC);

pyRadNgldmV = [];
for i = 1:length(pyradNgldmNamC)
if isfield(teststruct,pyradNgldmNamC{i})
pyRadNgldmV(i) = teststruct.(pyradNgldmNamC{i});
else
pyRadNgldmV(i) = NaN;
end
end

%% Compare

ngldmDiffV = (cerrNgldmV - pyRadNgldmV) ./ cerrNgldmV * 100
% %% Calculate features using pyradiomics
%
% testM = single(planC{indexS.scan}(scanNum).scanArray) - ...
% single(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
% mask3M = zeros(size(testM),'logical');
% [rasterSegments, planC, isError] = getRasterSegments(strNum,planC);
% [maskBoundBox3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
% mask3M(:,:,uniqueSlices) = maskBoundBox3M;
%
% scanType = 'original';
%
% dx = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dy = planC{indexS.scan}(scanNum).scanInfo(1).grid1Units;
% dz = mode(diff([planC{indexS.scan}(scanNum).scanInfo(:).zValue]));
% pixelSize = [dx dy dz]*10;
%
% teststruct = PyradWrapper(testM, mask3M, pixelSize, scanType, dirString);
%
% %teststruct = PyradWrapper(testM, mask3M, scanType);
%
% pyradNgldmNamC = {'SmallDependenceEmphasis', 'LargeDependenceEmphasis',...
% 'LowGrayLevelCountEmphasis', 'HighGrayLevelCountEmphasis', 'SmallDependenceLowGrayLevelEmphasis', ...
% 'SmallDependenceHighGrayLevelEmphasis', 'LargeDependenceLowGrayLevelEmphasis', ...
% 'LargeDependenceHighGrayLevelEmphasis', 'GrayLevelNonUniformity', 'GrayLevelNonUniformityNorm', ...
% 'DependenceNonUniformity', 'DependenceNonUniformityNormalized', ...
% 'DependencePercentage', 'GrayLevelVariance', 'DependenceVariance', ...
% 'DependenceEntropy', 'DependenceEnergy'};
%
%
% pyradNgldmNamC = strcat(['original', '_gldm_'],pyradNgldmNamC);
%
% pyRadNgldmV = [];
% for i = 1:length(pyradNgldmNamC)
% if isfield(teststruct,pyradNgldmNamC{i})
% pyRadNgldmV(i) = teststruct.(pyradNgldmNamC{i});
% else
% pyRadNgldmV(i) = NaN;
% end
% end
%
% %% Compare
% ngldmDiffV = (cerrNgldmV - pyRadNgldmV) ./ cerrNgldmV * 100

%% Compare with previously calculated pyradiomics features
saved_pyRadNgldmV = [0.177896163786327,80.7514028586554,NaN,NaN,0.00117423243903225,380.447741305586,3.84549682517831,136232.842773954,1250.78930651138,NaN,634.668395976707,0.0671962303839817,NaN,122.726429080550,26.6319748251348,7.44897904303236,NaN];
ngldmDiffV = (cerrNgldmV - saved_pyRadNgldmV) ./ cerrNgldmV * 100
Oops, something went wrong.

0 comments on commit d2a2383

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