Add Nikon uCT data loader

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'));
 

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

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

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 [email protected].');
+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');
+