function prepareCtSeries(indir, outdir)
% Usage:copy and paste the line belwo in MATLAB command line; change the
% directory and name as needed
%   prepareCtSeries('C:\WiscPhotonkV\ratDICOMs', 'C:\WiscPhotonkV\ratDICOMs_Downsampled');

% outdir = 'C:\WiscPhotonkV\ratDICOMs_2Downsampled';
% indir = 'C:\WiscPhotonkV\ratDICOMs';


ctvol = readCtSeries(indir);
% reset orientation (dirty fix)
ctvol.info.ImageOrientationPatient = [1;0;0;0;1;0];
ctvol.info.PatientPosition = 'HFS';
dsvol = simpleDownsampleCtVol(ctvol, [1 1 1]);
if size(dsvol.data,1) > size(dsvol.data,2)
    fprintf('Padding unequal dimension CT\n');
    dsvol.data = padarray(dsvol.data, [0, size(dsvol.data,1) - size(dsvol.data,2), 0], -1000, 'post');
elseif size(dsvol.data,1) < size(dsvol.data,2)
    fprintf('Padding unequal dimension CT\n');
    dsvol.data = padarray(dsvol.data, [size(dsvol.data,2) - size(dsvol.data,1), 0, 0], -1000, 'post');
end
writeCtVol(dsvol, outdir)

function s = readCtSlice(filename)
s.info = dicominfo(filename);
s.data = single(dicomread(filename)) * s.info.RescaleSlope + s.info.RescaleIntercept;

function ctvol = readCtSeries(indir)
filelist = dir(indir);
slices = arrayfun(@(x)readCtSlice(fullfile(indir, x.name)), filelist(3:end));
[~, ix] = sort(arrayfun(@(x)x.info.ImagePositionPatient(3), slices));
slices = slices(ix);

ctvol.info = slices(1).info;
ctvol.info.SliceThickness = slices(2).info.ImagePositionPatient(3) - slices(1).info.ImagePositionPatient(3);
ctvol.data = cell2mat(reshape({slices.data}, [1 1 numel(slices)]));

% ctinfs = arrayfun(@(x)dicominfo(fullfile(indir, x.name)), filelist(3:end), 'UniformOutput', false);
% ctimgs = arrayfun(@(x)dicomread(fullfile(indir, x.name)), filelist(3:end), 'UniformOutput', false);
% [~, ix] = sort(cellfun(@(x)x.ImagePositionPatient(3), ctinfs));
% ctinfs = ctinfs(ix);
% ctimgs = ctimgs(ix);

% ctvol.info = ctinfs{1};
% ctvol.info.SliceThickness = ctinfs{2}.ImagePositionPatient(3) - ctinfs{1}.ImagePositionPatient(3);
% ctvol.data = cell2mat(reshape(ctimgs, [1 1 numel(ctimgs)]));

function dsvol = simpleDownsampleCtVol(ctvol, dsfactor)
meanFilter = ones(dsfactor) / prod(dsfactor);
meanImg = imfilter(ctvol.data, meanFilter, 'same');
dsvol.data = meanImg(1:dsfactor(1):end, 1:dsfactor(2):end, 1:dsfactor(3):end);

dsvol.info = ctvol.info;
dsvol.info.PixelSpacing = ctvol.info.PixelSpacing .* reshape(dsfactor(1:2), [], 1);
dsvol.info.SliceThickness = ctvol.info.SliceThickness * dsfactor(3);

function writeCtVol(ctvol, outdir)
mkdir(outdir);
% changes common to all slices
ctvol.info.SeriesInstanceUID = dicomuid;
ctvol.info.RescaleIntercept = min(ctvol.data(:));
ctvol.info.RescaleSlope = (max(ctvol.data(:)) - min(ctvol.data(:))) / 65536;
ctvol.data = uint16((ctvol.data - ctvol.info.RescaleIntercept) / ctvol.info.RescaleSlope);
% ctvol.info
for k = 1:size(ctvol.data, 3)
    info = rmfield(ctvol.info, 'SliceThickness');
    info.SOPInstanceUID = dicomuid;
    info.ImagePositionPatient(3) = ctvol.info.ImagePositionPatient(3) + (k-1) * ctvol.info.SliceThickness;
    dicomwrite(ctvol.data(:,:,k), fullfile(outdir, sprintf('CT.%s.%d.dcm', info.SeriesInstanceUID, k)), info, 'CreateMode', 'copy');
end