1 files changed, 65 insertions, 0 deletions

diff --git a/SD-VBS/common/toolbox/toolbox_basic/affine/find_D.m b/SD-VBS/common/toolbox/toolbox_basic/affine/find_D.m
new file mode 100755
index 0000000..1e42cb2
--- /dev/null
+++ b/SD-VBS/common/toolbox/toolbox_basic/affine/find_D.m
@@ -0,0 +1,65 @@
+function D = find_D(I,J,mask,w)
+%
+%  function D = find_D(I,J,mask,w)
+%
+%  find the vector D such that it minimizes then
+%  difference between I(Ax-d)-J(x).
+%
+%  mask: the weight matrix,
+%  w: window size for estimating gradiant, use a large value
+%     when A,D are large.
+%
+
+%
+%  NOTE:  Because gradient values on the boarder regions of
+%         I and J can not be computed accuately when using
+%         a gaussian of large support, those boarder regions
+%         of width w are not used in computing D.
+%
+
+%
+%  Jianbo Shi
+%
+
+[gy1,gx1] = grad(I,w);
+[gy2,gx2] = grad(J,w);
+
+gx = 0.5*(gx1+gx2);
+gy = 0.5*(gy1+gy2);
+
+[size_y,size_x] = size(I);
+[center_y,center_x] = find_center(size_y,size_x);
+mask = mask(w+1:size_y-w,w+1:size_x-w);
+
+[x,y] = meshgrid(1:size_x,1:size_y);
+x = x- center_x;
+y = y-center_y;
+
+x = x(w+1:size_y-w,w+1:size_x-w);
+y = y(w+1:size_y-w,w+1:size_x-w);
+
+gx_sqr = gx.*mask.*gx;
+gx_gy = gx.*mask.*gy;
+gy_sqr = gy.*mask.*gy;
+
+
+T= zeros(2,2);
+
+T(1,1) = 0.5*trapz(trapz(gx_sqr));
+T(2,1) = trapz(trapz(gx_gy));
+T(2,2) = 0.5*trapz(trapz(gy_sqr));
+
+T = T+T';
+
+J = J(w+1:size_y-w,w+1:size_x-w);
+I = I(w+1:size_y-w,w+1:size_x-w);
+
+
+diff = (J-I).*mask;
+b(1) = trapz(trapz(diff.*gx));
+b(2) = trapz(trapz(diff.*gy));
+
+a = inv(T)*b';
+
+D= [a(1),a(2)];
+