How to remove open contours

7 visualizaciones (últimos 30 días)
Daniel
Daniel el 25 de Sept. de 2020
Respondida: Daniel el 29 de Sept. de 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.

Respuesta aceptada

Daniel
Daniel el 29 de Sept. de 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

Más respuestas (0)

Categorías

Más información sobre Contour Plots en Help Center y File Exchange.

Productos


Versión

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by