How can I transform the code of the ea_spherical_roi function to one that creates ellipses?

Dear all,

I was trying to transform the ea_spherical_roi function to one that creates ellipses. The original function is:

function roi = ea_spherical_roi(fname,center,radius,crop,ref,bg)

% Create sphere ROI based on specified center and radius (both in mm)

% Write out NIfTI or not

if isempty(fname)

writeoutNii = 0;

else

writeoutNii = 1;

end

% Expand radius in case multiple centers specified

if size(center,1)>1

if length(radius)==1

radius = repmat(radius, 1, size(center,1));

elseif size(center,1)~=length(radius)

error('Length of centers doesn''t match length of radius!');

end

end

% Crop the generate 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);

r = radius(i)./voxsize;

% Construct voxel grid for the sphere of cencter c and radius r

bboxlim = [max([1 1 1; ceil(c-r)]); min([dim; floor(c+r)])];

[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(:)];

% Find voxels within the sphere

xyz = xyz(sqrt(sum(((xyz - c) .* voxsize) .^ 2, 2)) <= radius(i), :);

ref.img(sub2ind(dim, xyz(:,1), xyz(:,2), xyz(:,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

My tranformed function is:

function roi = ea_elliptical_roi_test(fname, center, radii, orientation, crop, ref, bg)

% Create elliptical ROI based on specified center, radii, and orientation

% Write out NIfTI or not

if isempty(fname)

writeoutNii = 0;

else

writeoutNii = 1;

end

% Expand radii in case multiple centers specified

if size(center, 1) > 1

if size(radii, 1) == 1

radii = repmat(radii, size(center, 1), 1);

elseif size(center, 1) ~= size(radii, 1)

error('Length of centers doesn''t match length of radii!');

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

r = radii(i, :)./voxsize;

% Construct voxel grid for the ellipse of center c, radii r, and orientation

[xgrid, ygrid, zgrid] = meshgrid(1:dim(1), 1:dim(2), 1:dim(3));

xyz = rotate_coordinates([xgrid(:), ygrid(:), zgrid(:)], center(i, :), orientation);

xyz = xyz(sqrt((xyz(:, 1) - c(1)).^2 / r(1)^2 + (xyz(:, 2) - c(2)).^2 / r(2)^2 + (xyz(:, 3) - c(3)).^2 / r(1)^2) <= 1, :);

ref.img(sub2ind(dim, xyz(:, 1), xyz(:, 2), xyz(:, 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 rotatedCoords = rotate_coordinates(coords, center, orientation)

% Rotate 3D coordinates around the center

% Translate to the origin

coords(:, 1:3) = coords(:, 1:3) - center;

% Convert orientation to radians

orientationRad = deg2rad(orientation);

% Perform 2D rotation in the XY plane

X_rot = coords(:, 1) * cos(orientationRad) - coords(:, 2) * sin(orientationRad);

Y_rot = coords(:, 1) * sin(orientationRad) + coords(:, 2) * cos(orientationRad);

% Translate back to the original position

rotatedCoords = [X_rot + center(1), Y_rot + center(2), coords(:, 3) + center(3)];

end

I then plug in:

outputPath = '\path\to\output\ppX_left_ellipse.nii';

center = [-11.1, -10.24, -1.78];

radii = [12.5, 2, 4.8];

orientation = 45;

% Run the script

ea_elliptical_roi_test(outputPath, center, radii, orientation);

However, I keep getting error messages:

Arrays have incompatible sizes for this operation.

Error in ea_elliptical_roi_test (line 51)

r = radii(i, :)./voxsize;

Error in testscriptellipse (line 7)

ea_elliptical_roi_test(outputPath, center, radii, orientation);

Related documentation

Could someone help me and tell me whats wrong please?

Answer 1

[Paras Gupta]

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.