% Stephen R Bowen April 2009
% 
% Crops and resamples the planning CT data set and associated ROI masks,
% typically imported from Pinnacle using pinn2matlab.m.  CT data is
% linearly resampled and ROI masks are nearest-neighbor resampled so as to
% preserve binary values.  Bounding box coordinates are specified in the
% crop ranges user inputs [xmin xmax], [ymin ymax], [zmin zmax].  Be aware
% of half-voxel shifts in the boundary box used by Amira and remain
% consistent in defining your crop ranges.

%%%%% Start of user supplied inputs %%%%%

% 
infile = '../plan/HN003/HN003.mat';
outfile = '../plan/HN003/HN003final.mat';

% downsamples patient Geometry files, including their tissue masks

% crop ranges (if empty, no cropping is done)
xCrop = [-10.695 12.289];
yCrop = [-66.202 -39.682];
zCrop = [-6.801 8.669];

% resolution after resampling (if empty, resolution is unchanged
dxNew = 0.442;
dyNew = 0.442;
dzNew = 0.442;

x0New = -10.695;
y0New = -66.202;
z0New = -6.801;

%%%%% End of user supplied inputs %%%%%

load(infile);

% get dimensions of original data
[M,N,Q] = size(Geometry.data);

dx = Geometry.voxel_size(1);
dy = Geometry.voxel_size(2);
dz = Geometry.voxel_size(3);

x0 = Geometry.start(1);
y0 = Geometry.start(2);
z0 = Geometry.start(3);

% x,y,z coordinates of original data
x = single([0:M-1]*dx + x0);
y = single([0:N-1]*dy + y0);
z = single([0:Q-1]*dz + z0);

% If no cropping is desired, the resampling will be performed such that the
% bounding planes of the 3D grid data are unchanged.

% bounding plane coordinates:
xp1 = x(1) - dx/2;
xpM = x(end) + dx/2;
yp1 = y(1) - dy/2;
ypN = y(end) + dy/2;
zp1 = z(1) - dz/2;
zpQ = z(end) + dz/2;

% create the new coordinate system
if ~isempty(xCrop)
    xp1 = xCrop(1);
    xpM = xCrop(2);
end

if ~isempty(yCrop)
    yp1 = yCrop(1);
    ypN = yCrop(2);
end

if ~isempty(zCrop)
    zp1 = zCrop(1);
    zpQ = zCrop(2);
end

if isempty(dxNew)
    dxNew = dx;
end

if isempty(dyNew)
    dyNew = dy;
end

if isempty(dzNew)
    dzNew = dz;
end

% create new coordinate systems
if isempty(x0New)
    xNew = single(xp1+dxNew/2:dxNew:xpM-dxNew/2);
else
    xNew = x0New:dxNew:xpM;
    xNew = xNew(xNew >= xp1 & xNew <= xpM);
end

if isempty(y0New)
    yNew = single(yp1+dyNew/2:dyNew:ypN-dyNew/2);
else
    yNew = y0New:dyNew:ypN;
    yNew = yNew(yNew >= yp1 & yNew <= ypN);
end

if isempty(z0New)
    zNew = single(zp1+dzNew/2:dzNew:zpQ-dzNew/2);
else
    zNew = z0New:dzNew:zpQ;
    zNew = zNew(zNew >= zp1 & zNew <= zpQ);
end

clear Geometry_ds;
Geometry_ds.voxel_size = [dxNew dyNew dzNew];
Geometry_ds.start = [xNew(1) yNew(1) zNew(1)];

fprintf('Resampling CT data:\n');
Geometry_ds.data = gridResample3D(x,y,z,Geometry.data,xNew,yNew,zNew,'linear');
fprintf('Finished resampling CT data.\n');

Geometry_ds.ROIS = Geometry.ROIS;

% resample the ROIs
for k=1:length(Geometry.ROIS)    
    % expand the original ROI mask
    roimask = zeros(M,N,Q,'int8');  % create the original mask
    roimask(Geometry.ROIS{k}.ind) = 1;

    % resample the ROI mask
    fprintf('Resampling ROI mask: %s',Geometry.ROIS{k}.name);
    roimaskNew = gridResample3D(x,y,z,roimask,xNew,yNew,zNew,'nearest');
    fprintf('  Finished resampling %s.\n',Geometry.ROIS{k}.name);
    
    % save the new roi mask data
    Geometry_ds.ROIS{k}.ind = int32(find(roimaskNew));
    Geometry_ds.ROIS{k}.dim = size(Geometry_ds.data);
    Geometry_ds.ROIS{k}.pixdim = [dxNew dyNew dzNew];
    Geometry_ds.ROIS{k}.start = [xNew(1) yNew(1) zNew(1)];
    clear roimask roimaskNew;
end

if isfield(Geometry,'Trial');
    Geometry_ds.Trial = Geometry.Trial;
end

Geometry = Geometry_ds;

save(outfile,'Geometry');