i need determine parameters of illumintaion change, defined continuous piece-wise polynomial c(t), f(t) cubic curve defined 2 boundary points (t1,c) , (t2,0), f'(t1)=0 , f'(t2)=0. original paper: texture-consistent shadow removal
intensity curve sampled normal on boundary of shadow , looks this:
each row sample, displaying illumintaion change.so x number of column , y intensity of pixel.
i have real data (one sample avaraged samples): 
at have n samples , need determine parameters (c,t1,t2)
how can it?
i tried solving linear equation in matlab:
avr_curve average curve, obtained averaging on samples.
f(x)= x^3+a2*x^2+a1*x1+a0 cubic function
%t1,t2 selected hand t1= 10; t2= 15; offset=10; avr_curve= [41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105]; %gradx= convn(avr_curve,[-1 1],'same'); a= zeros(2*offset+1,3); %b= zeros(2*offset+1,1); b= avr_curve'; %for i= 1:2*offset+1 i=t1:t2 x= i-offset-1 a(i,1)= x^2; %a2 a(i,2)= x; %a1 a(i,3)= 1; %a0 b(i,1)= b(i,1)-x^3; end u= a\b; figure,plot(avr_curve(t1:t2)) %estimated cubic curve i= 1:2*offset+1 x= i-offset-1; fx(i)=x^3+u(1)*x^2+u(2)*x+u(3); end figure,plot(fx(t1:t2)) part of avr_curve on [t1 t2] 
cubic curve got (don't looks avr_curve) 
so i'm doing wrong?
update: seems error due model cubic polynomial using 3 variables this:
f(x)= x^3+a2*x^2+a1*x1+a0 - 3 variables but use 4 variables seems ok:
f(x)= a3*x^3+a2*x^2+a1*x1+a0 - 4 variables here code in matlab:
%defined hand t1= 10; t2= 14; avr_curve= [41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105]; x= [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]; %x= [-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; %real x axis %%%model 1 %%f(x)= x^3+a2*x^2+a1*x1+a0 - 3 variables %a= zeros(4,3); %b= [43 104]'; %%cubic equation @ t1 %a(1,1)= t1^2; %a2 %a(1,2)= t1; %a1 %a(1,3)= 1; %a0 %b(1,1)= b(1,1)-t1^3; %%cubic equation @ t2 %a(2,1)= t2^2; %a2 %a(2,2)= t2; %a1 %a(2,3)= 1; %a0 %b(2,1)= b(2,1)-t1^3; %%1st derivative @ t1 %a(3,1)= 2*t1; %a2 %a(3,2)= 1; %a1 %a(3,3)= 0; %a0 %b(3,1)= -3*t1^2; %%1st derivative @ t2 %a(4,1)= 2*t2; %a2 %a(4,2)= 1; %a1 %a(4,3)= 0; %a0 %b(4,1)= -3*t2^2; %model 2 %f(x)= a3*x^3+a2*x^2+a1*x1+a0 - 4 variables a= zeros(4,4); b= [43 104]'; %cubic equation @ t1 a(1,1)= t1^3; %a3 a(1,2)= t1^2; %a2 a(1,3)= t1; %a1 a(1,4)= 1; %a0 b(1,1)= b(1,1); %cubic equation @ t2 a(2,1)= t2^3; %a3 a(2,2)= t2^2; %a2 a(2,3)= t2; %a1 a(2,4)= 1; %a0 b(2,1)= b(2,1); %1st derivative @ t1 a(3,1)= 3*t1^2; %a3 a(3,2)= 2*t1; %a2 a(3,3)= 1; %a1 a(3,4)= 0; %a0 b(3,1)= 0; %1st derivative @ t2 a(4,1)= 3*t2^2; %a3 a(4,2)= 2*t2; %a2 a(4,3)= 1; %a1 a(4,4)= 0; %a0 b(4,1)= 0; size(a) size(b) u= a\b; u %estimated cubic curve %dx=[1:21]; % global view dx=t1-1:t2+1; % local view in [t1 t2] x= dx %fx(x)=x^3+u(1)*x^2+u(2)*x+u(3); % model 1 fx(x)= u(1)*x^3+u(2)*x^2+u(3)*x+u(4); % model 2 end err= 0; x= dx err= err+(fx(x)-avr_curve(x))^2; end err figure,plot(dx,avr_curve(dx),dx,fx(dx)) spline on interval [t1-1 t2+1] 
and on full interval
disclaimer
i cannot give guarantees on correctness of code or methods given below, use critical sense before using of that.
0. define problem
you have piecewise defined function
where f(t) cubic function, in order uniquely identify it, following additional conditions given
you want recover best values of parameters t1, t2 , sigma minimize error given set of points.
this curve fitting in least squares sense.
1 parametrize f(t) cubic function
in order compute error between candidate cl(t) function , set of points need compute f(t), general form (being cubic) is
so seems have 4 additional parameters consider. indeed parameters totally defined free 3 parameters t1, t2 , sigma.
important not confuse free parameters dependent ones.
given additional conditions on f(t) can set linear system
which has 1 solution (as expected) given by
this tell how compute parameters of cubic given 3 free parameters.
way cl(t) determined, it's time find best candidate.
2 minimize error
i go least squares now.
since not linear function, there no closed form computing sigma, t1 , t2.
there numerical methods, gauss-newton one.
however 1 way or required compute partial derivatives respect of 3 parameters.
i don't know how compute derivative respect of separation parameter t1.
i've searched mathse , found this question address same problem, nobody answered.
without partial derivatives least squares methods over.
so take more practical road , implemented brute force function in c try every possible triplet of parameter , return best match.
3 brute force function
given nature of problem, turned out o(n^2) in number of sample.
the algorithm proceeds follow: divide sample set in 3 parts, first 1 part of point before t1, second 1 of points between t1 , t2 , last 1 of points after t2.
the first part used compute sigma, sigma arithmetic average of points in part 1.
t1 , t2 computed through cycle, t1 set every possible point in original points set, starting second , going forward.
for every choice of t1, t2 set every possible point after t1.
at each iteration error computed , if minimum ever seen, parameters used saved.
the error computer absolute value of residuals since absolute value should fast (surely faster square) , fits purpose of metric.
4 code
#include <stdio.h> #include <math.h> float point_on_curve(float sigma, float t1, float t2, float t) { float a,b,c,d, k; if (t <= t1) return sigma; if (t >= t2) return 0; k = (t1-t2)*(t1-t2)*(t1-t2); = -2*sigma/k; b = 3*sigma*(t1+t2)/k; c = -6*sigma*t1*t2/k; d = sigma*t2*t2*(3*t1-t2)/k; return a*t*t*t + b*t*t + c*t + d; } float compute_error(float sigma, float t1, float t2, int s, int dx, int* data, int len) { float error=0; unsigned int i; (i = 0; < len; i++) error += fabs(point_on_curve(sigma, t1, t2, s+i*dx)- data[i]); return error; } /* * s starting time of samples set * dx separation in time between 2 sample (a.k.a. sampling period) * data array of samples * len number of samples * sigma, t1, t2 pointers output parameters computed * * return 1 if not enough (3) samples, 0 if went ok. */ int curve_fit(int s, int dx, int* data, unsigned int len, float* sigma, float* t1, float* t2) { float l_sigma = 0; float l_t1, l_t2; float sum = 0; float min_error, cur_error; char error_valid = 0; unsigned int i, j; if (len < 3) return 1; (i = 0; < len; i++) { /* compute sigma average of points <= */ sum += data[i]; l_sigma = sum/(i+1); /* set t1 point i+1 */ l_t1 = s+(i+1)*dx; (j = i+2; j < len; j++) { /* set t2 points i+2, i+3, i+4, ... */ l_t2 = s+j*dx; /* compute error */ cur_error = compute_error(l_sigma, l_t1, l_t2, s, dx, data, len); if (cur_error < min_error || !error_valid) { error_valid = 1; min_error = cur_error; *sigma = l_sigma; *t1 = l_t1; *t2 = l_t2; } } } return 0; } int main() { float sigma, t1, t2; int data[]={41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105}; unsigned int len = sizeof(data)/sizeof(int); unsigned int i; (i = 0; < len; i++) data[i] -= 107; /* subtract max */ if (curve_fit(1,1,data, len, &sigma, &t1, &t2)) printf("not enough data!\n"); else printf("parameters: sigma = %.3f, t1 = %.3f, t2 = %.3f\n", sigma, t1, t2); return 0; } please note cl(t) defined having 0 right limit, code assume case.
this why max value (107) subtracted every sample, have worked definition of cl(t) given @ beginning , late noted sample biased.
by i'm not going adapt code, can add parameter in problem , redo steps 1 if needed, or translate samples using maximum value.
the output of code is
parameters: sigma = -65.556, t1 = 10.000, t2 = 14.000 which match points set given, considering vertically translated -107.
Comments
Post a Comment