function RtomoAll = tomoBeamOnOffSinogram(GeomROI,xpmin,xpmax,Mxp,Nphi,Nrot)

% Figures out which beams should be on or off for tomotherapy treatment.


ptv = [];
ptv.x_count = GeomROI.dim(1);
ptv.y_count = GeomROI.dim(2);
ptv.z_count = GeomROI.dim(3);
ptv.non_zero_indices = GeomROI.ind;
ptv.non_zero_values = single(GeomROI);

ptvFull = full3D(ptv);

tumorMask = [];
tumorMask.start = GeomROI.start;
tumorMask.voxel_size = GeomROI.pixdim;
tumorMask.data = ptvFull;

[M,N,Q] = size(tumorMask.data);

% create axes
x = [0:M-1]*tumorMask.voxel_size(1) + tumorMask.start(1);
y = [0:N-1]*tumorMask.voxel_size(2) + tumorMask.start(2);
z = [0:Q-1]*tumorMask.voxel_size(3) + tumorMask.start(3);

% Calculate the maximum tumor cross section
bigSlice = double(~~sum(tumorMask.data,3));

% grow the tumor by 5 voxels
BW = bwdist(bigSlice);
bigSlice(BW < 5) = 1;

bigSlice = rot90(bigSlice);

% Specify beam parameters:  must match those in helicalDosecalcSetup.m
% exactly!
% xpmin = -20;
% xpmax = 20;
% Mxp = 64;

xpTomo = [-(Mxp-1)/2:(Mxp-1)/2]*(xpmax - xpmin)/Mxp;

% Nphi = 51;
phiTomo = [0:Nphi-1]*360/Nphi;

phi = [0:359];
% calculate the radon transform of the mask
[R,xp] = radon(bigSlice,phi);
xp = xp*tumorMask.voxel_size(1);

% create downsampling matrices
[PHI,XP] = meshgrid(phi,xp);
[PHITOMO,XPTOMO] = meshgrid(phiTomo,xpTomo);

% resample the sinogram
Rinterp = interp2(PHI,XP,R,PHITOMO,XPTOMO);
Rtomo = ~~interp2(PHI,XP,R,PHITOMO,XPTOMO);

% make copies of sinogram for 3D dose calculation

RtomoDouble = double(Rtomo);

% convert 2D sinograms to 3D for dose calc

RtomoDouble = reshape(RtomoDouble,Mxp*Nrot,1);
rowVector = ones(1,Nrot);
RtomoAll = RtomoDouble*rowVector;
RtomoAll = reshape(RtomoAll,Mxp,Nphi*Nrot);

fprintf('Unused beams turned off.\n')