decsg fails to decompose a few ellipses out of a set of 19 for unknown reason

Hello all, I am attempting to generate a decomposed geometry matrix of various ellipses in a square region using decsg in order to then extrude the geometry into 3D and mesh for use in a static structural model in the pde toolbox.

However, some of the ellipses (#4 when included with #s 1-3 or #19 when included with the rest) cause decsg to fail when included in the decsg inputs as so:

Error using decsg

Unable to decompose objects in the geometry description matrix into minimal regions.

Error in scatteredEllipses (line 42)

[dl, ~] = decsg(gd,sf,ns);

The scatteredEllipses function is presented below in its entirety, I have also attached the force_x, force_y, and force_z variables within forces.mat which are 512x512 double matrices of force magnitudes which represent desired ellipse locations.

clearvars -except force*; close all;
%% Fit ellipses to each traction island
tractionMask = abs(force_x) + abs(force_y) + abs(force_z); %combine x,y,z tractions to create mask of all nonzero elements
tractionMask(tractionMask ~= 0) = 1; %convert nonzero elements to 1
cc = bwconncomp(tractionMask); %get all connected ellipse regions
s = regionprops(cc,{'Centroid','Orientation','MajorAxisLength','MinorAxisLength'}); %fit ellipse parameters to each detected region
figure, imshow(tractionMask,[]), hold on
theta = linspace(0,2*pi); 
orientationList = zeros(size(s,1),2);
for i = 1:length(s)
    col = (s(i).MajorAxisLength/2)*cos(theta);
    row = (s(i).MinorAxisLength/2)*sin(theta);
    M = makehgtform('translate',[s(i).Centroid, 0],'zrotate',deg2rad(-1*s(i).Orientation));
    [orientationList(i,1),orientationList(i,2)] = pol2cart(deg2rad(-1*s(i).Orientation),1);
    D = M*[col;row;zeros(1,numel(row));ones(1,numel(row))];
    plot(D(1,:),D(2,:),'r','LineWidth',2); hold on
    s(i).Orientation = -1 * s(i).Orientation; %convert orientation to xy coordinate frame
end
%axis xy
hold off

%% Generate geometry
meshPtsFwdSol = 512;
halfSide = meshPtsFwdSol/2;
bound = [3; 4; -halfSide; halfSide; halfSide; -halfSide; halfSide; halfSide; -halfSide; -halfSide];
ns = char('bound');
sf = 'bound';

full = 1:length(s);
workingSubset = 4:length(s)-1;
g = 1;
for i = workingSubset %switch with full to cause error
    sf = [sf,'+elli',num2str(g)];
    ns = char(ns,['elli',num2str(g)]);
    elli(:,g) = [4; s(i).Centroid(1)-halfSide; s(i).Centroid(2)-halfSide; s(i).MajorAxisLength/2; s(i).MinorAxisLength/2; deg2rad(s(i).Orientation); 0; 0; 0; 0];
    g = g + 1;
end

gd = [bound,elli];
ns = ns';
[dl, ~] = decsg(gd,sf,ns);
%figure,pdegplot(dl,'EdgeLabels','on','FaceLabels','on')
pdem = createpde;
g = geometryFromEdges(pdem,dl);
%figure,pdegplot(pdem,'EdgeLabels','on','FaceLabels','on')
facets = facetAnalyticGeometry(pdem,g,0);
gm = analyticToDiscrete(facets);
pdem.Geometry = gm;
figure,pdegplot(pdem,'EdgeLabels','on','FaceLabels','on')

Answer 1

[Jaynik]

Hi Samuel,

After tinkering around with workingSubset, I found that ellipses #4 and #18 were causing the error when integrated with others and the code works for workingSubset = [1:3, 5;17, 19]. I started further analyzing #3 and #4 together by plotting the ellipses using the following code:

figure, hold on;
axis equal;
xlim([-halfSide, halfSide]);
ylim([-halfSide, halfSide]);
theta = linspace(0, 2*pi);
for idx = workingSubset
    a = s(idx).MajorAxisLength/2; % Semi-major axis
    b = s(idx).MinorAxisLength/2; % Semi-minor axis
    X0 = s(idx).Centroid(1)-halfSide; % Adjusted X centroid
    Y0 = s(idx).Centroid(2)-halfSide; % Adjusted Y centroid
    phi = deg2rad(s(idx).Orientation); % Orientation in radians
    
    % Parametric equation of the ellipse
    x = X0 + a*cos(theta)*cos(phi) - b*sin(theta)*sin(phi);
    y = Y0 + a*cos(theta)*sin(phi) + b*sin(theta)*cos(phi);
    
    plot(x, y, 'r', 'LineWidth', 2);
end

hold off;

On plotting the ellipses, I could see that the ellipses #3 and #4 are very close to each other. This could be a cause for the error of invalid geometry.

To resolve this:

• If possible, increase the distance between close geometries. If your geometries are very close to each other but not overlapping or intersecting, increasing the distance between them may help resolve the issue.
• You can also try simplyfying the geometries. This could involve reducing the number of points in your geometries, or simplifying the shapes themselves.

According to the documentation of decsg, specifically this part:

https://www.mathworks.com/help/pde/ug/decsg.html?s_tid=doc_ta#:~:text=Data%20Types%3A%20double-,Limitations,-In%20rare%20cases error can occur due to geometries being almost tangent.

Hope this help!