72 |
|
double DSFoutside = eval_DSFsurround(rbf, outvec, outside_rad); |
73 |
|
#define interp_rad inside_rad + (outside_rad-inside_rad) * \ |
74 |
|
(DSFtarget-DSFinside) / (DSFoutside-DSFinside) |
75 |
< |
/* interpolation search */ |
76 |
< |
while (outside_rad-inside_rad > rad_epsilon) { |
75 |
> |
/* Newton's method (sort of) */ |
76 |
> |
do { |
77 |
|
double test_rad = interp_rad; |
78 |
< |
double DSFtest = eval_DSFsurround(rbf, outvec, test_rad); |
79 |
< |
if (DSFtarget < DSFtest) { |
78 |
> |
double DSFtest; |
79 |
> |
if (test_rad >= outside_rad) |
80 |
> |
return(test_rad); |
81 |
> |
if (test_rad <= inside_rad) |
82 |
> |
return(test_rad*(test_rad>0)); |
83 |
> |
DSFtest = eval_DSFsurround(rbf, outvec, test_rad); |
84 |
> |
if (DSFtest > DSFtarget) { |
85 |
|
inside_rad = test_rad; |
86 |
|
DSFinside = DSFtest; |
87 |
|
} else { |
88 |
|
outside_rad = test_rad; |
89 |
|
DSFoutside = DSFtest; |
90 |
|
} |
91 |
< |
} |
91 |
> |
if (DSFoutside >= DSFinside) |
92 |
> |
return(test_rad); |
93 |
> |
} while (outside_rad-inside_rad > rad_epsilon); |
94 |
|
return(interp_rad); |
95 |
|
#undef interp_rad |
96 |
|
} |