Derivació Numèrica

Per aproximar numèricament la derivada d'una funció en un punt, s'utilitzen el que s'anomenen mètodes de diferències finites (que no són res més que la comparació amb valors propers).

Derivació en R

Es poden utilitzar diversos esquemes per calcular la derivada. Depenent del esquema escollit, la precissió serà millor.
% Considerem una funció f(x)=sin(x) de la qual la derivada, f'(x)=cos(x), és coneguda.
% Aproximem el seu valor numèricament per comparar la precisió del càlcul.
%
f = @(x) sin(x); %definició 'inline' de la funció
derf = @(x) cos(x); %definició de la derivada exacta
a = 1.2; %punt en el avaluem la derivada
dist = [1.e-2, 1.e-3, 1.e-4, 1.e-5]; %distancies a les que ens desplacem
for i = 1:length(dist)
h = dist(i);
derDavant = (f(a+h)-f(a))/h; %diferència endavant
error1 = abs(derDavant-derf(a)); %calculem l'error respecte la derivada exacta
derCentrada = (f(a+h)-f(a-h))/(2*h); %diferència centrada
error2 = abs(derCentrada-derf(a)); %calculem l'error respecte la derivada exacta
fprintf([' h= %.1e, derDavant= %.4e, er= %.2e, derCent= ' ...
'%.4e, er= %.2e \n'],h,derDavant,error1,derCentrada,error2); %escrivim els valors
end
h= 1.0e-02, derDavant= 3.5769e-01, er= 4.67e-03, derCent= 3.6235e-01, er= 6.04e-06 h= 1.0e-03, derDavant= 3.6189e-01, er= 4.66e-04, derCent= 3.6236e-01, er= 6.04e-08 h= 1.0e-04, derDavant= 3.6231e-01, er= 4.66e-05, derCent= 3.6236e-01, er= 6.04e-10 h= 1.0e-05, derDavant= 3.6235e-01, er= 4.66e-06, derCent= 3.6236e-01, er= 2.43e-12
%mostrem els resultats (fixeu-vos en la última columna!!)

Derivades Parcials

Si ara considerem una funció , els mateixos esquemes de diferències finites es poden aplicar a les derivades parcials, però fent els increments per cada component.
Equivalentment les seves versions de derivades endavant.
% Escollim ara la funció f(x,y)= sin(x)*cos(y)
fxy = @(x,y) sin(x)*cos(y);
parcX = @(x,y) cos(x)*cos(y); %derivada parcial resp. de x
parcY = @(x,y) -sin(x)*sin(y); %derivada parcial resp. de y
a=2;
b=3;
h=1.e-5;
derCentradaX=(fxy(a+h,b)-fxy(a-h,b))/(2*h); %aprox. derivada centrada resp. de x
derCentradaY=(fxy(a,b+h)-fxy(a,b-h))/(2*h); %aprox. derivada centrada resp. de y
grad=[derCentradaX,derCentradaY] %vector gradient en (a,b)
grad = 1×2
0.4120 -0.1283
errors=[abs(derCentradaX-parcX(a,b)),abs(derCentradaY-parcY(a,b))]
errors = 1×2
10-11 ×
0.2938 0.2899
%noteu que l'error és d'ordre h^2

Aplicació de la derivació: Detecció de Contorns en una Imatge

Imatge Digital

Una imatge digital no és rés més que una matriu els elements de la qual s'anomenen píxels. En cada element de la matriu es guarda el valor d'intensitat de llum de cada píxel (o tres en el cas de color RGB). D'aquesta forma, podem considerar una imatge com una funció de R que a cada pixel , li assigna un valor d'intensitat . Aquesta intensitat es mesura de a 255, si la intensitat s'expressa en enters, o bé de a 1 si és en valors reals. El mínim (= 0) correspon al negre i el màxim ( =1 o 255) correspon al blanc.
ImatgePixels.png
Fig 1. Idea de la informació d'una imatge digital
Com tota matriu, una imatge té files i columnes que representen els eixos i (vertical) i j (horitzontal) respectivament.

Derivada sobre una imatge

Considerant les coordenades d'un píxel i la seva funció intensitat , podem plantejar-nos calcular les seves derivades parcials a partir de les diferències finites. Notem que, en tractar-se de píxels, la distància entre ells és sempre
A les fotos amb contorns força ben delimitats hi ha zones amb intensitat semblant que es canvia significativament només quan passem d'una zona a l'altra. Podem fer servir això per detectar (extreure) els contorns: són conjunts de punts (píxels) on el gradient de la funció discreta aconsegueix valors grans.
Veiem un exemple sobre una imatge de monedes, coins.png, en que les derivades les farem per TOTS els pixels. Hem de tenir present que si la imatge és de color, la passarem primer a escala de grisos. Denotarem per A la matriu que conté la imatge per processar.
Imatge = imread('coins.png'); %carregar una imatge
[m,n,p] = size(Imatge); %les mides indiquen si es tracta d'una imatge de color (p=3 o de grisos p=1)
if ( p > 1) %en el cas d'una imatge RGB, la passem a grisos
Imatge = rgb2gray(Imatge);
end
A = im2double(Imatge); %passem els valors d'intensitat a doubles per poder fer càlculs
[m,n]=size(A);
 
figure()
imshow(Imatge); % mostrem la imatge
title('Imatge en grisos i doubles'); %posem títol

Contorns per files: Obtenim els horitzontals

Si calculem les derivades respecte la primera component (les files), la comparació es fa entre píxels que estan a la mateixa columna i per tant seria una comparació en 'vertical' que ens donarà els contorns horitzontals.
Per cada píxel , els elements que intervenen es poden veure a la figura:
vertical.png
contH=zeros(size(A)); %matriu de contorns horitzontals
%per la primera/última fila/columna no fem la derivada (es podria fer
%diferència endavant si fos necessari)
for i=2:m-1
for j=2:n-1
contH(i,j)=(A(i+1,j)-A(i-1,j))/2; %la distància entre pixels és h=1
end
end
figure()
imshow(contH)
title('Contorns Horitzontals');

Contorns per columnes: Obtenim els verticals

Anàlogament, si ho fem per columnes:
Per cada píxel , els elements que intervenen es poden veure a la figura:
horitzontal.png
contV=zeros(size(A)); %matriu de contorns en Y
for i=2:m-1
for j=2:n-1
contV(i,j)=(A(i,j+1)-A(i,j-1))/2; %la distància entre ixels és h=1
end
end
figure()
imshow(contV)
title('Contorns Verticals');

Contorns totals associats al gradient en cada pixel

Si tenim en compte els dos contorns anteriors, de fet és com tenir les derivades en y i en x respectivament, podem calcular per a cada píxel el gradient de la funció definit per , el seu mòdul del gradient serà el valor associat a cada píxel
gradient.png
% El podem calcular de dues formes: amb un doble bucle o directament
% aprofitant les dues matrius anteriors i fent la operació per elements que
% ofereix el Matlab
contA=zeros(size(A));
for i=2:m-1 %fem un doble bucle per recorre tots els píxels
for j=2:n-1
dAh=(A(i,j+1)-A(i,j-1))/2;
dAv=(A(i+1,j)-A(i-1,j))/2;
contA(i,j)=sqrt(dAh^2+dAv^2);
end
end
 
gradA=sqrt(contH.^2+contV.^2); %podem fer el càlcul directe per a tots els elements
max(gradA(:)-contA(:)) %veiem que el càlcul és el mateix
ans = 0
figure()
imshow(contA)
title('Imatge de gradient');

Mètode de Prewitt per millorar els contorns

Per intentar resaltar millor els contorns, Prewitt va proposar d'acumular les derivades tant en un pixel, com en els seus veïns a la dreta i l'esquerra (en vertical) i per sobre i sota (en horitzontal). Posteriorment es calcula el gradient resultant.
prewitt.png
contAPrew=zeros(size(A));
for i=2:m-1
for j=2:n-1
dAh1 = (A(i+1,j)-A(i-1,j))/2; %pels contorns horitzontals
dAh2 = (A(i+1,j+1)-A(i-1,j+1))/2;
dAh3 = (A(i+1,j-1)-A(i-1,j-1))/2;
dAh = dAh1 + dAh2 + dAh3 ;
 
dAv1 = (A(i,j+1)-A(i,j-1))/2; %pels contorns verticals
dAv2 = (A(i+1,j+1)-A(i+1,j-1))/2;
dAv3 = (A(i-1,j+1)-A(i-1,j-1))/2;
dAv = dAv1 + dAv2 + dAv3 ;
contAPrew(i,j)=sqrt(dAh^2+dAv^2); %gradient
end
end
figure()
imshow(contAPrew)
title('Filtre de Prewitt');

Exercici 1:

Feu el mateix estudi de contorns utilitzant la imatge BB8.jpg (contorns molt definits)

Exercici 2:

Apliqueu el càlcul de contorns a la imatge mandrill.jpg (contorns menys definits)
MandrillContorns.png
(c)Numerical Factory 2023