Derivació Numèrica

Per aproximar numèricament la derivada d'una funció en un punt, utilitzen el que s'anomenen diferències finites

Contents

Derivació en R

Es poden utilitzar diversos esquemes per calcular la derivada. Depenent del esquema escollit, la precissió serà millor.

$$f'(a) \approx \frac{f(a+h)-f(a)}{h}$

$$f'(a) \approx \frac{f(a+h)-f(a-h)}{2h}$

% Considerem una funció de la qual la derivada és coneguda i aproximem el
% seu valor
f=@(x) sin(x);
derf=@(x) cos(x);
a=1.2;  %punt en el avaluem la derivada
resultat=[];
for h=[1.e-2, 1.e-3, 1.e-4, 1.e-5] %distancia a la que ens desplacem
    derDavant=(f(a+h)-f(a))/h;
    error1=abs(derDavant-derf(a)); %calculem l'error
    derCentrada=(f(a+h)-f(a-h))/(2*h);
    error2=abs(derCentrada-derf(a)); %calculem l'error
    resultat=[resultat;derDavant,error1,derCentrada,error2];
end
resultat %mostrem els resultats (fixeu-vos en la última columna!!)
resultat =

   3.5769e-01   4.6662e-03   3.6235e-01   6.0393e-06
   3.6189e-01   4.6608e-04   3.6236e-01   6.0393e-08
   3.6231e-01   4.6603e-05   3.6236e-01   6.0362e-10
   3.6235e-01   4.6602e-06   3.6236e-01   2.4330e-12

Derivades Parcials

Si ara considerem una funció f(x,y), els mateixos esquemes de diferències finites es poden aplicar a les derivades parcials, però fent els increments per cada component.

$$\frac{\partial f}{\partial x}(a,b) \approx \frac{f(a+h,b)-f(a-h,b)}{2h}$

$$\frac{\partial f}{\partial y}(a,b) \approx \frac{f(a,b+h)-f(a,b-h)}{2h}$

Equivalentment les seves versions de derivades endavant.

fxy=@(x,y)   sin(x)*cos(y);
parcX=@(x,y) cos(x)*cos(y);
parcY=@(x,y) -sin(x)*sin(y);
a=2;
b=3;
h=1.e-5;
derCentradaX=(fxy(a+h,b)-fxy(a-h,b))/(2*h);
derCentradaY=(fxy(a,b+h)-fxy(a,b-h))/(2*h);
grad=[derCentradaX,derCentradaY] %vector gradient en (a,b)
errors=[abs(derCentradaX-parcX(a,b)),abs(derCentradaY-parcY(a,b))]
grad =

   4.1198e-01  -1.2832e-01


errors =

   2.9377e-12   2.8991e-12

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

Si considerem una imatge com una funció de R^2 que a cada pixel, li assigna un valor d'intensitat, podem derivar numèricament en cada pixel i aquesta derivada serà zero per zones del mateix valor i diferent de zero quan notem un canvi de zona (un contorn d'un objecte!!)

Les derivades les farem per TOTS els pixels. Si la imatge és de color, la passarem primer a escala de grisos.

%Imatge = imread('BB8.jpg');
Imatge = imread('coins.png');
figure('Name','Imatge Original','NumberTitle','on'); %posem títol
imshow(Imatge); % mostrem la imatge

[m,n,p]=size(Imatge);
if ( p > 1) %en el cas d'una imatge RGB, la passem a grisos
    Imatge = rgb2gray(Imatge);
    figure('Name','Imatge en Gris','NumberTitle','on');
    imshow(Imatge);
end
A = im2double(Imatge); %passem els valors d'intensitat a doubles per poder fer càlculs
[m,n]=size(A);

Contorns en X (són els 'horitzontals')

contAx=zeros(size(A)); %matriu de contorns en X
%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
        contAx(i,j)=(A(i+1,j)-A(i-1,j))/2; %la distància entre ixels és h=1
    end
end
figure('Name','Parcial respecte X','NumberTitle','on');
imshow(contAx)

Contorns en Y (són els 'verticals')

contAy=zeros(size(A)); %matriu de contorns en Y
for i=2:m-1
    for j=2:n-1
        contAy(i,j)=(A(i,j+1)-A(i,j-1))/2; %la distància entre ixels és h=1
    end
end
figure('Name','Parcial respecte Y','NumberTitle','on');
imshow(contAy)

Contorns totals associats al gradient en cada pixel

gradA=sqrt(contAx.^2+contAy.^2);
contA=zeros(size(A));
for i=2:m-1
    for j=2:n-1
        dAx=(A(i,j+1)-A(i,j-1))/2;
        dAy=(A(i+1,j)-A(i-1,j))/2;
        contA(i,j)=sqrt(dAx^2+dAy^2);
    end
end
max(gradA(:)-contA(:))
figure('Name','Gradient','NumberTitle','on');
imshow(contA)
ans =

     0

Prewitt filtre per contorns

Per intentar resaltar millor els contorns, Prewitt va proposar de 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.

contAPrew=zeros(size(A));
for i=2:m-1
    for j=2:n-1
        dAx= (A(i+1,j)-A(i-1,j))/2;
        dAx= dAx+(A(i+1,j+1)-A(i-1,j+1))/2;
        dAx= dAx+(A(i+1,j-1)-A(i-1,j-1))/2;

        dAy= (A(i,j+1)-A(i,j-1))/2;
        dAy= dAy+(A(i+1,j+1)-A(i+1,j-1))/2;
        dAy= dAy+(A(i-1,j+1)-A(i-1,j-1))/2;

        contAPrew(i,j)=sqrt(dAx^2+dAy^2);
    end
end
figNumber=figure('Name','Prewitt','NumberTitle','on');
imshow(contAPrew)

Exercici:

Feu el mateix estudi de contorns utilitzant la imatge BB8.jpg (també podeu fer-ho per qualsevol altra imatge)

(c)Numerical Factory 2018