How to remove open contours

3 views (last 30 days)
Daniel
Daniel on 25 Sep 2020
Answered: Daniel on 29 Sep 2020
I have a volume of data. At each x location, I'd like to process some contours in the y-z plane. As part of that, I'd like to only consider the closed contours. I'm trying to do this by filtering out any contours that touch the edge of the plane (they could technically be tangent, but I'm willing to risk that). Here's my code with more explanations in comments.
for i = 1:length(xinds) % looping through x locations
yzconlev{i} = contourc(zdim,ydim,squeeze(uavgdef(xinds(i),:,:))',levs); % creating contours
[yzczdim{i},yzcydim{i},yzcuavgdef{i}] = C2xyz(yzconlev{i}); % Using C2xyz to get interpolated
% grid data of each contour
for j = 1:length(yzcuavgdef{i}) % looping through contours at each x
coniny{i}{j} = find(yzcydim{i}{j}~=ydim(1) & yzcydim{i}{j}~=ydim(end)); % find contours
% that touch the edge as defined by ydim
coninz{i}{j} = find(yzczdim{i}{j}~=zdim(1) & yzczdim{i}{j}~=zdim(end)); % same, but for
% other dimension
if size(coniny{i}{j}) == size(yzcydim{i}{j})
iny{i}{j} = yzcydim{i}{j};
end
if size(coninz{i}{j}) == size(yzczdim{i}{j})
inz{i}{j} = yzczdim{i}{j};
end
if isempty(iny{i}{j}) == 0 && isempty(inz{i}{j}) == 0
iny2{i}{j} = iny{i}{j}; % can't index with j if I want to skip some...
inz2{i}{j} = inz{i}{j};
end
% would then continue processing only with closed contours
end
end
This is close. I get only closed contours in iny and inz, but it leaves empty cells where the open contours were, and now I need to find the intersection of the two. I have the start of a solution, but I'm struggling to figure out how to index it so it works. As is, I still get empty cells as a result of skipping some js based on the condition.

Accepted Answer

Daniel
Daniel on 29 Sep 2020
Here's what I finally figured out. I'm certain that improvements could be made as I'm not an "efficient" coder. I welcome input.
for i = 1:length(xinds) % loop through the planes
k = 1;
yzconlev{i} = contourc(zdim,ydim,squeeze(uavgdef(xinds(i),:,:))',levs); % create contours
[yzczdim{i},yzcydim{i},yzcuavgdef{i}] = C2xyz(yzconlev{i});
for j = 1:length(yzcuavgdef{i}) % loop through the contour levels in each plane
conin{i}{j} = find(yzcydim{i}{j}~=ydim(1) & yzcydim{i}{j}~=ydim(end)...
& yzczdim{i}{j}~=zdim(1) & yzczdim{i}{j}~=zdim(end)); % find contours that don't touch
% edges of domain to exlude open contours
if size(conin{i}{j}) == size(yzcydim{i}{j})
if size(conin{i}{j}) == size(yzczdim{i}{j})
iny{i}{k} = yzcydim{i}{j}; % create sets of y and z coords for only open contours
inz{i}{k} = yzczdim{i}{j};
k = k+1;
end
end
end
for j = 1:length(iny{i})
concenter{i}(j,:) = [mean(inz{i}{j}),mean(iny{i}{j})]; % find center of open contours
end
incen{i} = find(concenter{i}(:,1)<1 & concenter{i}(:,1)>-1 ...
& concenter{i}(:,2)<1 & concenter{i}(:,2)>-1); % find centers that are within "middle"
inconcen{i} = concenter{i}(incen{i},:); % subset of centers in "middle"
for j = 1:length(incen{i})
inycen{i}{j} = iny{i}{incen{i}(j)}; % subset of contours with centers in "middle"
inzcen{i}{j} = inz{i}{incen{i}(j)};
conradii{i}{j} = sqrt((inconcen{i}(j,1)-inzcen{i}{j}).^2+...
(inconcen{i}(j,2)-inycen{i}{j}).^2); % R, radii of contours with centers in "middle"
conavgradius{i}(j) = mean(conradii{i}{j}); % average radius of each contour
end
gt1{i} = find(conavgradius{i} >= 1); % find average radii greater than one
radgtone{i} = conavgradius{i}(gt1{i}); % subset of average radii greater than one
if isempty(radgtone{i}) == 1
radgtone{i} = NaN; % need this for max to work below
end
avgradgtone(i) = nanmean(radgtone{i});
medradgtone(i) = nanmedian(radgtone{i});
maxradgtone(i) = max(radgtone{i});
inconcengt1{i} = inconcen{i}(gt1{i},:); % subset of centers with average radii greater than one
for j = 1:length(gt1{i})
inycengt1{i}{j} = inycen{i}{gt1{i}(j)}; % subset of contours with average radii greater than one
inzcengt1{i}{j} = inzcen{i}{gt1{i}(j)};
end
end

More Answers (0)

Categories

Find more on Contour Plots in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!