matlab math image-processing geometry particle particles

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



math image-processing (5)

¡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.

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.

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).

Me gustaría proponer otra solución basada en una búsqueda de cuadrícula con refinamiento. No es tan avanzado como el de Ander o tan corto como el de rahnema1, pero debería ser muy fácil de seguir y comprender. Además, funciona bastante rápido.

El algoritmo contiene varias etapas:

  1. Generamos una cuadrícula uniformemente espaciada.
  2. Encontramos las distancias mínimas de puntos en la cuadrícula a todos los puntos proporcionados.
  3. Descartamos todos los puntos cuyas distancias están por debajo de un determinado percentil (por ejemplo, 95º).
  4. Elegimos la región que contiene la mayor distancia (esta debería contener el centro correcto si mi grilla inicial es lo suficientemente fina).
  5. Creamos una nueva malla de malla alrededor de la región elegida y encontramos distancias nuevamente (esta parte es claramente subóptima, porque las distancias se calculan en todos los puntos, incluidos los lejanos y los irrelevantes).
  6. Repetimos el refinamiento dentro de la región, mientras observamos la varianza del 5% superior de los valores -> si cae por debajo de un umbral predeterminado, lo rompemos.

Varias notas:

  • He asumido que los círculos no pueden ir más allá de la extensión de los puntos dispersos (es decir, el cuadrado delimitador de la dispersión actúa como un "muro invisible").
  • El percentil apropiado depende de qué tan fina sea la cuadrícula inicial. Esto también afectará la cantidad de iteraciones while y el valor inicial óptimo para cnt .

function [xBest,yBest,R] = q42806059 rng(1) x=rand(1,100)*5; y=rand(1,100)*5; %% Find the approximate region(s) where there exists a point farthest from all the rest: xExtent = linspace(min(x),max(x),numel(x)); yExtent = linspace(min(y),max(y),numel(y)).''; % Create a grid: [XX,YY] = meshgrid(xExtent,yExtent); % Compute pairwise distance from grid points to free points: D = reshape(min(pdist2([XX(:),YY(:)],[x(:),y(:)]),[],2),size(XX)); % Intermediate plot: % figure(); plot(x,y,''.k''); hold on; contour(XX,YY,D); axis square; grid on; % Remove irrelevant candidates: D(D<prctile(D(:),95)) = NaN; D(D > xExtent | D > yExtent | D > yExtent(end)-yExtent | D > xExtent(end)-xExtent) = NaN; %% Keep only the region with the largest distance L = bwlabel(~isnan(D)); [~,I] = max(table2array(regionprops(''table'',L,D,''MaxIntensity''))); D(L~=I) = NaN; % surf(XX,YY,D,''EdgeColor'',''interp'',''FaceColor'',''interp''); %% Iterate until sufficient precision: xExtent = xExtent(~isnan(min(D,[],1,''omitnan''))); yExtent = yExtent(~isnan(min(D,[],2,''omitnan''))); cnt = 1; % increase or decrease according to the nature of the problem while true % Same ideas as above, so no explanations: xExtent = linspace(xExtent(1),xExtent(end),20); yExtent = linspace(yExtent(1),yExtent(end),20).''; [XX,YY] = meshgrid(xExtent,yExtent); D = reshape(min(pdist2([XX(:),YY(:)],[x(:),y(:)]),[],2),size(XX)); D(D<prctile(D(:),95)) = NaN; I = find(D == max(D(:))); xBest = XX(I); yBest = YY(I); if nanvar(D(:)) < 1E-10 || cnt == 10 R = D(I); break end xExtent = (1+[-1 +1]*10^-cnt)*xBest; yExtent = (1+[-1 +1]*10^-cnt)*yBest; cnt = cnt+1; end % Finally: % rectangle(''Position'',[xBest-R,yBest-R,2*R,2*R],''Curvature'',[1 1],''EdgeColor'',''r'');

El resultado que obtengo para los datos de ejemplo de Ander es [x,y,r] = [0.7832, 2.0694, 0.7815] (que es lo mismo). El tiempo de ejecución es aproximadamente la mitad de la solución de Ander.

Aquí están las parcelas intermedias:

Contorno de la mayor distancia (clara) desde un punto al conjunto de todos los puntos proporcionados:

Después de considerar la distancia desde el límite, manteniendo solo el 5% superior de los puntos distantes, y considerando solo la región que contiene la mayor distancia (el pedazo de superficie representa los valores guardados):

Y finalmente:

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'');