Ajustement du plus grand cercle dans la zone libre de l'image avec une particule distribuée

Question

Je travaille sur des images permettant de détecter et d'ajuster le plus grand cercle possible dans l'une des zones libres d'une image contenant des particules distribuées :

(capable de détecter l'emplacement de la particule).

L'une des solutions consiste à définir un cercle touchant toute combinaison de 3 points, à vérifier si le cercle est vide, puis à trouver le plus grand cercle parmi tous les cercles vides. Cependant, cela conduit à un nombre énorme de combinaisons, à savoir C(n,3) où n est le nombre total de particules dans l'image.

J'apprécierais si quelqu'un pouvait me donner un indice ou une méthode alternative que je pourrais explorer.

Answer 1

Faisons un peu de maths mon ami, car les maths mènent toujours à la fin !

Wikipedia :

En mathématiques, un diagramme de Voronoï est un partitionnement d'un plan en régions basées sur la distance aux points dans un sous-ensemble spécifique du plan.

Par exemple :

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

voronoi(x,y);

Ce qui est bien dans ce diagramme, c'est que si vous remarquez, tous les bords/vertices de ces zones bleues sont tous à égale distance des points qui les entourent. Ainsi, si nous connaissons l'emplacement des sommets et que nous calculons les distances aux points les plus proches, nous pouvons choisir le sommet présentant la plus grande distance comme centre du cercle.

Il est intéressant de noter que les bords d'une région de Voronoï sont également définis comme les circoncentres des triangles générés par une triangulation de Delaunay.

Donc si nous calculons la triangulation de Delaunay de la zone, et leurs circoncentres

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

Et calculez les distances entre les circonférences et l'un des points qui définissent chaque triangle :

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

Nous avons alors le centre ( cc ) et le rayon ( distance ) de tous les cercles possibles qui n'ont pas de point à l'intérieur. Il nous faut juste le plus grand !

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

Maintenant, nous allons tracer

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

Remarquez que le centre du cercle se trouve sur un sommet du diagramme de Voronoï.

---

NOTE Cette méthode trouvera le centre dans [0-5], [0-5]. Vous pouvez facilement la modifier pour changer cette contrainte. Vous pouvez également essayer de trouver le cercle qui s'inscrit dans sa totalité dans la zone intéressée (et non pas seulement le centre). Cela nécessiterait un petit ajout à la fin où le maximum est obtenu.

Answer 2

J'aimerais proposer une autre solution basée sur une recherche en grille avec raffinement. Elle n'est pas aussi avancée que celle d'Ander ni aussi courte que celle de rahnema1, mais elle devrait être très facile à suivre et à comprendre. En outre, il fonctionne assez rapidement.

L'algorithme comporte plusieurs étapes :

1. Nous générons une grille à espacement régulier.
2. Nous trouvons les distances minimales des points de la grille à tous les points fournis.
3. Nous éliminons tous les points dont la distance est inférieure à un certain percentile (par exemple 95e).
4. Nous choisissons la région qui contient la plus grande distance (celle-ci devrait contenir le centre correct si ma grille initiale est suffisamment fine).
5. Nous créons une nouvelle grille de maillage autour de la région choisie et trouvons à nouveau les distances (cette partie est clairement sous-optimale, car les distances sont calculées pour tous les points, y compris les points éloignés et non pertinents).
6. Nous itérons le raffinement à l'intérieur de la région, tout en gardant un œil sur la variance des 5 % de valeurs les plus élevées -> si elle tombe en dessous d'un seuil prédéfini, nous rompons.

Plusieurs notes :

• Je suis parti du principe que les cercles ne peuvent pas dépasser l'étendue des points dispersés (c'est-à-dire que le carré de délimitation de la dispersion agit comme un "mur invisible").

• Le percentile approprié dépend de la finesse de la grille initiale. Cela affectera également la quantité de while itérations, et la valeur initiale optimale de 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,2R,2R],'Curvature',[1 1],'EdgeColor','r');

Le résultat que j'obtiens pour les données de l'exemple d'Ander est le suivant [x,y,r] = [0.7832, 2.0694, 0.7815] (ce qui revient au même). Le temps d'exécution est environ la moitié de celui de la solution d'Ander.

Voici les tracés intermédiaires :

Contour de la plus grande distance (nette) entre un point et l'ensemble de tous les points fournis :

Après avoir considéré la distance par rapport à la frontière, en ne gardant que les 5% de points les plus éloignés, et en ne considérant que la région qui contient la plus grande distance (le morceau de surface représente les valeurs gardées) :

Et enfin :

Answer 3

Vous pouvez utiliser bwdist de Image Processing Toolbox pour calculer la transformée de distance de l'image. Ceci peut être considéré comme une méthode pour créer un diagramme de Voronoï qui est bien expliqué dans la réponse 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');

Answer 4

Le fait que ce problème puisse être résolu à l'aide d'une "recherche directe" (comme on peut le voir dans l'exemple suivant) est une bonne chose. une autre réponse ) signifie que l'on peut considérer cela comme un optimisation globale problème. Il existe plusieurs façons de résoudre de tels problèmes, chacune étant adaptée à certains scénarios. Par curiosité personnelle, j'ai décidé de résoudre ce problème à l'aide d'une algorithme génétique .

D'une manière générale, un tel algorithme nous oblige à considérer la solution comme un ensemble de "gènes" soumis à une "évolution" sous une certaine "fonction d'adéquation". Il se trouve qu'il est assez facile d'identifier les gènes et la fonction d'aptitude dans ce problème :

• Gènes : x , y , r .
• Fonction d'aptitude : techniquement, il s'agit de l'aire maximale du cercle, mais cela équivaut à la valeur maximale de la fonction d'aptitude. r (ou minimum -r puisque l'algorithme requiert une fonction pour minimiser ).
• Contrainte spéciale - si r est plus grande que la distance euclidienne au plus proche des points fournis (c'est-à-dire que le cercle contient un point), l'organisme "meurt".

Vous trouverez ci-dessous une implémentation de base d'un tel algorithme (" base "parce qu'il n'est pas du tout optimisé et qu'il y a beaucoup de place pour l'optimisation. sans mauvais jeu de mots dans ce problème).

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
%}

Et voici un tracé "time-lapse" de 47 générations d'un parcours typique :

(Les points bleus représentent la génération actuelle, les croix rouges les organismes "insta-killed", les hexagrammes verts les organismes "non-insta-killed" et le cercle rouge la destination).

Answer 5

Je ne suis pas habitué au traitement de l'image, donc c'est juste une idée :

Implémentez quelque chose comme un filtre gaussien (flou) qui transforme chaque particule (pixels) en un gradient rond avec r=taille_de_l'image (tous se chevauchant). De cette façon, vous devriez obtenir une image où les pixels les plus blancs devraient donner les meilleurs résultats. Malheureusement, la démonstration dans gimp a échoué car le flou extrême a fait disparaître les points.

Vous pourriez aussi étendre de manière incrémentielle tous les pixels existants en marquant tous les pixels voisins dans une zone (exemple : r=4), les pixels restants seraient le même résultat (ceux avec la plus grande distance à tout pixel)