Área do Lago do Parque do Ingá (Maringá) - Monte Carlo

Como podemos usar um monte de pequenas pedras para calcular a área de um lago?

Primeiramente devemos ter uma área A conhecida onde o lago está inserido (A = 300m x 500m = 150000 m*m):

Agora, jogamos N pedras aleatoriamente de tal forma que estas pedras caiam dentro da área A. Então, temos um número n de pedras que caem dentro do lago.

Assim, temos uma relação Área do lago = A*n/N.

Após jogarmos mil pedras (N = 1000), notamos que 354 pedras cairam dentro do lago ( n = 354):

Utilizando a relação, temos que a área do lago é: Área do lago = 150000*354/1000 = 53100 m*m.

O método estatístico utilizado para calcular a área do lago, é conhecido como Método de Monte Carlo. Quanto maior o número de pedras lançadas, melhor será a precisão da sua resposta. Para N muito grande, temos que a área do lago do parque do Ingá é igual à 55322 m*m.

A imagem do lago foi retirado do Google Map. O programa desenvolvido para simular as pedras foi escrito na linguagem C utilizando a biblioteca SDL para trabalhar com a imagem.

Algoritmos de Preenchimento

Area Open

No área open, é aplicado uma abertura somente às regiões planas ontem tem a área maior que um valor determinado.

Imagem original:

Image aplicada abertura somente para áreas maiores que 250px:

Algoritmo:

a = imread('ao.gif');

e = imread('e.png');

tomMax = max(max(a));

novo = a-a;

parametro = 250;

for j=1:tomMax,

imagem = logical(a==j);

[label, tamanho] = bwlabel(imagem,8);

novo2 = imagem-imagem;

for i = 1:tamanho,

aux = logical(label == i);

soma = sum(sum(aux));

if (soma >parametro),

novo2 = novo2+aux;

end

end

novo2 = uint8(novo2>0);

novo2 = imopen(novo2,e);

novo = max(novo, (novo2*j));

end;

Inf-Reconstrução

A inf-reconstrução, consiste em reconstruir um blob a partir de um marcador.

Aplicamos várias iterações das operações: intersecção da imagem original com o marcador dilatado.

O termo inf é utilizado porque paramos as iterações quando a nova imagem gerada não é alterada em comparação com a imagem da iteração anterior.

Imagem original:

Marcador:

Blog reconstruído:

Algoritmo:

A = imread('a.gif');

B = imread('b.gif');

E = imread('e.png');

Bnovo = min(imdilate(B,E),A);

while sum(sum(Bnovo))!=sum(sum(B))

___B = Bnovo;

___Bnovo = min(imdilate(B,E),A);

end;

Granulometría

Para calcular a granulometría, aumentamos o elemento estruturante em uma unidade a cada iteração e aplicamos uma operação de abertura à imagem original. Paramos as iterações quando for elimitados todos os pixeis brancos da imagem. Ou seja, a soma de pixeis brancos são iguais a zero.

Imagem a ser aplicado granulometria:

O elemento estruturante começa com uma matriz E = [1] de dimensões (1x1). A cada iteração é aumentado a matriz em uma dimensão na linha e na coluna.

O gráfico abaixo apresenta quantos pixeis brancos foram removidos, a cada iteração, em relação ao número total de pixeis brancos da imagem original:

O gráfico abaixo apresenta quantos pixeis brancos foram eliminados em relação a quantos pixeis brancos tinham na iteração anterior.

Algoritmo:

d = 1;

E = ones(d,d);

A = imread('barbara.png');

elim1 = [];

elim2 = [];

total = num = sum(sum(A));

while num!=0

____A = imopen(A,E);

____num_novo = sum(sum(A));

____elim1 = [elim1, total-num_novo];

____elim2 = [elim2, num - num_novo];

____num = num_novo;

____d = d + 1;

____E = ones(d,d);

end;

Fechamento - Propriedades

Algoritmo:

A = imread('a.png');

A == imclose(A,[1]);

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

B = imread('b.png');

E = imread('e.png');

AfE = imclose(A,E);

BfE = imclose(B,E);

imshow(AfE);

imshow(BfE);

Algoritmo:

A = imread('a.png');

E = imread('e.png');

AfE = imclose(A,E);

imshow(AfE);

Algoritmo:

A = imread('a.png');

E = imread('e.png');

AfE = imclose(A,E);

AfEfE = imclose(AfE,E);

AfE == AfEfE

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

E = imread('e.png);

AfE = imclose(A,E);

AfE == not(imopen(not(A),E'));

Resposta:

ans = 1;

Abertura - Propriedades

Algoritmo:

A = imread('a.png');

imopen(A,[1]) == A;

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

B = imread('b.png');

E = imread('e.png');

AoE = imopen(A,E);

BoE = imopen(B,E);

imshow(AoE);

imshow(BoE);

Algoritmo:

A = imread('a.png');

E = imread('e.png');

AoE = imopen(A,E);

imshow(AoE);

imshow(A);

Algoritmo:

A = imread('a.png');

E = imread('e.png');

AoE = imopen(A,E);

AoEoE = imopen(AoE,E);

AoE == AoEoE;

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

E = imread('b.png');

AoE = imopen(A,E);

AoE == not(imclose(not(A),E'));

Resposta:

ans = 1;

Erosão - Propriedades

Algoritmo:

A = imread('a.png');

E = imread('e.png');

AeE = imerode(A,E);

EeA = imerode(E,A);

AeE != EeA

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

E = imread('e.png');

AeE = imerode(A,E);

imshow(AeE);

imshow(A);

Algoritmo:

A = imread('a.png');

E1 = imread('e1.png');

E2 = imread('e2.png');

AeE1 = imerode(A,E1);

AeE2 = imerode(A,E2);

imshow(AeE1);

imshow(AeE2);

Algoritmo:

A = imread('a.png');

B = imread('b.png');

AiB = min(A,B);

AiBeE = imerode(AiB,E);

AeE = imerode(A,E);

BeE = imerode(B,E);

AeEiBeE = min(AeE,BeE);

AeEiBeE == AiBeE

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

E = imread('b.png');

AeE = imerode(A,E);

dual == not(imdilate(not(A),E'));

Resposta:

ans = 1;

Dilatação - Propriedades

Algoritmo:

A = imread('a.png');

B = imread('b.png);

imdilate(A,B) == imdilate(B,A);

Resposta:

asc = 1;

Algoritmo:

A = imread('a.png');

E1 = imread('e1.png');

E2 = imread('e2.png');

AdE1 = imdilate(A,E1);

AdE2 = imdilate(A,E2);

imshow(AdE1);

imshow(AdE2);

Algoritmo:

A = imread('a.png');

B = imread('b.png');

E = imread('e.png);

AdE = imdilate(A,E);

BdE = imdilate(B,E);

AdEuniaoBdE = max(AdE,BdE);

AuniaoB = max(A,B);

AuniaoBdE = imdilate(AuniaoB,E);

AdEuniaoBdE == AuniaoBdE

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

E = imread('e.png');

AdE = imdilate(A,E);

imshow(A);

imshow(AdE);

Algoritmo:

A = imread('a.png');

B = imread('b.png');

C = imread('c.png');

AdB = imdilate(A,B);

AdB_dC = imdilate(AdB,C);

BdC = imdilate(B,C);

Ad_BdC = imdilate(A,BdC);

AdB_dC == Ad_BdC;

Resposta:

ans = 1;

Algoritmo:

A = imread('a.png');

B = imread('b.png');

AdB = imdilate(A,B);

AdB == not(imerode(not(A),B'));

Resposta:

ans = 1;

União, Intersecção e Complemento com Imagens

Conjuntos

Seja A e B conjuntos:

• União (U): A U B = {x є E: x є A ou x є B}
• Intersecção(∩): A ∩ B : {x є E: x є A e x є B}
• Complemento ('): A' = {x є E : x ɇ A}

Imagens

Imagem A:

Imagem B:

A união B:

A intersecção B:

Complemento de (A intersecção B):

Algoritmo

%LENDO IMAGENS

A = imread('a.png');

B = imread('b.png');

%UNIÃO ENTRE A IMAGEM A E B.

U = max(A,B);

imwrite(U, 'u.png');

%INTERSECÇÃO ENTRE A IMAGEM A E B

I = min(A,B);

imwrite(I, 'i.png');

%COMPLEMENTO DA IMAGEM INTERSECÇÃO

N = not(I);

imwrite(N,'n.png');

vídeos antigos

Modo curioso para o cálculo de π

Um modo curioso para o cálculo de π é o seguinte:

Desenhe em uma folha de papel uma série de retas paralelas, de modo que a distância entre duas retas consecutivas seja de 5 cm. Corte dez palitos de 2.5 cm. De uma altura de 30 cm solte os 10 palitos de uma vez sobre a folha e conte quantos deles tocam alguma das retas desenhadas. Repita várias vezes. Dividindo o número total de palitos lançados pelo número de palitos que cruzam alguma das retas, obtém­se um número próximo de π. Quanto maior número de lançamentos, mais próximo de π será a resposta.

Este método é conhecido como Agulha de Buffon, ele foi elaborado em 1777, e envolvia o lançamento aleatório de uma agulha em um plano com infinitas linhas paralelas e a determinação da probabilidade de que a agulha cruze uma das linhas. Dizem que Buffon descobriu esse método em um período em que esteve com a perna quebrada, entediado, lançava o palito do seu cachimbo no chão de tábuas do seu quarto. Então, descobriu que a probabilidade dos lançamentos cruzarem as linhas está intimamente relacionada com o valor de π.

Mas por que funciona???

O caso mais simples, é quando o espaço entre duas linhas tem a distância de 1 unidade e a agulha também tem o tamanho de 1 unidade.

Neste problema tem mais duas variáveis: o ângulo (Ө) em que a agulha cai em relação as linhas e a distância (D) do centro da agulha até a linha mais próxima. Como pode ser verificado na imagem abaixo:

Podemos notar que o ângulo Ө pertence ao intervalo [0,π] e que d = sen(Ө)/2. O d é o cateto oposto (CO) e a hipotenusa (H) vale 1/2 unidade, então, sen(Ө) = CO/H .

Dados essas condições, devemos tentar predizer, baseado somente em Ө, quando que a agulha cruza uma das linhas.

Nota-se que a agulha cruza uma das linhas quando D < d, ou seja, D < sen(Ө)/2. Agora, basta analizar quando D < sen(Ө)/2. Vamos analizar o gráfico de sen(Ө)/2 no intervalo [0,π]:

Podemos notar que o valor máximo da função, coincide com o valor máximo de D, que é igual à 1/2 unidade. Ou seja, conseguimos tirar duas informações sobre um ponto aleatório dentro desse gráfico: D e d. Assim, quando este ponto esta abaixo da linha sen(Ө)/2, D < d, então, a agulha cruza a linha. De modo contrário, quando o ponto esta acima da linha sen(Ө)/2, D > d, então, a agulha não cruza a linha.

Logo, a probabilidade (p) de que a agulha cruze uma linha será dada pela razão entre a área do retângulo (0.5 x π) e a área abaixo da função d = sen(Ө)/2.

Temos que as áreas são:

Aretangulo = π/2

Ad = 1 (basta calcular a integral de d no intervalo [0, π]).

Então, temos:

probabilidade = Ad / Aretangulo = 1 / ( π/2) = 2/π.

probabilidade = Número de sucessos / Número de Eventos. Ou seja, número de agulhas que cruzam as linhas (NC) / número total de agulhas lançadas (TA).

Igualando as probabilidade's temos:

NC / TA = 2/π

Enfim, isolando π:

π = (número total de agulhas lançadas) * 2 / (número de agulhas que cruzam as linhas)