matlab - Ajuste del círculo más grande en el área libre en la imagen con partícula distribuida




math image-processing (4)

Estoy trabajando en imágenes para detectar y ajustar el círculo más grande posible en cualquiera de las áreas libres de una imagen que contiene partículas distribuidas:

(capaz de detectar la ubicación de la partícula).

Una dirección es definir un círculo tocando cualquier combinación de 3 puntos, verificando si el círculo está vacío y luego encontrando el círculo más grande entre todos los círculos vacíos. Sin embargo, conduce a una gran cantidad de combinación, es decir, C(n,3) , donde n es el número total de partículas en la imagen.

Agradecería si alguien me puede proporcionar alguna pista o método alternativo que pueda explorar.


¡Hagamos algunas matemáticas, amigo, ya que las matemáticas siempre llegarán al final!

Wikipedia:

En matemáticas, un diagrama de Voronoi es una partición de un plano en regiones basadas en la distancia a puntos en un subconjunto específico del plano.

Por ejemplo:

rng(1)
x=rand(1,100)*5;
y=rand(1,100)*5;


voronoi(x,y);

Lo bueno de este diagrama es que si notas que todos los bordes / vértices de esas áreas azules están a la misma distancia de los puntos a su alrededor. Por lo tanto, si conocemos la ubicación de los vértices y calculamos las distancias a los puntos más cercanos, entonces podemos elegir el vértice con la distancia más alta como nuestro centro del círculo.

Curiosamente, los bordes de las regiones de Voronoi también se definen como los circuncentros de los triángulos generados por una triangulación de Delaunay.

Entonces, si calculamos la triangulación de Delaunay del área y sus circuncentros

dt=delaunayTriangulation([x;y].');
cc=circumcenter(dt); %voronoi edges

Y calcule las distancias entre los circuncentros y cualquiera de los puntos que definen cada triángulo:

for ii=1:size(cc,1)
    if cc(ii,1)>0 && cc(ii,1)<5 && cc(ii,2)>0 && cc(ii,2)<5
    point=dt.Points(dt.ConnectivityList(ii,1),:); %the first one, or any other (they are the same distance)
    distance(ii)=sqrt((cc(ii,1)-point(1)).^2+(cc(ii,2)-point(2)).^2);
    end
end

Luego tenemos el centro ( cc ) y el radio ( distance ) de todos los círculos posibles que no tienen punto dentro de ellos. ¡Solo necesitamos el más grande!

[r,ind]=max(distance); %Tada!

Ahora vamos a trazar

hold on

ang=0:0.01:2*pi; 
xp=r*cos(ang);
yp=r*sin(ang);

point=cc(ind,:);

voronoi(x,y)
triplot(dt,'color','r','linestyle',':')
plot(point(1)+xp,point(2)+yp,'k');
plot(point(1),point(2),'g.','markersize',20);

Observe cómo el centro del círculo está en un vértice del diagrama de Voronoi.

NOTA : esto encontrará el centro dentro de [0-5], [0-5]. puede modificarlo fácilmente para cambiar esta restricción. También puede intentar encontrar el círculo que encaje en su totalidad dentro del área interesada (en lugar de solo el centro). Esto requeriría una pequeña adición al final donde se obtiene el máximo.


El hecho de que este problema se pueda resolver mediante una "búsqueda directa" (como se puede ver en otra respuesta ) significa que uno puede ver esto como un problema de optimización global . Existen varias formas de resolver tales problemas, cada una apropiada para ciertos escenarios. Por curiosidad personal, decidí resolver esto usando un algoritmo genético .

En términos generales, dicho algoritmo requiere que pensemos en la solución como un conjunto de "genes" sujetos a "evolución" bajo una determinada "función de aptitud". Como sucede, es bastante fácil identificar los genes y la función de aptitud en este problema:

  • Genes: x , y , r .
  • Función de aptitud: técnicamente, área máxima del círculo, pero esto es equivalente al máximo r (o mínimo -r , ya que el algoritmo requiere una función para minimizar ).
  • Restricción especial: si r es mayor que la distancia euclidiana al punto más cercano (es decir, el círculo contiene un punto), el organismo "muere".

A continuación se muestra una implementación básica de dicho algoritmo (" básico " porque no está optimizado por completo, y hay mucho espacio para la optimización sin juego de palabras en este problema).

function [x,y,r] = q42806059b(cloudOfPoints)
% Problem setup
if nargin == 0
  rng(1)
  cloudOfPoints = rand(100,2)*5; % equivalent to Ander's initialization.
end
%{
figure(); plot(cloudOfPoints(:,1),cloudOfPoints(:,2),'.w'); hold on; axis square;
set(gca,'Color','k'); plot(0.7832,2.0694,'ro'); plot(0.7832,2.0694,'r*');
%}
nVariables = 3;
options = optimoptions(@ga,'UseVectorized',true,'CreationFcn',@gacreationuniform,...
                       'PopulationSize',1000);

S = max(cloudOfPoints,[],1); L = min(cloudOfPoints,[],1); % Find geometric bounds:
% In R2017a: use [S,L] = bounds(cloudOfPoints,1);

% Here we also define distance-from-boundary constraints.
g = ga(@(g)vectorized_fitness(g,cloudOfPoints,[L;S]), nVariables,...
       [],[], [],[], [L 0],[S min(S-L)], [], options);
x = g(1); y = g(2); r = g(3);
%{
plot(x,y,'ro'); plot(x,y,'r*'); 
rectangle('Position',[x-r,y-r,2*r,2*r],'Curvature',[1 1],'EdgeColor','r'); 
%}

function f = vectorized_fitness(genes,pts,extent)
% genes = [x,y,r]
% extent = [Xmin Ymin; Xmax Ymax]
% f, the fitness, is the largest radius.
f = min(pdist2(genes(:,1:2), pts, 'euclidean'), [], 2);
% Instant death if circle contains a point:
f( f < genes(:,3) ) = Inf;
% Instant death if circle is too close to boundary:
f( any( genes(:,3) > genes(:,1:2) - extent(1,:) | ...
        genes(:,3) > extent(2,:) - genes(:,1:2), 2) ) = Inf;
% Note: this condition may possibly be specified using the A,b inputs of ga().
f(isfinite(f)) = -genes(isfinite(f),3);
%DEBUG: 
%{
     scatter(genes(:,1),genes(:,2),10 ,[0, .447, .741] ,'o'); % All
 z = ~isfinite(f); scatter(genes(z,1),genes(z,2),30,'r','x'); % Killed
 z =  isfinite(f); scatter(genes(z,1),genes(z,2),30,'g','h'); % Surviving
 [~,I] = sort(f); scatter(genes(I(1:5),1),genes(I(1:5),2),30,'y','p'); % Elite
%}

Y aquí hay un diagrama de "lapso de tiempo" de 47 generaciones de una carrera típica:

(Donde los puntos azules son la generación actual, las cruces rojas son organismos "muertos instantáneamente", los hexagramas verdes son los organismos "no muertos instantáneamente" y el círculo rojo marca el destino).


No estoy acostumbrado al procesamiento de imágenes, así que es solo una idea:

Implemente algo como un filtro gaussiano (desenfoque) que transforma cada partícula (píxeles) en un gradiente redondo con r = image_size (todos superpuestos). De esta manera, debería obtener una imagen donde la mayoría de los píxeles blancos deberían ser los mejores resultados. Desafortunadamente, la demostración en gimp falló porque el desenfoque extremo hizo que los puntos desaparecieran.

Alternativamente, puede extender incrementalmente todos los píxeles existentes marcando todos los píxeles vecinos en un área (ejemplo: r = 4), los píxeles restantes serían el mismo resultado (aquellos con la mayor distancia a cualquier píxel)


Puede usar bwdist de Image Processing Toolbox para calcular la transformación de distancia de la imagen. Esto puede considerarse como un método para crear un diagrama de voronoi que bien explicado en la respuesta de @ AnderBiguri.

img = imread('AbmxL.jpg');
%convert the image to a binary image
points = img(:,:,3)<200;
%compute the distance transform of the binary image
dist = bwdist(points);
%find the circle that has maximum radius
radius = max(dist(:));
%find position of the circle
[x y] = find(dist == radius);
imshow(dist,[]);
hold on
plot(y,x,'ro');







particles