Hi Tatyana,
I understand that you are facing errors when tranforming the code in 'ea_spherical_roi.m' of the Lead-DBS MATLAB toolbox to generate ellipses instead of spheres. The latest version of the toolbox does not seem to have the reference template image 't1.nii'. Moreover, executing the provided code for 'ea_elliptical_roi_test' resulted in a different error than the one specified in the question.
You can instead refer to the following working code to generate ellipses.
function roi = ea_elliptical_roi_test(fname, center, semiaxes, orientation, crop, ref, bg)
% Create ellipse ROI based on specified center, semiaxes, and orientation (all in mm)
% Write out NIfTI or not
if isempty(fname)
writeoutNii = 0;
else
writeoutNii = 1;
end
% Expand semiaxes in case multiple centers specified
if size(center,1) > 1
if size(semiaxes,1) == 1
semiaxes = repmat(semiaxes, size(center,1), 1);
elseif size(center,1) ~= size(semiaxes,1)
error('Number of centers doesn''t match number of semiaxes sets!');
end
end
% Crop the generated ROI image or not
if ~exist('crop','var')
crop = 1;
end
% Reference template image, use MNI t1 by default
if exist('ref','var')
ref = ea_load_nii(ref);
else
ref = ea_load_nii([ea_space,'t1.nii']);
end
% Preset background
if ~exist('bg','var')
ref.img(:) = 0;
else
if isscalar(bg)
ref.img(:) = bg;
elseif isequal(size(ref.img), size(bg))
ref.img = bg;
else
error('Background should be either a scalar value or of the same size as the reference image!');
end
end
voxsize = ref.voxsize;
dim = ref.dim;
for i = 1:size(center,1)
% mm to voxel conversion
c = ea_mm2vox(center(i,:), ref.mat);
a = semiaxes(i,:) ./ voxsize; % Semi-axis lengths in voxels
% Construct rotation matrix from orientation
R = ea_orient2mat(orientation(i,:));
% Calculate the maximum extent of the ellipse in each dimension
% after applying the rotation
rotated_extents = abs(R * diag(a));
max_extent = max(rotated_extents, [], 2)';
% Construct conservative bounding box limits
bboxlim = [max([1 1 1; ceil(c-max_extent)]); min([dim; floor(c+max_extent)])];
[xgrid, ygrid, zgrid] = meshgrid(bboxlim(1,1):bboxlim(2,1),...
bboxlim(1,2):bboxlim(2,2),...
bboxlim(1,3):bboxlim(2,3));
% Flatten voxel grid to x, y and z subscripts
xyz = [xgrid(:), ygrid(:), zgrid(:)];
% Apply rotation to the grid points
xyz_rot = (R * (xyz - c)')' + c;
% Find voxels within the ellipse
inside_ellipse = sum(((xyz_rot - c) ./ a) .^ 2, 2) <= 1;
ref.img(sub2ind(dim, xyz(inside_ellipse,1), xyz(inside_ellipse,2), xyz(inside_ellipse,3))) = 1;
end
% Adapt ROI NIfTI structure
ref.dt(1) = 16;
ref.fname = fname;
roi = ref;
% Write out NIfTI
if writeoutNii
ea_write_nii(ref);
% Crop ROI image
if crop
ea_autocrop(fname);
end
end
end
function R = ea_orient2mat(angle)
% Convert angle from degrees to radians
angle_rad = deg2rad(angle);
% Rotation matrix for rotation about the z-axis
R = [cos(angle_rad), -sin(angle_rad), 0;
sin(angle_rad), cos(angle_rad), 0;
0, 0, 1];
end
Please run the following script to generate the ellipse using the above code.
outputPath = 'ppX_left_ellipse.nii';
center = [-11.1, -10.24, -1.78];
radii = [12.5, 2, 4.8];
orientation = 45;
% Adding the path to a different reference image than the default t1.nii
ea_elliptical_roi_test(outputPath, center, radii, orientation, [ea_space,'brainmask.nii\brainmask.nii']);
Hope this helps with your query.
