¿Cuál es una forma simple de calcular la superposición entre una imagen y un polígono?

Question

Tengo un polígono cerrado que no se intersecta a sí mismo. Sus vértices se guardan en dos vectores X e Y. Finalmente, los valores de X e Y están limitados entre 0 y 22.

Me gustaría construir una matriz de tamaño 22x22 y establecer el valor de cada contenedor igual a true si parte del polígono se superpone con ese bin, de lo contrario es falso.

Mi idea inicial era generar una grilla de puntos definidos con [a, b] = meshgrid(1:22) y luego usar inpolygon para determinar qué puntos de la grilla estaban en el polígono.

[a b] = meshgrid(1:22); 
inPoly1 = inpolygon(a,b,X,Y); 

Sin embargo, esto sólo devuelve verdadero si si el centro de la bandeja está contenida en el polígono, es decir, que devuelve la forma rojo en la imagen a continuación. Sin embargo, qué necesidad es más a lo largo de las líneas de la forma verde (aunque sigue siendo una solución incompleta).

Para obtener el blob verde realicé cuatro llamadas al inpolygon. Para cada comparación, cambié la cuadrícula de puntos NE, NW, SE o SW por 1/2. Esto es equivalente a probar si las esquinas de un contenedor están en el polígono.

inPoly2 = inpolygon(a-.5,b-.5,X,Y) | inpolygon(a+.5,b-.5,X,Y) | inpolygon(a-.5,b+5,X,Y) | inpolygon(a+.5,b+.5,X,Y); 

Si bien esto me proporciona con una solución parcial se produce un error en el caso en que un vértice es contener en un contenedor pero ninguna de las esquinas de basura son.

¿Hay una manera más directa de atacar este problema, preferiblemente con una solución que produce un código más legible?

Esta parcela fue dibujado con:

imagesc(inPoly1 + inPoly2); hold on; 
line(a, b, 'w.'); 
line(X, Y, 'y); 

Answer 1

Una sugerencia es usar la función polybool (no disponible en 2008b o anterior). Encuentra la intersección de dos polígonos y devuelve los vértices resultantes (o un vector vacío si no existen vértices). Para usarlo aquí, iteramos (usando arrayfun) sobre todos los cuadrados de su cuadrícula para ver si el argumento de salida a polybool está vacío (por ejemplo, no hay superposición).

N=22; 
sqX = repmat([1:N]',1,N); 
sqX = sqX(:); 
sqY = repmat(1:N,N,1); 
sqY = sqY(:); 

intersects = arrayfun((@(xs,ys) ... 
     (~isempty(polybool('intersection',X,Y,[xs-1 xs-1 xs xs],[ys-1 ys ys ys-1])))),... 
     sqX,sqY); 

intersects = reshape(intersects,22,22); 

Aquí es la imagen resultante:

Código para el trazado:

imagesc(.5:1:N-.5,.5:1:N-.5,intersects'); 
hold on; 
plot(X,Y,'w'); 
for x = 1:N 
    plot([0 N],[x x],'-k'); 
    plot([x x],[0 N],'-k'); 
end 
hold off; 

Answer 2

¿Qué hay de este algoritmo pseudocódigo:

For each pair of points p1=p(i), p2=p(i+1), i = 1..n-1 
    Find the line passing through p1 and p2 
    Find every tile this line intersects // See note 
    Add intersecting tiles to the list of contained tiles 

Find the red area using the centers of each tile, and add these to the list of contained tiles 

Nota: Esta línea tendrá un poco de esfuerzo para poner en práctica, pero Creo que hay un algoritmo bastante sencillo y conocido para eso.

Además, si estaba usando .NET, simplemente definiría un rectángulo correspondiente a cada mosaico de cuadrícula, y luego ver cuáles se cruzan con el polígono. Sin embargo, no sé si controlar la intersección es fácil en Matlab.

Answer 3

Primero definir un círculo de baja resolución para este ejemplo

X=11+cos(linspace(0,2*pi,10))*5; 
Y=11+sin(linspace(0,2.01*pi,10))*5; 

Al igual que su ejemplo se ajusta con una cuadrícula de ~ 22 unidades. Luego, siguiendo su ejemplo, declaramos una malla de malla y verificamos si hay puntos en el polígono.

stepSize=0.1; 
[a b] = meshgrid(1:stepSize:22); 
inPoly1 = inpolygon(a,b,X,Y); 

La única diferencia es que cuando su solución original tomó medidas de una, esta cuadrícula puede tomar pasos más pequeños. Y, por último, para incluir cualquier cosa dentro de los "bordes" de las plazas

inPolyFull=unique(round([a(inPoly1) b(inPoly1)]) ,'rows'); 

El round simplemente toma nuestra red de alta resolución y rondas de los puntos adecuadamente a sus equivalentes de baja resolución más cercanas. A continuación, eliminamos todos los duplicados en un estilo vectorial o en forma de pares llamando al unique con el calificador 'rows'.Y eso es todo

Para ver el resultado,

[aOrig bOrig] = meshgrid(1:22); 
imagesc(1:stepSize:22,1:stepSize:22,inPoly1); hold on; 
plot(X,Y,'y'); 
plot(aOrig,bOrig,'k.'); 
plot(inPolyFull(:,1),inPolyFull(:,2),'w.'); hold off; 

Cambio de la stepSize tiene el efecto esperado de mejorar el resultado a costa de la velocidad y la memoria.

Si necesita que el resultado sea en el mismo formato que el inPoly2 en su ejemplo, puede utilizar

inPoly2=zeros(22); 
inPoly2(inPolyFull(:,1),inPolyFull(:,2))=1 

Espero que ayude. Puedo pensar en otras formas de hacerlo, pero esto parece ser el más directo.

Answer 4

se recomienda usar poly2mask en el procesamiento de imágenes caja de herramientas, hace más o menos lo usted quiere, creo, y también más o menos lo que usted y Salain sugirieron.

Answer 5

Ligera mejora

En primer lugar, para simplificar su "solución parcial" - lo que estás haciendo es simplemente buscando en las esquinas. Si en lugar de considerar la cuadrícula 22x22 de puntos, podría considerar la cuadrícula 23x23 de esquinas (que se compensará de la cuadrícula más pequeña por (-0.5, -0.5). Una vez que tenga eso, puede marcar los puntos en la grilla 22x22 que tienen al menos una esquina en el polígono

solución completa:..

sin embargo, lo que realmente está buscando es si el polígono se cruza con la caja de 1x1 que rodea cada píxel Esto no significa necesariamente incluya cualquiera de las esquinas, pero requiere requiere que el polígono intersecte uno de los cuatro lados de la caja.

Una forma usted podría encontrar los píxeles en el que el polígono se cruza con la caja que contiene es con el siguiente algoritmo:

For each pair of adjacent points in the polygon, calling them pA and pB: 
    Calculate rounded Y-values: Round(pA.y) and Round(pB.y) 
    For each horizontal pixel edge between these two values: 
     * Solve the simple linear equation to find out at what X-coordinate 
      the line between pA and pB crosses this edge 
     * Round the X-coordinate 
     * Use the rounded X-coordinate to mark the pixels above and below 
      where it crosses the edge 
    Do a similar thing for the other axis 

Así, por ejemplo, decimos que estamos mirando pA = (1, 1) y pB = (2, 3).

• Primero, calculamos los valores Y redondeados: 1 y 3.
• A continuación, nos fijamos en los bordes de los píxeles entre estos valores: y = 1.5 y y = 2.5 (bordes de los píxeles están a mitad de desplazamiento desde píxeles
• Para cada uno de éstos, se resuelve la ecuación lineal para encontrar donde pA ->pB se cruza con nuestra . bordes Esto nos da: x = 1.25, y = 1.5, y x = 1.75, y = 2.5
• para cada una de estas intersecciones, tomamos el X-valor redondeado, y lo utilizan para marcar los píxeles a ambos lados del borde
  • x = 1.25 se redondea a 1.. (para el borde y = 1.5). Por lo tanto, podemos marcar el píxeles en (1, 1) y (1, 2) como parte de nuestro conjunto.
  • x = 1.75 se redondea a 2 (para el borde y = 2.5). Por lo tanto, podemos marcar los píxeles en (2, 2) y (2, 3).

Así que por los bordes horizontales atendidos. A continuación, vamos a ver las verticales:

• En primer lugar, calcular los X-valores redondeados: 1 y 2
• A continuación, nos fijamos en los bordes de los píxeles. Aquí, solo hay uno: x = 1.5.
• Para este borde, encontramos el que coincide con la línea pA ->pB. Esto nos da x = 1.5, y = 2.
• Por esta intersección, tomamos la Y-valor redondeado, y lo utilizan para marcar píxeles cada lado de la orilla:
  • y = 2 se redondea a 2. Por tanto, podemos marcar los píxeles a (1, 2) y (2, 2).

hecho!

Bueno, más o menos. Esto le dará los bordes, pero no completará el cuerpo del polígono. Sin embargo, puede combinarlos con sus resultados previos (rojos) para obtener el conjunto completo.

Answer 6

Bueno, creo que llego tarde, aunque estrictamente hablando, el tiempo de la recompensa fue hasta mañana;). Pero aquí va mi intento. Primero, una función que marca celdas que contienen/tocan un punto. Dada una rejilla estructurado con lx espaciamiento, Ly, y un conjunto de puntos con coordenadas (xp, yp), establecer células que contienen:

function cells = mark_cells(lx, ly, Xp, Yp, cells) 

% Find cell numbers to which points belong. 
% Search by subtracting point coordinates from 
% grid coordinates and observing the sign of the result. 
% Points lying on edges/grid points are assumed 
% to belong to all surrounding cells. 

sx=sign(bsxfun(@minus, lx, Xp')); 
sy=sign(bsxfun(@minus, ly, Yp')); 
cx=diff(sx, 1, 2); 
cy=diff(sy, 1, 2); 

% for every point, mark the surrounding cells 
for i=1:size(cy, 1) 
    cells(find(cx(i,:)), find(cy(i,:)))=1; 
end 
end 

Ahora, el resto del código. Para cada segmento en el polígono (debe recorrer los segmentos uno por uno), interseque el segmento con las líneas de la cuadrícula. La intersección se realiza con cuidado, para líneas horizontales y verticales por separado, utilizando las coordenadas del punto de la cuadrícula para evitar imprecisiones numéricas. Para los puntos de intersección encontrado que llamo mark_cells para marcar las células que rodean a 1:

% example grid 
nx=21; 
ny=51; 
lx = linspace(0, 1, nx); 
ly = linspace(0, 1, ny); 
dx=1/(nx-1); 
dy=1/(ny-1); 
cells = zeros(nx-1, ny-1); 

% for every line in the polygon... 
% Xp and Yp contain start-end points of a single segment 
Xp = [0.15 0.61]; 
Yp = [0.1 0.78]; 

% line equation 
slope = diff(Yp)/diff(Xp); 
inter = Yp(1) - (slope*Xp(1)); 

if isinf(slope) 
    % SPECIAL CASE: vertical polygon segments 
    % intersect horizontal grid lines 
    ymax = 1+floor(max(Yp)/dy); 
    ymin = 1+ceil(min(Yp)/dy); 
    x=repmat(Xp(1), 1, ymax-ymin+1); 
    y=ly(ymin:ymax); 
    cells = mark_cells(lx, ly, x, y, cells); 
else 
    % SPECIAL CASE: not horizontal polygon segments 
    if slope ~= 0 
     % intersect horizontal grid lines 
     ymax = 1+floor(max(Yp)/dy); 
     ymin = 1+ceil(min(Yp)/dy); 
     xmax = (ly(ymax)-inter)/slope; 
     xmin = (ly(ymin)-inter)/slope; 
     % interpolate in x... 
     x=linspace(xmin, xmax, ymax-ymin+1); 
     % use exact grid point y-coordinates! 
     y=ly(ymin:ymax); 
     cells = mark_cells(lx, ly, x, y, cells); 
    end 

    % intersect vertical grid lines 
    xmax = 1+floor(max(Xp)/dx); 
    xmin = 1+ceil(min(Xp)/dx); 
    % interpolate in y... 
    ymax = inter+slope*lx(xmax); 
    ymin = inter+slope*lx(xmin); 
    % use exact grid point x-coordinates! 
    x=lx(xmin:xmax); 
    y=linspace(ymin, ymax, xmax-xmin+1); 
    cells = mark_cells(lx, ly, x, y, cells); 
end 

de salida para el ejemplo de malla/segmento:

Caminando a través de todos los segmentos del polígono le da el polígono 'halo'.Finalmente, el interior del polígono se obtiene utilizando la función estándar de polígono. Avíseme si necesita más detalles sobre el código.