Prevenire il patch se l'area è già occupata

Question

Questa è una sfida abbastanza difficile che ho con il mio codice. Prima di tutto il codice che sto mettendo qui non è eseguibile perché sto usando un foglio Excel (ma sono felice di mandarlo via email se la gente vuole provare ad usare il mio codice).

Quello che ho è un foglio Excel con dati sulle fibre trasversali in un'immagine microscopica che ho scattato. L'informazione è fondamentalmente: location della sezione, area, angle di rotazione.

Da questo calcolo l'angolo di orientamento Phi e Gamma. Successivamente utilizzo la funzione scatter per tracciare un punto di colori diversi per ciascun valore di angolo Phi. Uso un colore costante per un intervallo di 10 gradi. Il che mi dà una foto come questa:

Ora il mio obiettivo di è calcolare l'area di ogni regione omogenea. Quindi cerco un modo per tracciare diciamo tutti i punti all'interno della regione -10 +10 (sto facendo 20 gradi per ora, ma ne farò 10 dopo). Ho usato uno sguardo e ottengo una foto come questa:

I corrisponde bianchi in cui i punti sono all'interno della gamma, turisti. Dopo di ciò, uso la toolbox in MATLAB per convertire ogni punto in un pixel. Quindi otterrò uno sfondo nero con un sacco di pixel bianchi, quindi uso imdilate per creare cerchi, riempire buchi e isolare ogni regione con un colore specifico. Finalmente utilizzo le funzioni boundary e patch, per creare ogni contorno e riempirli con un colore. E ho una foto come questa:

Che è quello che voglio e posso ottenere l'area di ogni regione e l'area totale (ho usato una soglia per scartare le piccole aree). Quindi eseguo il codice più volte per ogni regione, e io uso imfuse per rimetterli insieme e vedere come appare.

IL PROBLEMA è, si sovrappongono molto, e questo perché ci sono alcuni errori nei miei dati, e quindi alcuni punti blu saranno in rosso e così via. Quindi voglio eseguire il codice una volta, poi quando lo eseguo nuovamente con un altro intervallo, fa la stessa cosa ma non prende in considerazione il valore quando c'è già qualcosa tracciato prima.

Ho provato a farlo, dopo aver eseguito una volta, salvando la matrice bw4 e aggiungendo una condizione quando si traccia la foto in bianco e nero, dicendo se Phi è nel mio raggio E non c'è bianco qui, quindi puoi mettere bianco, altrimenti è nero. Ma non sembra funzionare.

Capisco che è una cosa abbastanza complicata da spiegare, ma apprezzerei qualsiasi idea e aprirò alla chat via email o in altro modo. Sto inserendo il codice completo ora e posso inviarti il ​​mio foglio Excel se vuoi eseguirlo sul tuo computer e vedere di persona.

clearvars -except data colheaders bw4 
close all 
clc 
%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% CHANGE DATA FOR EACH SAMPLE %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

cd 'C:\Users\dkarta\Desktop\Sample 12\12.6' 
data=xlsread('Sample12_6res.xlsx'); 

cd 'C:\Users\dkarta\Documents\MATLAB' 

%data=Sample121res; % Data name 
imax=length(data); % Numbers of rows in data sheet 
y=11900; % Number of pixels in the y on image j 
%% 

data(:,15)=data(:,9)*pi/180; % Convers Column 9 (angle of rotation) in rads 
data(:,16)=y-data(:,6); % Reset the Y coordinate axis to bottom left 
delta = 0 : 0.01 : 2*pi; % Angle in paramteric equations 
theta=45*pi/180; % Sample cutting angle in rads 

%AA=[data(:,5)' data(:,16)' phi'] 
% Define colors 
beta=acos(data(1:imax,8)./data(1:imax,7));%./acos(0); 
phi=atan(sin(beta).*cos(data(1:imax,15))./(sin(theta)*sin(beta).*sin(data(1:imax,15))+cos(theta)*cos(beta)))/(pi/2); 
phi2=phi/2+1/2; % Scales in plane angle phi between 0 and 1 
gamma=atan((cos(theta)*sin(beta).*sin(data(1:imax,15))-sin(theta)*cos(beta))./... 
    (sin(theta)*sin(beta).*sin(data(1:imax,15))+cos(theta)*cos(beta)))/(pi/2); 
gamma2=gamma+1/2; % Scales out of plane angle gamma between 0 and 1 
%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% MESHGRID AND COLOURMAP %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%% 
x1=data(1:imax,5); 
y1=data(1:imax,16); 
z1=phi*90; 
z2=gamma*90; 
n=300; 
%Create regular grid across data space 
[X,Y] = meshgrid(linspace(min(x1),max(x1),n), linspace(min(y1),max(y1),n)); 

% Creating a colormap with 10 degree constant colors 
map4=[0 0 1;0 1/3 1;0 2/3 1; 0 1 1;0 1 2/3;0 1 1/3;0 1 0;1/3 1 0;2/3 1 0;1 1 0;1 0.75 0;1 0.5 0;1 0.25 0;1 0 0;0.75 0 0.25;0.5 0 0.5;0.25 0 0.75; 0 0 1]; 
Colormap4=colormap(map4); 
h=colorbar; 
caxis([-90 90]) 
set(h, 'YTick', [-90:10:90]) 
%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% PLOT USING SCATTER - ISOLATE SOME REGIONS %%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
a=-10; % Lower boundary for angle interval 
b=10; % Upper boundary for angle interval 
c=z1>a & z1 < b; 
c=c.*1; 
%j=1; 

y1=(y1-min(y1)+1); 
y2=max(y1)-y1+1; 
[X1,Y1]=meshgrid(1:500,1:500); 
griddata(x1,y2,c,X1,Y1); 
clear c1 


for i=1:imax 
    if z1(i)< b && z1(i)> a %&& bw4(round(y1(i)),round(x1(i))) == 0 
     c(i) = 1; 
     c1(round(y2(i)),round(x1(i)))=1; 
    else 
     c(i)= 0; 
     c1(round(y2(i)),round(x1(i)))=0; 
    end 

end 
C=[c c c]; 

%c(find(c==0)) = NaN; 
%contourf(X,Y,griddata(x1,y1,c,X,Y),100,'EdgeColor', 'None') 
figure(1), scatter(x1,y1,3,z1,'filled'); 
axis equal 
axis ([0 8000 0 12000]) 
axis off 
figure(2), scatter(x1,y1,3,C,'filled'); 
axis equal 
axis ([0 8000 0 12000]) 
axis off 


se=strel('disk',50,8); 
bw2=imdilate(c1,se); 
bw4=bwlabel(bw2); 
bw3=imfill(bw4,'holes'); 
max(bw4(:)); 
figure(3),imshow(c1,'InitialMagnification', 10); 
figure(4), imshow(bw2,'InitialMagnification', 10); 
figure(5), imshow(bw3,'InitialMagnification', 10); 
figure(6),imshow(label2rgb(bw4),'InitialMagnification', 10); 

k=ones(max(bw4(:)),1); 
clear bw5 
for i=1:length(x1) 
    if bw3(round(y2(i)),round(x1(i))) ~= 0 
     m=bw3(round(y2(i)),round(x1(i))); 
     bw5{m}(k(m),1)=x1(i); bw5{m}(k(m),2)=y2(i); 
     k(m)=k(m)+1; 
    end 
end 

figure(7), imshow(~c1,'InitialMagnification', 10); 
hold on 
for i=1:max(bw4(:)) 
    %scatter(bw5{i}(:,1),bw5{i}(:,2)) 
    j = boundary(bw5{i}(:,1),bw5{i}(:,2),0.5); 
    %poly=convhull(bw5{i}(:,1),bw5{i}(:,2)); 
    %plot(bw5{i}(poly,1),bw5{i}(poly,2)), title('convhull') 
    if polyarea(bw5{i}(j,1),bw5{i}(j,2))> 10^5; 
     patch(bw5{i}(j,1),bw5{i}(j,2),'r'), title('boundary') 
     indexminy(i)=find(min(bw5{i}(:,2)) == bw5{i}(:,2)); 
     indexminx(i)=find(min(bw5{i}(:,1)) == bw5{i}(:,1)); 
     indexmaxy(i)=find(max(bw5{i}(:,2)) == bw5{i}(:,2)); 
     indexmaxx(i)=find(max(bw5{i}(:,1)) == bw5{i}(:,1)); 
     %xmin = bw5{i}(indexminx); xmax = bw5{i}(indexmaxx); 
     %ymin = bw5{i}(indexminy); ymax = bw5{i}(indexmaxy); 
     str=[(indexminx(i)+indexmaxx(i))/2,(indexminy(i)+indexmaxy(i))/2,'Region no.',num2str(i)]; 
     text((min(x1(i))+max(x1(i)))/2,(min(y1(i))+max(y1(i)))/2,str) 
     polya(i)=polyarea(bw5{i}(j,1),bw5{i}(j,2)); 
    end 
end 
spolya=sum(polya(:)) 
print -dpng -r500 B 

solo per mostrarvi più foto di quando ho fondermi alcuni di loro:

E quando ho fusibile:

Come si può vedere si sovrappongono, che non voglio, quindi voglio ogni immagine che creo per 'sapere' quello che sto facendo nelle esecuzioni precedenti in modo che non si sovrapponga. Voglio ottenere l'area percentuale di ogni regione e se si sovrappongono non posso usare l'area totale effettiva del mio campione ei risultati sono sbagliati.

Answer 1

Non ho il mio MATLAB funzionante, ma qui è quello che devi fare.

Per la prima corsa fanno una matrice di zeri uguale alla dimensione dell'immagine

already_taken = zeros(size(bw3)); 

Poi ogni esecuzione, è possibile riempire le regioni adottate da questa iterazione. Così, alla fine del codice, in cui si salva l'output in un png, leggere di nuovo in qualcosa di simile

this_png = rgb2gray(imread(current_png_path))>threshold; 

convertire questo in una matrice logica facendo qualche thresholding e aggiungere questi valori in già preso. Così, alla fine del codice, fare un

already_taken = already_taken | this_png; % You might need to check if you need a single | or a double || 

Così ora avete un'immagine di pixel già prese, malato cuocere che io non permetto BW2 di prendere questi valori al primo posto

bw2(already_taken) = 0; 

E alla fine del codice, quando voglio scrivere il mio png, la mia creazione di confine intelligente potrebbe essere di nuovo entrata nell'area già_taken, così ancora una volta dovrò mettere un po 'di controllo. Per quanto ho capito, questo limite viene creato in base al tuo bw5. Quindi, se mai riempi questa matrice, prova a mettere un controllo simile come ho fatto sopra per bw2.

Spero che questo aiuti.