From 0298d7bcab059cf48f6026038e85db8bab7895bd Mon Sep 17 00:00:00 2001 From: Ander Biguri Date: Mon, 30 Nov 2020 13:06:39 +0000 Subject: [PATCH] Add Nikon uCT data loader --- MATLAB/InitTIGRE.m | 1 + MATLAB/Utilities/IO/Nikon/NikonDataLoader.m | 28 ++++ .../Utilities/IO/Nikon/loadNikonProjections.m | 144 ++++++++++++++++++ .../Utilities/IO/Nikon/readXtekctGeometry.m | 102 +++++++++++++ 4 files changed, 275 insertions(+) create mode 100644 MATLAB/Utilities/IO/Nikon/NikonDataLoader.m create mode 100644 MATLAB/Utilities/IO/Nikon/loadNikonProjections.m create mode 100644 MATLAB/Utilities/IO/Nikon/readXtekctGeometry.m diff --git a/MATLAB/InitTIGRE.m b/MATLAB/InitTIGRE.m index 0e63fa00..08e214f7 100644 --- a/MATLAB/InitTIGRE.m +++ b/MATLAB/InitTIGRE.m @@ -23,6 +23,7 @@ addpath('./Utilities'); addpath('./Utilities/Quality_measures'); addpath('./Utilities/IO/VarianCBCT'); +addpath('./Utilities/IO/Nikon'); addpath(genpath('./Test_data')); diff --git a/MATLAB/Utilities/IO/Nikon/NikonDataLoader.m b/MATLAB/Utilities/IO/Nikon/NikonDataLoader.m new file mode 100644 index 00000000..b32684ad --- /dev/null +++ b/MATLAB/Utilities/IO/Nikon/NikonDataLoader.m @@ -0,0 +1,28 @@ +function [proj,geo,angles]=NikonDataLoader(filepath,varargin) +%NikonDataLoader(filepath) Loads Nikon uCT datasets into TIGRE standard +% +% NikonDataLoader(filepath OPT,VAL,...) uses options and values. +% These are options in case you don't want to load the entire +% dataset, but only particular sets of projections. +% The possible options in OPT are: +% 'sampling': type of sampling. default 'equidistant' Can be: +% 'equidistant': equidistantly sample the entire set +% of angles. 'num_angles' should be set for partial +% loading of the data. +% 'step': sample the entire set of projections every +% 'sampling_step' angles. +% 'continous': Load the first 'num_angles' amount of +% angles only. +% +% 'num_angles': Number of total angles to load. Default all of +% them. Useful for 'equidistant' and 'continous' loading +% +% 'sampling_step': step to load when loading projections. +% Default=1. Useful for 'step' loading. +%% + +[geo,angles]=readXtekctGeometry(filepath); +[proj,geo,angles]=loadNikonProjections(filepath,geo,angles,varargin{:}); + + +end \ No newline at end of file diff --git a/MATLAB/Utilities/IO/Nikon/loadNikonProjections.m b/MATLAB/Utilities/IO/Nikon/loadNikonProjections.m new file mode 100644 index 00000000..78890de0 --- /dev/null +++ b/MATLAB/Utilities/IO/Nikon/loadNikonProjections.m @@ -0,0 +1,144 @@ +function [proj,geo,angles]=loadNikonProjections(filepath,geo,angles,varargin) +%[proj,angles]=loadNikonProjections(filepath,geo,angles,varargin) +% loads Nikon uCT machine projections +% +% loadNikonData(filepath,geo,angles) Loads a dataset given its FILEPATH, +% a geometry GEO and angles ANGLES (loaded from readXtekctGeometry()) +% +% loadNikonData(filepath,geo,angles, OPT,VAL,...) uses options and values. +% These are options in case you don't want to load the entire +% dataset, but only particular sets of projections. +% The possible options in OPT are: +% 'sampling': type of sampling. default 'equidistant' Can be: +% 'equidistant': equidistantly sample the entire set +% of angles. 'num_angles' should be set for partial +% loading of the data. +% 'step': sample the entire set of projections every +% 'sampling_step' angles. +% 'continous': Load the first 'num_angles' amount of +% angles only. +% +% 'num_angles': Number of total angles to load. Default all of +% them. Useful for 'equidistant' and 'continous' loading +% +% 'sampling_step': step to load when loading projections. Default +% 1. Useful for 'step' loading. +% + +% developed by A. Biguri and W. Sun 06.07.2020 + + +%% Parse inputs. + +[angles_to_load,index]=parse_inputs(geo,angles,varargin); + +%% get filename +% assuming TIF and 4 digits. +firstfile = dir([filepath,'/*.tif']); % +firstfile = firstfile(1).name; +filename = firstfile(1:end-8); + +fprintf("Dataset in: %s \n", filepath); + +%% load images +%proj=[]; +l = length(angles_to_load); +proj = single(zeros(geo.nDetector(1),geo.nDetector(2),l)); +for ii=1:length(angles_to_load) + if(~mod(ii,50)) + fprintf("Loading: %d / %d \n", ii, length(angles_to_load)); + end + proj(:,:,ii)=single(imread([filepath,filename,num2str(index(ii),'%04d'),'.tif'])); +end +%% Beer lambert +if any(proj(:))>single(geo.whitelevel) + warning('Changing the Whitelevel value as projection data has higher value than specified in the file'); + geo.whitelevel=max(proj(:)+1); +end +proj=-log(proj/single(geo.whitelevel)); +geo=rmfield(geo,'whitelevel'); +%% +angles=angles_to_load; +end + +function [angles,indices]=parse_inputs(geo,angles,argin) +opts= {'sampling','num_angles','sampling_step'}; +defaults=ones(length(opts),1); +% Check inputs +nVarargs = length(argin); +if mod(nVarargs,2) + error('Invalid number of inputs') +end + +% check if option has been passed as input +for ii=1:2:nVarargs + ind=find(ismember(opts,lower(argin{ii}))); + if ~isempty(ind) + defaults(ind)=0; + else + error(['Optional parameter "' argin{ii} '" does not exist' ]); + end +end + +for ii=1:length(opts) + opt=opts{ii}; + default=defaults(ii); + % if one option is not default, then extract value from input + if default==0 + ind=double.empty(0,1);jj=1; + while isempty(ind) + ind=find(isequal(opt,lower(argin{jj}))); + jj=jj+1; + end + if isempty(ind) + error(['Optional parameter "' argin{jj} '" does not exist' ]); + end + val=argin{jj}; + end + + switch opt + case 'sampling' + if default + sampling='equidistant'; + else + if strcmpi(val,'equidistant') || strcmpi(val,'step')|| strcmpi(val,'continuous') + sampling=val; + else + error('sampling type not understood, must be equidistant, step or continuous'); + end + end + case 'num_angles' + if default + num_angles=length(angles); + else + num_angles=val; + end + case 'sampling_step' + if default + sampling_step=1; + else + sampling_step=val; + end + otherwise + error(['Invalid input name:', num2str(opt),'\n No such option']); + end +end +% now lets build the angles + +if strcmpi(sampling,'equidistant') + sampling_step=round(length(angles)/num_angles); + indices=1:sampling_step:length(angles); + angles=angles(indices); + return; +end +if strcmpi(sampling,'continuous') + indices=1:num_angles; + angles=angles(indices); + return; +end +if strcmpi(sampling,'step') + indices=1:sampling_step:length(angles); + angles=angles(indices); + return; +end +end \ No newline at end of file diff --git a/MATLAB/Utilities/IO/Nikon/readXtekctGeometry.m b/MATLAB/Utilities/IO/Nikon/readXtekctGeometry.m new file mode 100644 index 00000000..65270310 --- /dev/null +++ b/MATLAB/Utilities/IO/Nikon/readXtekctGeometry.m @@ -0,0 +1,102 @@ +function [geo,angles]=readXtekctGeometry(folder_or_file) + +% Developed by A. Biguri and W. Sun +% W. Sun edited on 06.10.2018 for 3D pro version 5.2.6809.15380 (2018) + +if endsWith(folder_or_file,'.xtekct') + % is the xtekct file itself + fid=fopen(folder_or_file); + isfolder=false; +else + % is the folder where it lives + isfolder=true; + file = dir([folder_or_file,'/*.xtekct']); % + if isempty(file) + error(['No .xtekct file found in folder: ', folder_or_file]); + end + filename=file(1).name; + fid=fopen([folder_or_file,'/',filename]); +end +% check if it was oppened right +if(fid==-1) + error('Wrong file path'); +end +xtekctText = textscan(fid, '%s %s', 'Delimiter', '=', 'HeaderLines', 1, 'CommentStyle', '['); +fclose(fid); + +%% Detector information +% Number of pixel in the detector +geo.nDetector=[str2double(xtekctText{2}(strcmp('DetectorPixelsX', xtekctText{1}))); + str2double(xtekctText{2}(strcmp('DetectorPixelsY', xtekctText{1})))]; +% Size of pixels in the detector +geo.dDetector=[str2double(xtekctText{2}(strcmp('DetectorPixelSizeX', xtekctText{1}))); + str2double(xtekctText{2}(strcmp('DetectorPixelSizeY', xtekctText{1})))]; +% Total size of the detector +geo.sDetector=geo.nDetector.*geo.dDetector; + +%% Offset of the detector: +geo.offDetector=[str2double(xtekctText{2}(strcmp('DetectorOffsetX', xtekctText{1}))); + str2double(xtekctText{2}(strcmp('DetectorOffsetY', xtekctText{1}))) ]; + +%% Image information +% Number of pixel in the detector +geo.nVoxel=[str2double(xtekctText{2}(strcmp('VoxelsX', xtekctText{1}))); + str2double(xtekctText{2}(strcmp('VoxelsY', xtekctText{1}))); + str2double(xtekctText{2}(strcmp('VoxelsZ', xtekctText{1}))) ]; +% Size of each pixel +geo.dVoxel=[str2double(xtekctText{2}(strcmp('VoxelSizeX', xtekctText{1}))); + str2double(xtekctText{2}(strcmp('VoxelSizeY', xtekctText{1}))); + str2double(xtekctText{2}(strcmp('VoxelSizeZ', xtekctText{1}))) ]; +% Size of the image in mm +geo.sVoxel=geo.nVoxel.*geo.dVoxel; +geo.offOrigin=[0;0;0]; +%% Global geometry +geo.DSO=str2double(xtekctText{2}(strcmp('SrcToObject', xtekctText{1}))); +geo.DSD=str2double(xtekctText{2}(strcmp('SrcToDetector', xtekctText{1}))); +geo.COR=-str2double(xtekctText{2}(strcmp('CentreOfRotationTop', xtekctText{1}))); +if (geo.COR==0) + warning('Centre of Rotation seems to be zero. Make sure that it is true and that the machine did not omit that information.'); +else + warning('TIGRE doesnt know if the sign of COR is the right one. Consider triying both and reporting to tigre.toolbox@gmail.com.'); +end +%% whitelevel + +geo.whitelevel=str2double(xtekctText{2}(strcmp('WhiteLevel', xtekctText{1}))); + +%% +%% angles +% It can be either an .ang or .txt file +% .ang +if ~isfolder + % rewrite it. + [folder_or_file,~,~] = fileparts(folder_or_file); +end + + +file_angles=dir([folder_or_file,'/','*.ang']); +if ~isempty(file_angles) + fid=fopen([folder_or_file,'/',file_angles.name]); + xtekctText = textscan(fid, '%s %s', 'Delimiter', '\t', 'HeaderLines', 1); + fclose(fid); + angles=str2double(xtekctText{2})*pi/180; + angles=angles.'; + return +end +% ._ctdata.txt +file_angles=dir([folder_or_file,'/','_ctdata.txt']); +if ~isempty(file_angles) + fid=fopen([folder_or_file,'/',file_angles.name]); + xtekctText = textscan(fid, '%s %s %s', 'Delimiter', '\t', 'HeaderLines', 3); + fclose(fid); + angles=str2double(xtekctText{2})*pi/180; + angles=angles.'; + return +end + +warning('File with definition of angles not found, estimating them from geometri info.') +anlge_step=str2double(xtekctText{2}(strcmp('AngularStep', xtekctText{1}))); +initial_angle=str2double(xtekctText{2}(strcmp('InitialAngle', xtekctText{1}))); +n_angles=str2double(xtekctText{2}(strcmp('Projections', xtekctText{1}))); +angles=(initial_angle:anlge_step:(initial_angle-anlge_step+360))*pi/180; +assert(size(angles,2)==n_angles,'Assertion failed: Inconsistent data detected. Number of projections and angle information do not match\n'); +