[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/multi_convolution.hxx | ![]() |
00001 //-- -*- c++ -*- 00002 /************************************************************************/ 00003 /* */ 00004 /* Copyright 2003 by Christian-Dennis Rahn */ 00005 /* and Ullrich Koethe */ 00006 /* */ 00007 /* This file is part of the VIGRA computer vision library. */ 00008 /* The VIGRA Website is */ 00009 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00012 /* vigra@informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 #ifndef VIGRA_MULTI_CONVOLUTION_H 00039 #define VIGRA_MULTI_CONVOLUTION_H 00040 00041 #include "separableconvolution.hxx" 00042 #include "array_vector.hxx" 00043 #include "multi_array.hxx" 00044 #include "accessor.hxx" 00045 #include "numerictraits.hxx" 00046 #include "navigator.hxx" 00047 #include "metaprogramming.hxx" 00048 #include "multi_pointoperators.hxx" 00049 #include "functorexpression.hxx" 00050 #include "tinyvector.hxx" 00051 #include "algorithm.hxx" 00052 00053 namespace vigra 00054 { 00055 00056 namespace detail 00057 { 00058 00059 struct DoubleYielder 00060 { 00061 const double value; 00062 DoubleYielder(double v, unsigned, const char *const) : value(v) {} 00063 DoubleYielder(double v) : value(v) {} 00064 void operator++() {} 00065 double operator*() const { return value; } 00066 }; 00067 00068 template <typename X> 00069 struct IteratorDoubleYielder 00070 { 00071 X it; 00072 IteratorDoubleYielder(X i, unsigned, const char *const) : it(i) {} 00073 IteratorDoubleYielder(X i) : it(i) {} 00074 void operator++() { ++it; } 00075 double operator*() const { return *it; } 00076 }; 00077 00078 template <typename X> 00079 struct SequenceDoubleYielder 00080 { 00081 typename X::const_iterator it; 00082 SequenceDoubleYielder(const X & seq, unsigned dim, 00083 const char *const function_name = "SequenceDoubleYielder") 00084 : it(seq.begin()) 00085 { 00086 if (seq.size() == dim) 00087 return; 00088 std::string msg = "(): Parameter number be equal to the number of spatial dimensions."; 00089 vigra_precondition(false, function_name + msg); 00090 } 00091 void operator++() { ++it; } 00092 double operator*() const { return *it; } 00093 }; 00094 00095 template <typename X> 00096 struct WrapDoubleIterator 00097 { 00098 typedef 00099 typename IfBool< IsConvertibleTo<X, double>::value, 00100 DoubleYielder, 00101 typename IfBool< IsIterator<X>::value || IsArray<X>::value, 00102 IteratorDoubleYielder<X>, 00103 SequenceDoubleYielder<X> 00104 >::type 00105 >::type type; 00106 }; 00107 00108 template <class Param1, class Param2, class Param3> 00109 struct WrapDoubleIteratorTriple 00110 { 00111 typename WrapDoubleIterator<Param1>::type sigma_eff_it; 00112 typename WrapDoubleIterator<Param2>::type sigma_d_it; 00113 typename WrapDoubleIterator<Param3>::type step_size_it; 00114 WrapDoubleIteratorTriple(Param1 sigma_eff, Param2 sigma_d, Param3 step_size) 00115 : sigma_eff_it(sigma_eff), sigma_d_it(sigma_d), step_size_it(step_size) {} 00116 void operator++() 00117 { 00118 ++sigma_eff_it; 00119 ++sigma_d_it; 00120 ++step_size_it; 00121 } 00122 double sigma_eff() const { return *sigma_eff_it; } 00123 double sigma_d() const { return *sigma_d_it; } 00124 double step_size() const { return *step_size_it; } 00125 static double sqr(double x) { return x * x; } 00126 static void sigma_precondition(double sigma, const char *const function_name) 00127 { 00128 if (sigma < 0.0) 00129 { 00130 std::string msg = "(): Scale must be positive."; 00131 vigra_precondition(false, function_name + msg); 00132 } 00133 } 00134 double sigma_scaled(const char *const function_name = "unknown function ") const 00135 { 00136 sigma_precondition(sigma_eff(), function_name); 00137 sigma_precondition(sigma_d(), function_name); 00138 double sigma_squared = sqr(sigma_eff()) - sqr(sigma_d()); 00139 if (sigma_squared > 0.0) 00140 { 00141 return std::sqrt(sigma_squared) / step_size(); 00142 } 00143 else 00144 { 00145 std::string msg = "(): Scale would be imaginary or zero."; 00146 vigra_precondition(false, function_name + msg); 00147 return 0; 00148 } 00149 } 00150 }; 00151 00152 template <unsigned dim> 00153 struct multiArrayScaleParam 00154 { 00155 typedef TinyVector<double, dim> p_vector; 00156 typedef typename p_vector::const_iterator return_type; 00157 p_vector vec; 00158 00159 template <class Param> 00160 multiArrayScaleParam(Param val, const char *const function_name = "multiArrayScaleParam") 00161 { 00162 typename WrapDoubleIterator<Param>::type in(val, dim, function_name); 00163 for (unsigned i = 0; i != dim; ++i, ++in) 00164 vec[i] = *in; 00165 } 00166 return_type operator()() const 00167 { 00168 return vec.begin(); 00169 } 00170 static void precondition(unsigned n_par, const char *const function_name = "multiArrayScaleParam") 00171 { 00172 char n[3] = "0."; 00173 n[0] += dim; 00174 std::string msg = "(): dimension parameter must be "; 00175 vigra_precondition(dim == n_par, function_name + msg + n); 00176 } 00177 multiArrayScaleParam(double v0, double v1, const char *const function_name = "multiArrayScaleParam") 00178 { 00179 precondition(2, function_name); 00180 vec = p_vector(v0, v1); 00181 } 00182 multiArrayScaleParam(double v0, double v1, double v2, const char *const function_name = "multiArrayScaleParam") 00183 { 00184 precondition(3, function_name); 00185 vec = p_vector(v0, v1, v2); 00186 } 00187 multiArrayScaleParam(double v0, double v1, double v2, double v3, const char *const function_name = "multiArrayScaleParam") 00188 { 00189 precondition(4, function_name); 00190 vec = p_vector(v0, v1, v2, v3); 00191 } 00192 multiArrayScaleParam(double v0, double v1, double v2, double v3, double v4, const char *const function_name = "multiArrayScaleParam") 00193 { 00194 precondition(5, function_name); 00195 vec = p_vector(v0, v1, v2, v3, v4); 00196 } 00197 }; 00198 00199 } // namespace detail 00200 00201 #define VIGRA_CONVOLUTION_OPTIONS(function_name, default_value, member_name) \ 00202 template <class Param> \ 00203 ConvolutionOptions & function_name(const Param & val) \ 00204 { \ 00205 member_name = ParamVec(val, "ConvolutionOptions::" #function_name); \ 00206 return *this; \ 00207 } \ 00208 ConvolutionOptions & function_name() \ 00209 { \ 00210 member_name = ParamVec(default_value, "ConvolutionOptions::" #function_name); \ 00211 return *this; \ 00212 } \ 00213 ConvolutionOptions & function_name(double v0, double v1) \ 00214 { \ 00215 member_name = ParamVec(v0, v1, "ConvolutionOptions::" #function_name); \ 00216 return *this; \ 00217 } \ 00218 ConvolutionOptions & function_name(double v0, double v1, double v2) \ 00219 { \ 00220 member_name = ParamVec(v0, v1, v2, "ConvolutionOptions::" #function_name); \ 00221 return *this; \ 00222 } \ 00223 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3) \ 00224 { \ 00225 member_name = ParamVec(v0, v1, v2, v3, "ConvolutionOptions::" #function_name); \ 00226 return *this; \ 00227 } \ 00228 ConvolutionOptions & function_name(double v0, double v1, double v2, double v3, double v4) \ 00229 { \ 00230 member_name = ParamVec(v0, v1, v2, v3, v4, "ConvolutionOptions::" #function_name); \ 00231 return *this; \ 00232 } 00233 00234 00235 /** \brief Options class template for convolutions. 00236 00237 <b>\#include</b> <vigra/multi_convolution.hxx> 00238 00239 This class enables the calculation of scale space convolutions 00240 such as \ref gaussianGradientMultiArray() on data with anisotropic 00241 discretization. For these, the result of the ordinary calculation 00242 has to be multiplied by factors of \f$1/w^{n}\f$ for each dimension, 00243 where \f$w\f$ is the step size of the grid in said dimension and 00244 \f$n\f$ is the differential order of the convolution, e.g., 1 for 00245 gaussianGradientMultiArray(), and 0 for gaussianSmoothMultiArray(), 00246 respectively. Also for each dimension in turn, the convolution's scale 00247 parameter \f$\sigma\f$ has to be replaced by 00248 \f$\sqrt{\sigma_\mathrm{eff}^2 - \sigma_\mathrm{D}^2}\Big/w\f$, 00249 where \f$\sigma_\mathrm{eff}\f$ is the resulting effective filtering 00250 scale. The data is assumed to be already filtered by a 00251 gaussian smoothing with the scale parameter \f$\sigma_\mathrm{D}\f$ 00252 (such as by measuring equipment). All of the above changes are 00253 automatically employed by the convolution functions for <tt>MultiArray</tt>s 00254 if a corresponding options object is provided. 00255 00256 The <tt>ConvolutionOptions</tt> class must be parameterized by the dimension 00257 <tt>dim</tt> 00258 of the <tt>MultiArray</tt>s on which it is used. The actual per-axis 00259 options are set by (overloaded) member functions explained below, 00260 or else default to neutral values corresponding to the absence of the 00261 particular option. 00262 00263 All member functions set <tt>dim</tt> values of the respective convolution 00264 option, one for each dimension. They may be set explicitly by multiple 00265 arguments for up to five dimensions, or by a single argument to the same 00266 value for all dimensions. For the general case, a single argument that is 00267 either a C-syle array, an iterator, or a C++ standard library style 00268 sequence (such as <tt>std::vector</tt>, with member functions <tt>begin()</tt> 00269 and <tt>size()</tt>) supplies the option values for any number of dimensions. 00270 00271 Note that the return value of all member functions is <tt>*this</tt>, which 00272 provides the mechanism for concatenating member function calls as shown below. 00273 00274 <b>usage with explicit parameters:</b> 00275 00276 \code 00277 ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(1, 2.3); 00278 \endcode 00279 00280 <b>usage with arrays:</b> 00281 00282 \code 00283 const double step_size[3] = { x_scale, y_scale, z_scale }; 00284 ConvolutionOptions<3> opt = ConvolutionOptions<3>().stepSize(step_size); 00285 \endcode 00286 00287 <b>usage with C++ standard library style sequences:</b> 00288 00289 \code 00290 TinyVector<double, 4> step_size(1, 1, 2.0, 1.5); 00291 TinyVector<double, 4> r_sigmas(1, 1, 2.3, 3.2); 00292 ConvolutionOptions<4> opt = ConvolutionOptions<4>().stepSize(step_size).resolutionStdDev(r_sigmas); 00293 \endcode 00294 00295 <b>usage with iterators:</b> 00296 00297 \code 00298 ArrayVector<double> step_size; 00299 step_size.push_back(0); 00300 step_size.push_back(3); 00301 step_size.push_back(4); 00302 ArrayVector<double>::iterator i = step_size.begin(); 00303 ++i; 00304 ConvolutionOptions<2> opt = ConvolutionOptions<2>().stepSize(i); 00305 \endcode 00306 00307 <b>general usage in a convolution function call:</b> 00308 00309 \code 00310 MultiArray<3, double> test_image; 00311 MultiArray<3, double> out_image; 00312 gaussianSmoothMultiArray(srcMultiArrayRange(test_image), 00313 destMultiArray(out_image), 00314 5.0, 00315 ConvolutionOptions<3>() 00316 .stepSize (1, 1, 3.2) 00317 .resolutionStdDev(1, 1, 4) 00318 ); 00319 \endcode 00320 00321 */ 00322 template <unsigned dim> 00323 class ConvolutionOptions 00324 { 00325 public: 00326 typedef typename MultiArrayShape<dim>::type Shape; 00327 typedef detail::multiArrayScaleParam<dim> ParamVec; 00328 typedef typename ParamVec::return_type ParamIt; 00329 00330 ParamVec sigma_eff; 00331 ParamVec sigma_d; 00332 ParamVec step_size; 00333 ParamVec outer_scale; 00334 double window_ratio; 00335 Shape from_point, to_point; 00336 00337 ConvolutionOptions() 00338 : sigma_eff(0.0), 00339 sigma_d(0.0), 00340 step_size(1.0), 00341 outer_scale(0.0), 00342 window_ratio(0.0) 00343 {} 00344 00345 typedef typename detail::WrapDoubleIteratorTriple<ParamIt, ParamIt, ParamIt> 00346 ScaleIterator; 00347 typedef typename detail::WrapDoubleIterator<ParamIt>::type 00348 StepIterator; 00349 00350 ScaleIterator scaleParams() const 00351 { 00352 return ScaleIterator(sigma_eff(), sigma_d(), step_size()); 00353 } 00354 StepIterator stepParams() const 00355 { 00356 return StepIterator(step_size()); 00357 } 00358 00359 ConvolutionOptions outerOptions() const 00360 { 00361 ConvolutionOptions outer = *this; 00362 // backward-compatible values: 00363 return outer.stdDev(outer_scale()).resolutionStdDev(0.0); 00364 } 00365 00366 // Step size per axis. 00367 // Default: dim values of 1.0 00368 VIGRA_CONVOLUTION_OPTIONS(stepSize, 1.0, step_size) 00369 #ifdef DOXYGEN 00370 /** Step size(s) per axis, i.e., the distance between two 00371 adjacent pixels. Required for <tt>MultiArray</tt> 00372 containing anisotropic data. 00373 00374 Note that a convolution containing a derivative operator 00375 of order <tt>n</tt> results in a multiplication by 00376 \f${\rm stepSize}^{-n}\f$ for each axis. 00377 Also, the above standard deviations 00378 are scaled according to the step size of each axis. 00379 Default value for the options object if this member function is not 00380 used: A value of 1.0 for each dimension. 00381 */ 00382 ConvolutionOptions<dim> & stepSize(...); 00383 #endif 00384 00385 // Resolution standard deviation per axis. 00386 // Default: dim values of 0.0 00387 VIGRA_CONVOLUTION_OPTIONS(resolutionStdDev, 0.0, sigma_d) 00388 #ifdef DOXYGEN 00389 /** Resolution standard deviation(s) per axis, i.e., a supposed 00390 pre-existing gaussian filtering by this value. 00391 00392 The standard deviation actually used by the convolution operators 00393 is \f$\sqrt{{\rm sigma}^{2} - {\rm resolutionStdDev}^{2}}\f$ for each 00394 axis. 00395 Default value for the options object if this member function is not 00396 used: A value of 0.0 for each dimension. 00397 */ 00398 ConvolutionOptions<dim> & resolutionStdDev(...); 00399 #endif 00400 00401 // Standard deviation of scale space operators. 00402 // Default: dim values of 0.0 00403 VIGRA_CONVOLUTION_OPTIONS(stdDev, 0.0, sigma_eff) 00404 VIGRA_CONVOLUTION_OPTIONS(innerScale, 0.0, sigma_eff) 00405 00406 #ifdef DOXYGEN 00407 /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray(). 00408 00409 Usually not 00410 needed, since a single value for all axes may be specified as a parameter 00411 <tt>sigma</tt> to the call of 00412 an convolution operator such as \ref gaussianGradientMultiArray(), and 00413 anisotropic data requiring the use of the stepSize() member function. 00414 Default value for the options object if this member function is not 00415 used: A value of 0.0 for each dimension. 00416 */ 00417 ConvolutionOptions<dim> & stdDev(...); 00418 00419 /** Standard deviation(s) of scale space operators, or inner scale(s) for \ref structureTensorMultiArray(). 00420 00421 Usually not 00422 needed, since a single value for all axes may be specified as a parameter 00423 <tt>sigma</tt> to the call of 00424 an convolution operator such as \ref gaussianGradientMultiArray(), and 00425 anisotropic data requiring the use of the stepSize() member function. 00426 Default value for the options object if this member function is not 00427 used: A value of 0.0 for each dimension. 00428 */ 00429 ConvolutionOptions<dim> & innerScale(...); 00430 #endif 00431 00432 // Outer scale, for structure tensor. 00433 // Default: dim values of 0.0 00434 VIGRA_CONVOLUTION_OPTIONS(outerScale, 0.0, outer_scale) 00435 #ifdef DOXYGEN 00436 /** Standard deviation(s) of the second convolution of the 00437 structure tensor. 00438 00439 Usually not needed, since a single value for 00440 all axes may be specified as a parameter <tt>outerScale</tt> to 00441 the call of \ref structureTensorMultiArray(), and 00442 anisotropic data requiring the use of the stepSize() member 00443 function. 00444 Default value for the options object if this member function is not 00445 used: A value of 0.0 for each dimension. 00446 */ 00447 ConvolutionOptions<dim> & outerScale(...); 00448 #endif 00449 00450 /** Size of the filter window as a multiple of the scale parameter. 00451 00452 This option is only used for Gaussian filters and their derivatives. 00453 By default, the window size of a Gaussian filter is automatically 00454 determined such that the error resulting from restricting the 00455 infinitely large Gaussian function to a finite size is minimized. 00456 In particular, the window radius is determined as 00457 <tt>radius = round(3.0 * sigma + 0.5 * order)</tt>, where 'order' is the 00458 desired derivative order. In some cases, it is desirable to trade off 00459 accuracy for speed, and this function can be used to request a smaller 00460 window radius. 00461 00462 Default: <tt>0.0</tt> (i.e. determine the window size automatically) 00463 */ 00464 ConvolutionOptions<dim> & filterWindowSize(double ratio) 00465 { 00466 vigra_precondition(ratio >= 0.0, 00467 "ConvolutionOptions::filterWindowSize(): ratio must not be negative."); 00468 window_ratio = ratio; 00469 return *this; 00470 } 00471 00472 /** Restrict the filter to a subregion of the input array. 00473 00474 This is useful for speeding up computations by ignoring irrelevant 00475 areas in the array. <b>Note:</b> It is assumed that the output array 00476 of the convolution has the size given in this function. 00477 00478 Default: <tt>from = Shape(), to = Shape()</tt> (i.e. use entire array) 00479 */ 00480 ConvolutionOptions<dim> & subarray(Shape const & from, Shape const & to) 00481 { 00482 from_point = from; 00483 to_point = to; 00484 return *this; 00485 } 00486 }; 00487 00488 namespace detail 00489 { 00490 00491 /********************************************************/ 00492 /* */ 00493 /* internalSeparableConvolveMultiArray */ 00494 /* */ 00495 /********************************************************/ 00496 00497 template <class SrcIterator, class SrcShape, class SrcAccessor, 00498 class DestIterator, class DestAccessor, class KernelIterator> 00499 void 00500 internalSeparableConvolveMultiArrayTmp( 00501 SrcIterator si, SrcShape const & shape, SrcAccessor src, 00502 DestIterator di, DestAccessor dest, KernelIterator kit) 00503 { 00504 enum { N = 1 + SrcIterator::level }; 00505 00506 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType; 00507 typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor; 00508 00509 // temporary array to hold the current line to enable in-place operation 00510 ArrayVector<TmpType> tmp( shape[0] ); 00511 00512 typedef MultiArrayNavigator<SrcIterator, N> SNavigator; 00513 typedef MultiArrayNavigator<DestIterator, N> DNavigator; 00514 00515 TmpAcessor acc; 00516 00517 { 00518 // only operate on first dimension here 00519 SNavigator snav( si, shape, 0 ); 00520 DNavigator dnav( di, shape, 0 ); 00521 00522 for( ; snav.hasMore(); snav++, dnav++ ) 00523 { 00524 // first copy source to tmp for maximum cache efficiency 00525 copyLine(snav.begin(), snav.end(), src, tmp.begin(), acc); 00526 00527 convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc), 00528 destIter( dnav.begin(), dest ), 00529 kernel1d( *kit ) ); 00530 } 00531 ++kit; 00532 } 00533 00534 // operate on further dimensions 00535 for( int d = 1; d < N; ++d, ++kit ) 00536 { 00537 DNavigator dnav( di, shape, d ); 00538 00539 tmp.resize( shape[d] ); 00540 00541 for( ; dnav.hasMore(); dnav++ ) 00542 { 00543 // first copy source to tmp since convolveLine() cannot work in-place 00544 copyLine(dnav.begin(), dnav.end(), dest, tmp.begin(), acc); 00545 00546 convolveLine(srcIterRange(tmp.begin(), tmp.end(), acc), 00547 destIter( dnav.begin(), dest ), 00548 kernel1d( *kit ) ); 00549 } 00550 } 00551 } 00552 00553 /********************************************************/ 00554 /* */ 00555 /* internalSeparableConvolveSubarray */ 00556 /* */ 00557 /********************************************************/ 00558 00559 template <class SrcIterator, class SrcShape, class SrcAccessor, 00560 class DestIterator, class DestAccessor, class KernelIterator> 00561 void 00562 internalSeparableConvolveSubarray( 00563 SrcIterator si, SrcShape const & shape, SrcAccessor src, 00564 DestIterator di, DestAccessor dest, KernelIterator kit, 00565 SrcShape const & start, SrcShape const & stop) 00566 { 00567 enum { N = 1 + SrcIterator::level }; 00568 00569 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType; 00570 typedef MultiArray<N, TmpType> TmpArray; 00571 typedef typename TmpArray::traverser TmpIterator; 00572 typedef typename AccessorTraits<TmpType>::default_accessor TmpAcessor; 00573 00574 SrcShape sstart, sstop, axisorder, tmpshape; 00575 TinyVector<double, N> overhead; 00576 for(int k=0; k<N; ++k) 00577 { 00578 sstart[k] = start[k] - kit[k].right(); 00579 if(sstart[k] < 0) 00580 sstart[k] = 0; 00581 sstop[k] = stop[k] - kit[k].left(); 00582 if(sstop[k] > shape[k]) 00583 sstop[k] = shape[k]; 00584 overhead[k] = double(sstop[k] - sstart[k]) / (stop[k] - start[k]); 00585 } 00586 00587 indexSort(overhead.begin(), overhead.end(), axisorder.begin(), std::greater<double>()); 00588 00589 SrcShape dstart, dstop(sstop - sstart); 00590 dstop[axisorder[0]] = stop[axisorder[0]] - start[axisorder[0]]; 00591 00592 // temporary array to hold the current line to enable in-place operation 00593 MultiArray<N, TmpType> tmp(dstop); 00594 00595 typedef MultiArrayNavigator<SrcIterator, N> SNavigator; 00596 typedef MultiArrayNavigator<TmpIterator, N> TNavigator; 00597 typedef MultiArrayNavigator<DestIterator, N> DNavigator; 00598 00599 TmpAcessor acc; 00600 00601 { 00602 // only operate on first dimension here 00603 SNavigator snav( si, sstart, sstop, axisorder[0]); 00604 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[0]); 00605 00606 ArrayVector<TmpType> tmpline(sstop[axisorder[0]] - sstart[axisorder[0]]); 00607 00608 int lstart = start[axisorder[0]] - sstart[axisorder[0]]; 00609 int lstop = lstart + (stop[axisorder[0]] - start[axisorder[0]]); 00610 00611 for( ; snav.hasMore(); snav++, tnav++ ) 00612 { 00613 // first copy source to tmp for maximum cache efficiency 00614 copyLine(snav.begin(), snav.end(), src, tmpline.begin(), acc); 00615 00616 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc), 00617 destIter(tnav.begin(), acc), 00618 kernel1d( kit[axisorder[0]] ), lstart, lstop); 00619 } 00620 } 00621 00622 // operate on further dimensions 00623 for( int d = 1; d < N; ++d) 00624 { 00625 TNavigator tnav( tmp.traverser_begin(), dstart, dstop, axisorder[d]); 00626 00627 ArrayVector<TmpType> tmpline(dstop[axisorder[d]] - dstart[axisorder[d]]); 00628 00629 int lstart = start[axisorder[d]] - sstart[axisorder[d]]; 00630 int lstop = lstart + (stop[axisorder[d]] - start[axisorder[d]]); 00631 00632 for( ; tnav.hasMore(); tnav++ ) 00633 { 00634 // first copy source to tmp because convolveLine() cannot work in-place 00635 copyLine(tnav.begin(), tnav.end(), acc, tmpline.begin(), acc ); 00636 00637 convolveLine(srcIterRange(tmpline.begin(), tmpline.end(), acc), 00638 destIter( tnav.begin() + lstart, acc ), 00639 kernel1d( kit[axisorder[d]] ), lstart, lstop); 00640 } 00641 00642 dstart[axisorder[d]] = lstart; 00643 dstop[axisorder[d]] = lstop; 00644 } 00645 00646 copyMultiArray(tmp.traverser_begin()+dstart, stop-start, acc, di, dest); 00647 } 00648 00649 00650 template <class K> 00651 void 00652 scaleKernel(K & kernel, double a) 00653 { 00654 for(int i = kernel.left(); i <= kernel.right(); ++i) 00655 kernel[i] = detail::RequiresExplicitCast<typename K::value_type>::cast(kernel[i] * a); 00656 } 00657 00658 00659 } // namespace detail 00660 00661 /** \addtogroup MultiArrayConvolutionFilters Convolution filters for multi-dimensional arrays. 00662 00663 These functions realize a separable convolution on an arbitrary dimensional 00664 array that is specified by iterators (compatible to \ref MultiIteratorPage) 00665 and shape objects. It can therefore be applied to a wide range of data structures 00666 (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.). 00667 */ 00668 //@{ 00669 00670 /********************************************************/ 00671 /* */ 00672 /* separableConvolveMultiArray */ 00673 /* */ 00674 /********************************************************/ 00675 00676 /** \brief Separated convolution on multi-dimensional arrays. 00677 00678 This function computes a separated convolution on all dimensions 00679 of the given multi-dimensional array. Both source and destination 00680 arrays are represented by iterators, shape objects and accessors. 00681 The destination array is required to already have the correct size. 00682 00683 There are two variants of this functions: one takes a single kernel 00684 of type \ref vigra::Kernel1D which is then applied to all dimensions, 00685 whereas the other requires an iterator referencing a sequence of 00686 \ref vigra::Kernel1D objects, one for every dimension of the data. 00687 Then the first kernel in this sequence is applied to the innermost 00688 dimension (e.g. the x-dimension of an image), while the last is applied to the 00689 outermost dimension (e.g. the z-dimension in a 3D image). 00690 00691 This function may work in-place, which means that <tt>siter == diter</tt> is allowed. 00692 A full-sized internal array is only allocated if working on the destination 00693 array directly would cause round-off errors (i.e. if 00694 <tt>typeid(typename NumericTraits<typename DestAccessor::value_type>::RealPromote) 00695 != typeid(typename DestAccessor::value_type)</tt>. 00696 00697 If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent 00698 a valid subarray of the input array. The convolution is then restricted to that 00699 subarray, and it is assumed that the output array only refers to the 00700 subarray (i.e. <tt>diter</tt> points to the element corresponding to 00701 <tt>start</tt>). 00702 00703 <b> Declarations:</b> 00704 00705 pass arguments explicitly: 00706 \code 00707 namespace vigra { 00708 // apply the same kernel to all dimensions 00709 template <class SrcIterator, class SrcShape, class SrcAccessor, 00710 class DestIterator, class DestAccessor, class T> 00711 void 00712 separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 00713 DestIterator diter, DestAccessor dest, 00714 Kernel1D<T> const & kernel, 00715 SrcShape const & start = SrcShape(), 00716 SrcShape const & stop = SrcShape()); 00717 00718 // apply each kernel from the sequence 'kernels' in turn 00719 template <class SrcIterator, class SrcShape, class SrcAccessor, 00720 class DestIterator, class DestAccessor, class KernelIterator> 00721 void 00722 separableConvolveMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 00723 DestIterator diter, DestAccessor dest, 00724 KernelIterator kernels, 00725 SrcShape const & start = SrcShape(), 00726 SrcShape const & stop = SrcShape()); 00727 } 00728 \endcode 00729 00730 use argument objects in conjunction with \ref ArgumentObjectFactories : 00731 \code 00732 namespace vigra { 00733 // apply the same kernel to all dimensions 00734 template <class SrcIterator, class SrcShape, class SrcAccessor, 00735 class DestIterator, class DestAccessor, class T> 00736 void 00737 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 00738 pair<DestIterator, DestAccessor> const & dest, 00739 Kernel1D<T> const & kernel, 00740 SrcShape const & start = SrcShape(), 00741 SrcShape const & stop = SrcShape()); 00742 00743 // apply each kernel from the sequence 'kernels' in turn 00744 template <class SrcIterator, class SrcShape, class SrcAccessor, 00745 class DestIterator, class DestAccessor, class KernelIterator> 00746 void 00747 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 00748 pair<DestIterator, DestAccessor> const & dest, 00749 KernelIterator kernels, 00750 SrcShape const & start = SrcShape(), 00751 SrcShape const & stop = SrcShape()); 00752 } 00753 \endcode 00754 00755 <b> Usage:</b> 00756 00757 <b>\#include</b> <vigra/multi_convolution.hxx> 00758 00759 \code 00760 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 00761 MultiArray<3, unsigned char> source(shape); 00762 MultiArray<3, float> dest(shape); 00763 ... 00764 Kernel1D<float> gauss; 00765 gauss.initGaussian(sigma); 00766 // create 3 Gauss kernels, one for each dimension 00767 ArrayVector<Kernel1D<float> > kernels(3, gauss); 00768 00769 // perform Gaussian smoothing on all dimensions 00770 separableConvolveMultiArray(srcMultiArrayRange(source), destMultiArray(dest), 00771 kernels.begin()); 00772 \endcode 00773 00774 <b> Required Interface:</b> 00775 00776 see \ref separableConvolveMultiArray(), in addition: 00777 00778 \code 00779 int dimension = 0; 00780 VectorElementAccessor<DestAccessor> elementAccessor(0, dest); 00781 \endcode 00782 00783 \see vigra::Kernel1D, convolveLine() 00784 */ 00785 doxygen_overloaded_function(template <...> void separableConvolveMultiArray) 00786 00787 template <class SrcIterator, class SrcShape, class SrcAccessor, 00788 class DestIterator, class DestAccessor, class KernelIterator> 00789 void 00790 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src, 00791 DestIterator d, DestAccessor dest, 00792 KernelIterator kernels, 00793 SrcShape const & start = SrcShape(), 00794 SrcShape const & stop = SrcShape()) 00795 { 00796 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType; 00797 00798 if(stop != SrcShape()) 00799 { 00800 enum { N = 1 + SrcIterator::level }; 00801 for(int k=0; k<N; ++k) 00802 vigra_precondition(0 <= start[k] && start[k] < stop[k] && stop[k] <= shape[k], 00803 "separableConvolveMultiArray(): invalid subarray shape."); 00804 00805 detail::internalSeparableConvolveSubarray(s, shape, src, d, dest, kernels, start, stop); 00806 } 00807 else if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult) 00808 { 00809 // need a temporary array to avoid rounding errors 00810 MultiArray<SrcShape::static_size, TmpType> tmpArray(shape); 00811 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, 00812 tmpArray.traverser_begin(), typename AccessorTraits<TmpType>::default_accessor(), kernels ); 00813 copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest)); 00814 } 00815 else 00816 { 00817 // work directly on the destination array 00818 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels ); 00819 } 00820 } 00821 00822 template <class SrcIterator, class SrcShape, class SrcAccessor, 00823 class DestIterator, class DestAccessor, class KernelIterator> 00824 inline 00825 void separableConvolveMultiArray( 00826 triple<SrcIterator, SrcShape, SrcAccessor> const & source, 00827 pair<DestIterator, DestAccessor> const & dest, 00828 KernelIterator kit, 00829 SrcShape const & start = SrcShape(), 00830 SrcShape const & stop = SrcShape()) 00831 { 00832 separableConvolveMultiArray( source.first, source.second, source.third, 00833 dest.first, dest.second, kit, start, stop ); 00834 } 00835 00836 template <class SrcIterator, class SrcShape, class SrcAccessor, 00837 class DestIterator, class DestAccessor, class T> 00838 inline void 00839 separableConvolveMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src, 00840 DestIterator d, DestAccessor dest, 00841 Kernel1D<T> const & kernel, 00842 SrcShape const & start = SrcShape(), 00843 SrcShape const & stop = SrcShape()) 00844 { 00845 ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel); 00846 00847 separableConvolveMultiArray( s, shape, src, d, dest, kernels.begin(), start, stop); 00848 } 00849 00850 template <class SrcIterator, class SrcShape, class SrcAccessor, 00851 class DestIterator, class DestAccessor, class T> 00852 inline void 00853 separableConvolveMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 00854 pair<DestIterator, DestAccessor> const & dest, 00855 Kernel1D<T> const & kernel, 00856 SrcShape const & start = SrcShape(), 00857 SrcShape const & stop = SrcShape()) 00858 { 00859 ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel); 00860 00861 separableConvolveMultiArray( source.first, source.second, source.third, 00862 dest.first, dest.second, kernels.begin(), start, stop); 00863 } 00864 00865 /********************************************************/ 00866 /* */ 00867 /* convolveMultiArrayOneDimension */ 00868 /* */ 00869 /********************************************************/ 00870 00871 /** \brief Convolution along a single dimension of a multi-dimensional arrays. 00872 00873 This function computes a convolution along one dimension (specified by 00874 the parameter <tt>dim</tt> of the given multi-dimensional array with the given 00875 <tt>kernel</tt>. Both source and destination arrays are represented by 00876 iterators, shape objects and accessors. The destination array is required to 00877 already have the correct size. 00878 00879 If <tt>start</tt> and <tt>stop</tt> have non-default values, they must represent 00880 a valid subarray of the input array. The convolution is then restricted to that 00881 subarray, and it is assumed that the output array only refers to the 00882 subarray (i.e. <tt>diter</tt> points to the element corresponding to 00883 <tt>start</tt>). 00884 00885 This function may work in-place, which means that <tt>siter == diter</tt> is allowed. 00886 00887 <b> Declarations:</b> 00888 00889 pass arguments explicitly: 00890 \code 00891 namespace vigra { 00892 template <class SrcIterator, class SrcShape, class SrcAccessor, 00893 class DestIterator, class DestAccessor, class T> 00894 void 00895 convolveMultiArrayOneDimension(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 00896 DestIterator diter, DestAccessor dest, 00897 unsigned int dim, vigra::Kernel1D<T> const & kernel, 00898 SrcShape const & start = SrcShape(), 00899 SrcShape const & stop = SrcShape()); 00900 } 00901 \endcode 00902 00903 use argument objects in conjunction with \ref ArgumentObjectFactories : 00904 \code 00905 namespace vigra { 00906 template <class SrcIterator, class SrcShape, class SrcAccessor, 00907 class DestIterator, class DestAccessor, class T> 00908 void 00909 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 00910 pair<DestIterator, DestAccessor> const & dest, 00911 unsigned int dim, vigra::Kernel1D<T> const & kernel, 00912 SrcShape const & start = SrcShape(), 00913 SrcShape const & stop = SrcShape()); 00914 } 00915 \endcode 00916 00917 <b> Usage:</b> 00918 00919 <b>\#include</b> <vigra/multi_convolution.hxx> 00920 00921 \code 00922 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 00923 MultiArray<3, unsigned char> source(shape); 00924 MultiArray<3, float> dest(shape); 00925 ... 00926 Kernel1D<float> gauss; 00927 gauss.initGaussian(sigma); 00928 00929 // perform Gaussian smoothing along dimensions 1 (height) 00930 convolveMultiArrayOneDimension(srcMultiArrayRange(source), destMultiArray(dest), 1, gauss); 00931 \endcode 00932 00933 \see separableConvolveMultiArray() 00934 */ 00935 doxygen_overloaded_function(template <...> void convolveMultiArrayOneDimension) 00936 00937 template <class SrcIterator, class SrcShape, class SrcAccessor, 00938 class DestIterator, class DestAccessor, class T> 00939 void 00940 convolveMultiArrayOneDimension(SrcIterator s, SrcShape const & shape, SrcAccessor src, 00941 DestIterator d, DestAccessor dest, 00942 unsigned int dim, vigra::Kernel1D<T> const & kernel, 00943 SrcShape const & start = SrcShape(), 00944 SrcShape const & stop = SrcShape()) 00945 { 00946 enum { N = 1 + SrcIterator::level }; 00947 vigra_precondition( dim < N, 00948 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller " 00949 "than the data dimensionality" ); 00950 00951 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType; 00952 typedef typename AccessorTraits<TmpType>::default_const_accessor TmpAccessor; 00953 ArrayVector<TmpType> tmp( shape[dim] ); 00954 00955 typedef MultiArrayNavigator<SrcIterator, N> SNavigator; 00956 typedef MultiArrayNavigator<DestIterator, N> DNavigator; 00957 00958 SrcShape sstart, sstop(shape), dstart, dstop(shape); 00959 00960 if(stop != SrcShape()) 00961 { 00962 sstart = start; 00963 sstop = stop; 00964 sstart[dim] = 0; 00965 sstop[dim] = shape[dim]; 00966 dstop = stop - start; 00967 } 00968 00969 SNavigator snav( s, sstart, sstop, dim ); 00970 DNavigator dnav( d, dstart, dstop, dim ); 00971 00972 for( ; snav.hasMore(); snav++, dnav++ ) 00973 { 00974 // first copy source to temp for maximum cache efficiency 00975 copyLine( snav.begin(), snav.end(), src, 00976 tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() ); 00977 00978 convolveLine(srcIterRange( tmp.begin(), tmp.end(), TmpAccessor()), 00979 destIter( dnav.begin(), dest ), 00980 kernel1d( kernel), start[dim], stop[dim]); 00981 } 00982 } 00983 00984 template <class SrcIterator, class SrcShape, class SrcAccessor, 00985 class DestIterator, class DestAccessor, class T> 00986 inline void 00987 convolveMultiArrayOneDimension(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 00988 pair<DestIterator, DestAccessor> const & dest, 00989 unsigned int dim, vigra::Kernel1D<T> const & kernel, 00990 SrcShape const & start = SrcShape(), 00991 SrcShape const & stop = SrcShape()) 00992 { 00993 convolveMultiArrayOneDimension(source.first, source.second, source.third, 00994 dest.first, dest.second, dim, kernel, start, stop); 00995 } 00996 00997 /********************************************************/ 00998 /* */ 00999 /* gaussianSmoothMultiArray */ 01000 /* */ 01001 /********************************************************/ 01002 01003 /** \brief Isotropic Gaussian smoothing of a multi-dimensional arrays. 01004 01005 This function computes an isotropic convolution of the given N-dimensional 01006 array with a Gaussian filter at the given standard deviation <tt>sigma</tt>. 01007 Both source and destination arrays are represented by 01008 iterators, shape objects and accessors. The destination array is required to 01009 already have the correct size. This function may work in-place, which means 01010 that <tt>siter == diter</tt> is allowed. It is implemented by a call to 01011 \ref separableConvolveMultiArray() with the appropriate kernel. 01012 01013 Anisotropic data should be passed with appropriate 01014 \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional 01015 unless the parameter <tt>sigma</tt> is left out. 01016 01017 <b> Declarations:</b> 01018 01019 pass arguments explicitly: 01020 \code 01021 namespace vigra { 01022 template <class SrcIterator, class SrcShape, class SrcAccessor, 01023 class DestIterator, class DestAccessor> 01024 void 01025 gaussianSmoothMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 01026 DestIterator diter, DestAccessor dest, 01027 double sigma, const ConvolutionOptions<N> & opt); 01028 } 01029 \endcode 01030 01031 use argument objects in conjunction with \ref ArgumentObjectFactories : 01032 \code 01033 namespace vigra { 01034 template <class SrcIterator, class SrcShape, class SrcAccessor, 01035 class DestIterator, class DestAccessor> 01036 void 01037 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01038 pair<DestIterator, DestAccessor> const & dest, 01039 double sigma, const ConvolutionOptions<N> & opt); 01040 } 01041 \endcode 01042 01043 <b> Usage:</b> 01044 01045 <b>\#include</b> <vigra/multi_convolution.hxx> 01046 01047 \code 01048 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 01049 MultiArray<3, unsigned char> source(shape); 01050 MultiArray<3, float> dest(shape); 01051 ... 01052 // perform isotropic Gaussian smoothing at scale 'sigma' 01053 gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma); 01054 \endcode 01055 01056 <b> Usage with anisotropic data:</b> 01057 01058 <b>\#include</b> <vigra/multi_convolution.hxx> 01059 01060 \code 01061 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 01062 MultiArray<3, unsigned char> source(shape); 01063 MultiArray<3, float> dest(shape); 01064 TinyVector<float, 3> step_size; 01065 TinyVector<float, 3> resolution_sigmas; 01066 ... 01067 // perform anisotropic Gaussian smoothing at scale 'sigma' 01068 gaussianSmoothMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma, 01069 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas)); 01070 \endcode 01071 01072 \see separableConvolveMultiArray() 01073 */ 01074 doxygen_overloaded_function(template <...> void gaussianSmoothMultiArray) 01075 01076 template <class SrcIterator, class SrcShape, class SrcAccessor, 01077 class DestIterator, class DestAccessor> 01078 void 01079 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src, 01080 DestIterator d, DestAccessor dest, 01081 const ConvolutionOptions<SrcShape::static_size> & opt, 01082 const char *const function_name = "gaussianSmoothMultiArray" ) 01083 { 01084 typedef typename DestAccessor::value_type DestType; 01085 01086 static const int N = SrcShape::static_size; 01087 01088 typename ConvolutionOptions<N>::ScaleIterator params = opt.scaleParams(); 01089 ArrayVector<Kernel1D<double> > kernels(N); 01090 01091 for (int dim = 0; dim < N; ++dim, ++params) 01092 kernels[dim].initGaussian(params.sigma_scaled(function_name), 1.0, opt.window_ratio); 01093 01094 separableConvolveMultiArray(s, shape, src, d, dest, kernels.begin(), opt.from_point, opt.to_point); 01095 } 01096 01097 template <class SrcIterator, class SrcShape, class SrcAccessor, 01098 class DestIterator, class DestAccessor> 01099 inline void 01100 gaussianSmoothMultiArray( SrcIterator s, SrcShape const & shape, SrcAccessor src, 01101 DestIterator d, DestAccessor dest, double sigma, 01102 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01103 { 01104 ConvolutionOptions<SrcShape::static_size> par = opt; 01105 gaussianSmoothMultiArray(s, shape, src, d, dest, par.stdDev(sigma)); 01106 } 01107 01108 template <class SrcIterator, class SrcShape, class SrcAccessor, 01109 class DestIterator, class DestAccessor> 01110 inline void 01111 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01112 pair<DestIterator, DestAccessor> const & dest, 01113 const ConvolutionOptions<SrcShape::static_size> & opt) 01114 { 01115 gaussianSmoothMultiArray( source.first, source.second, source.third, 01116 dest.first, dest.second, opt ); 01117 } 01118 01119 template <class SrcIterator, class SrcShape, class SrcAccessor, 01120 class DestIterator, class DestAccessor> 01121 inline void 01122 gaussianSmoothMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01123 pair<DestIterator, DestAccessor> const & dest, double sigma, 01124 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01125 { 01126 gaussianSmoothMultiArray( source.first, source.second, source.third, 01127 dest.first, dest.second, sigma, opt ); 01128 } 01129 01130 01131 /********************************************************/ 01132 /* */ 01133 /* gaussianGradientMultiArray */ 01134 /* */ 01135 /********************************************************/ 01136 01137 /** \brief Calculate Gaussian gradient of a multi-dimensional arrays. 01138 01139 This function computes the Gaussian gradient of the given N-dimensional 01140 array with a sequence of first-derivative-of-Gaussian filters at the given 01141 standard deviation <tt>sigma</tt> (differentiation is applied to each dimension 01142 in turn, starting with the innermost dimension). Both source and destination arrays 01143 are represented by iterators, shape objects and accessors. The destination array is 01144 required to have a vector valued pixel type with as many elements as the number of 01145 dimensions. This function is implemented by calls to 01146 \ref separableConvolveMultiArray() with the appropriate kernels. 01147 01148 Anisotropic data should be passed with appropriate 01149 \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional 01150 unless the parameter <tt>sigma</tt> is left out. 01151 01152 <b> Declarations:</b> 01153 01154 pass arguments explicitly: 01155 \code 01156 namespace vigra { 01157 template <class SrcIterator, class SrcShape, class SrcAccessor, 01158 class DestIterator, class DestAccessor> 01159 void 01160 gaussianGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 01161 DestIterator diter, DestAccessor dest, 01162 double sigma, const ConvolutionOptions<N> & opt); 01163 } 01164 \endcode 01165 01166 use argument objects in conjunction with \ref ArgumentObjectFactories : 01167 \code 01168 namespace vigra { 01169 template <class SrcIterator, class SrcShape, class SrcAccessor, 01170 class DestIterator, class DestAccessor> 01171 void 01172 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01173 pair<DestIterator, DestAccessor> const & dest, 01174 double sigma, const ConvolutionOptions<N> & opt); 01175 } 01176 \endcode 01177 01178 <b> Usage:</b> 01179 01180 <b>\#include</b> <vigra/multi_convolution.hxx> 01181 01182 \code 01183 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 01184 MultiArray<3, unsigned char> source(shape); 01185 MultiArray<3, TinyVector<float, 3> > dest(shape); 01186 ... 01187 // compute Gaussian gradient at scale sigma 01188 gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma); 01189 \endcode 01190 01191 <b> Usage with anisotropic data:</b> 01192 01193 <b>\#include</b> <vigra/multi_convolution.hxx> 01194 01195 \code 01196 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 01197 MultiArray<3, unsigned char> source(shape); 01198 MultiArray<3, TinyVector<float, 3> > dest(shape); 01199 TinyVector<float, 3> step_size; 01200 TinyVector<float, 3> resolution_sigmas; 01201 ... 01202 // compute Gaussian gradient at scale sigma 01203 gaussianGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma, 01204 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas)); 01205 \endcode 01206 01207 <b> Required Interface:</b> 01208 01209 see \ref separableConvolveMultiArray(), in addition: 01210 01211 \code 01212 int dimension = 0; 01213 VectorElementAccessor<DestAccessor> elementAccessor(0, dest); 01214 \endcode 01215 01216 \see separableConvolveMultiArray() 01217 */ 01218 doxygen_overloaded_function(template <...> void gaussianGradientMultiArray) 01219 01220 template <class SrcIterator, class SrcShape, class SrcAccessor, 01221 class DestIterator, class DestAccessor> 01222 void 01223 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01224 DestIterator di, DestAccessor dest, 01225 ConvolutionOptions<SrcShape::static_size> const & opt, 01226 const char *const function_name = "gaussianGradientMultiArray") 01227 { 01228 typedef typename DestAccessor::value_type DestType; 01229 typedef typename DestType::value_type DestValueType; 01230 typedef typename NumericTraits<DestValueType>::RealPromote KernelType; 01231 01232 static const int N = SrcShape::static_size; 01233 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType; 01234 01235 for(int k=0; k<N; ++k) 01236 if(shape[k] <=0) 01237 return; 01238 01239 vigra_precondition(N == (int)dest.size(di), 01240 "gaussianGradientMultiArray(): Wrong number of channels in output array."); 01241 01242 ParamType params = opt.scaleParams(); 01243 ParamType params2(params); 01244 01245 ArrayVector<Kernel1D<KernelType> > plain_kernels(N); 01246 for (int dim = 0; dim < N; ++dim, ++params) 01247 { 01248 double sigma = params.sigma_scaled(function_name); 01249 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio); 01250 } 01251 01252 typedef VectorElementAccessor<DestAccessor> ElementAccessor; 01253 01254 // compute gradient components 01255 for (int dim = 0; dim < N; ++dim, ++params2) 01256 { 01257 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels); 01258 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 1, 1.0, opt.window_ratio); 01259 detail::scaleKernel(kernels[dim], 1.0 / params2.step_size()); 01260 separableConvolveMultiArray(si, shape, src, di, ElementAccessor(dim, dest), kernels.begin(), 01261 opt.from_point, opt.to_point); 01262 } 01263 } 01264 01265 template <class SrcIterator, class SrcShape, class SrcAccessor, 01266 class DestIterator, class DestAccessor> 01267 void 01268 gaussianGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01269 DestIterator di, DestAccessor dest, double sigma, 01270 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01271 { 01272 ConvolutionOptions<SrcShape::static_size> par = opt; 01273 gaussianGradientMultiArray(si, shape, src, di, dest, par.stdDev(sigma)); 01274 } 01275 01276 template <class SrcIterator, class SrcShape, class SrcAccessor, 01277 class DestIterator, class DestAccessor> 01278 inline void 01279 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01280 pair<DestIterator, DestAccessor> const & dest, 01281 ConvolutionOptions<SrcShape::static_size> const & opt ) 01282 { 01283 gaussianGradientMultiArray( source.first, source.second, source.third, 01284 dest.first, dest.second, opt ); 01285 } 01286 01287 template <class SrcIterator, class SrcShape, class SrcAccessor, 01288 class DestIterator, class DestAccessor> 01289 inline void 01290 gaussianGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01291 pair<DestIterator, DestAccessor> const & dest, 01292 double sigma, 01293 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01294 { 01295 gaussianGradientMultiArray( source.first, source.second, source.third, 01296 dest.first, dest.second, sigma, opt ); 01297 } 01298 01299 /********************************************************/ 01300 /* */ 01301 /* symmetricGradientMultiArray */ 01302 /* */ 01303 /********************************************************/ 01304 01305 /** \brief Calculate gradient of a multi-dimensional arrays using symmetric difference filters. 01306 01307 This function computes the gradient of the given N-dimensional 01308 array with a sequence of symmetric difference filters a (differentiation is applied 01309 to each dimension in turn, starting with the innermost dimension). Both source and 01310 destination arrays are represented by iterators, shape objects and accessors. 01311 The destination array is required to have a vector valued pixel type with as many 01312 elements as the number of dimensions. This function is implemented by calls to 01313 \ref convolveMultiArrayOneDimension() with the symmetric difference kernel. 01314 01315 Anisotropic data should be passed with appropriate 01316 \ref ConvolutionOptions, the parameter <tt>opt</tt> is optional 01317 otherwise. 01318 01319 <b> Declarations:</b> 01320 01321 pass arguments explicitly: 01322 \code 01323 namespace vigra { 01324 template <class SrcIterator, class SrcShape, class SrcAccessor, 01325 class DestIterator, class DestAccessor> 01326 void 01327 symmetricGradientMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 01328 DestIterator diter, DestAccessor dest, 01329 const ConvolutionOptions<N> & opt); 01330 } 01331 \endcode 01332 01333 use argument objects in conjunction with \ref ArgumentObjectFactories : 01334 \code 01335 namespace vigra { 01336 template <class SrcIterator, class SrcShape, class SrcAccessor, 01337 class DestIterator, class DestAccessor> 01338 void 01339 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01340 pair<DestIterator, DestAccessor> const & dest, 01341 const ConvolutionOptions<N> & opt); 01342 } 01343 \endcode 01344 01345 <b> Usage:</b> 01346 01347 <b>\#include</b> <vigra/multi_convolution.hxx> 01348 01349 \code 01350 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 01351 MultiArray<3, unsigned char> source(shape); 01352 MultiArray<3, TinyVector<float, 3> > dest(shape); 01353 ... 01354 // compute gradient 01355 symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest)); 01356 \endcode 01357 01358 <b> Usage with anisotropic data:</b> 01359 01360 <b>\#include</b> <vigra/multi_convolution.hxx> 01361 01362 \code 01363 MultiArray<3, unsigned char>::size_type shape(width, height, depth); 01364 MultiArray<3, unsigned char> source(shape); 01365 MultiArray<3, TinyVector<float, 3> > dest(shape); 01366 TinyVector<float, 3> step_size; 01367 ... 01368 // compute gradient 01369 symmetricGradientMultiArray(srcMultiArrayRange(source), destMultiArray(dest), 01370 ConvolutionOptions<3>().stepSize(step_size)); 01371 \endcode 01372 01373 <b> Required Interface:</b> 01374 01375 see \ref convolveMultiArrayOneDimension(), in addition: 01376 01377 \code 01378 int dimension = 0; 01379 VectorElementAccessor<DestAccessor> elementAccessor(0, dest); 01380 \endcode 01381 01382 \see convolveMultiArrayOneDimension() 01383 */ 01384 doxygen_overloaded_function(template <...> void symmetricGradientMultiArray) 01385 01386 template <class SrcIterator, class SrcShape, class SrcAccessor, 01387 class DestIterator, class DestAccessor> 01388 void 01389 symmetricGradientMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01390 DestIterator di, DestAccessor dest, 01391 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01392 { 01393 typedef typename DestAccessor::value_type DestType; 01394 typedef typename DestType::value_type DestValueType; 01395 typedef typename NumericTraits<DestValueType>::RealPromote KernelType; 01396 01397 static const int N = SrcShape::static_size; 01398 typedef typename ConvolutionOptions<N>::StepIterator StepType; 01399 01400 for(int k=0; k<N; ++k) 01401 if(shape[k] <=0) 01402 return; 01403 01404 vigra_precondition(N == (int)dest.size(di), 01405 "symmetricGradientMultiArray(): Wrong number of channels in output array."); 01406 01407 Kernel1D<KernelType> filter; 01408 filter.initSymmetricGradient(); 01409 01410 StepType step_size_it = opt.stepParams(); 01411 01412 typedef VectorElementAccessor<DestAccessor> ElementAccessor; 01413 01414 // compute gradient components 01415 for (int d = 0; d < N; ++d, ++step_size_it) 01416 { 01417 Kernel1D<KernelType> symmetric(filter); 01418 detail::scaleKernel(symmetric, 1 / *step_size_it); 01419 convolveMultiArrayOneDimension(si, shape, src, 01420 di, ElementAccessor(d, dest), 01421 d, symmetric, opt.from_point, opt.to_point); 01422 } 01423 } 01424 01425 template <class SrcIterator, class SrcShape, class SrcAccessor, 01426 class DestIterator, class DestAccessor> 01427 inline void 01428 symmetricGradientMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01429 pair<DestIterator, DestAccessor> const & dest, 01430 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01431 { 01432 symmetricGradientMultiArray(source.first, source.second, source.third, 01433 dest.first, dest.second, opt); 01434 } 01435 01436 /********************************************************/ 01437 /* */ 01438 /* laplacianOfGaussianMultiArray */ 01439 /* */ 01440 /********************************************************/ 01441 01442 /** \brief Calculate Laplacian of a N-dimensional arrays using Gaussian derivative filters. 01443 01444 This function computes the Laplacian the given N-dimensional 01445 array with a sequence of second-derivative-of-Gaussian filters at the given 01446 standard deviation <tt>sigma</tt>. Both source and destination arrays 01447 are represented by iterators, shape objects and accessors. Both source and destination 01448 arrays must have scalar value_type. This function is implemented by calls to 01449 \ref separableConvolveMultiArray() with the appropriate kernels, followed by summation. 01450 01451 Anisotropic data should be passed with appropriate 01452 \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional 01453 unless the parameter <tt>sigma</tt> is left out. 01454 01455 <b> Declarations:</b> 01456 01457 pass arguments explicitly: 01458 \code 01459 namespace vigra { 01460 template <class SrcIterator, class SrcShape, class SrcAccessor, 01461 class DestIterator, class DestAccessor> 01462 void 01463 laplacianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 01464 DestIterator diter, DestAccessor dest, 01465 double sigma, const ConvolutionOptions<N> & opt); 01466 } 01467 \endcode 01468 01469 use argument objects in conjunction with \ref ArgumentObjectFactories : 01470 \code 01471 namespace vigra { 01472 template <class SrcIterator, class SrcShape, class SrcAccessor, 01473 class DestIterator, class DestAccessor> 01474 void 01475 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01476 pair<DestIterator, DestAccessor> const & dest, 01477 double sigma, const ConvolutionOptions<N> & opt); 01478 } 01479 \endcode 01480 01481 <b> Usage:</b> 01482 01483 <b>\#include</b> <vigra/multi_convolution.hxx> 01484 01485 \code 01486 MultiArray<3, float> source(shape); 01487 MultiArray<3, float> laplacian(shape); 01488 ... 01489 // compute Laplacian at scale sigma 01490 laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma); 01491 \endcode 01492 01493 <b> Usage with anisotropic data:</b> 01494 01495 <b>\#include</b> <vigra/multi_convolution.hxx> 01496 01497 \code 01498 MultiArray<3, float> source(shape); 01499 MultiArray<3, float> laplacian(shape); 01500 TinyVector<float, 3> step_size; 01501 TinyVector<float, 3> resolution_sigmas; 01502 ... 01503 // compute Laplacian at scale sigma 01504 laplacianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(laplacian), sigma, 01505 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas)); 01506 \endcode 01507 01508 <b> Required Interface:</b> 01509 01510 see \ref separableConvolveMultiArray(), in addition: 01511 01512 \code 01513 int dimension = 0; 01514 VectorElementAccessor<DestAccessor> elementAccessor(0, dest); 01515 \endcode 01516 01517 \see separableConvolveMultiArray() 01518 */ 01519 doxygen_overloaded_function(template <...> void laplacianOfGaussianMultiArray) 01520 01521 template <class SrcIterator, class SrcShape, class SrcAccessor, 01522 class DestIterator, class DestAccessor> 01523 void 01524 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01525 DestIterator di, DestAccessor dest, 01526 ConvolutionOptions<SrcShape::static_size> const & opt ) 01527 { 01528 using namespace functor; 01529 01530 typedef typename DestAccessor::value_type DestType; 01531 typedef typename NumericTraits<DestType>::RealPromote KernelType; 01532 typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor; 01533 01534 static const int N = SrcShape::static_size; 01535 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType; 01536 01537 ParamType params = opt.scaleParams(); 01538 ParamType params2(params); 01539 01540 ArrayVector<Kernel1D<KernelType> > plain_kernels(N); 01541 for (int dim = 0; dim < N; ++dim, ++params) 01542 { 01543 double sigma = params.sigma_scaled("laplacianOfGaussianMultiArray"); 01544 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio); 01545 } 01546 01547 SrcShape dshape(shape); 01548 if(opt.to_point != SrcShape()) 01549 dshape = opt.to_point - opt.from_point; 01550 01551 MultiArray<N, KernelType> derivative(dshape); 01552 01553 // compute 2nd derivatives and sum them up 01554 for (int dim = 0; dim < N; ++dim, ++params2) 01555 { 01556 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels); 01557 kernels[dim].initGaussianDerivative(params2.sigma_scaled(), 2, 1.0, opt.window_ratio); 01558 detail::scaleKernel(kernels[dim], 1.0 / sq(params2.step_size())); 01559 01560 if (dim == 0) 01561 { 01562 separableConvolveMultiArray( si, shape, src, 01563 di, dest, kernels.begin(), opt.from_point, opt.to_point); 01564 } 01565 else 01566 { 01567 separableConvolveMultiArray( si, shape, src, 01568 derivative.traverser_begin(), DerivativeAccessor(), 01569 kernels.begin(), opt.from_point, opt.to_point); 01570 combineTwoMultiArrays(di, dshape, dest, derivative.traverser_begin(), DerivativeAccessor(), 01571 di, dest, Arg1() + Arg2() ); 01572 } 01573 } 01574 } 01575 01576 template <class SrcIterator, class SrcShape, class SrcAccessor, 01577 class DestIterator, class DestAccessor> 01578 void 01579 laplacianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01580 DestIterator di, DestAccessor dest, double sigma, 01581 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01582 { 01583 ConvolutionOptions<SrcShape::static_size> par = opt; 01584 laplacianOfGaussianMultiArray(si, shape, src, di, dest, par.stdDev(sigma)); 01585 } 01586 01587 template <class SrcIterator, class SrcShape, class SrcAccessor, 01588 class DestIterator, class DestAccessor> 01589 inline void 01590 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01591 pair<DestIterator, DestAccessor> const & dest, 01592 ConvolutionOptions<SrcShape::static_size> const & opt ) 01593 { 01594 laplacianOfGaussianMultiArray( source.first, source.second, source.third, 01595 dest.first, dest.second, opt ); 01596 } 01597 01598 template <class SrcIterator, class SrcShape, class SrcAccessor, 01599 class DestIterator, class DestAccessor> 01600 inline void 01601 laplacianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01602 pair<DestIterator, DestAccessor> const & dest, 01603 double sigma, 01604 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01605 { 01606 laplacianOfGaussianMultiArray( source.first, source.second, source.third, 01607 dest.first, dest.second, sigma, opt ); 01608 } 01609 01610 /********************************************************/ 01611 /* */ 01612 /* hessianOfGaussianMultiArray */ 01613 /* */ 01614 /********************************************************/ 01615 01616 /** \brief Calculate Hessian matrix of a N-dimensional arrays using Gaussian derivative filters. 01617 01618 This function computes the Hessian matrix the given scalar N-dimensional 01619 array with a sequence of second-derivative-of-Gaussian filters at the given 01620 standard deviation <tt>sigma</tt>. Both source and destination arrays 01621 are represented by iterators, shape objects and accessors. The destination array must 01622 have a vector valued element type with N*(N+1)/2 elements (it represents the 01623 upper triangular part of the symmetric Hessian matrix). This function is implemented by calls to 01624 \ref separableConvolveMultiArray() with the appropriate kernels. 01625 01626 Anisotropic data should be passed with appropriate 01627 \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional 01628 unless the parameter <tt>sigma</tt> is left out. 01629 01630 <b> Declarations:</b> 01631 01632 pass arguments explicitly: 01633 \code 01634 namespace vigra { 01635 template <class SrcIterator, class SrcShape, class SrcAccessor, 01636 class DestIterator, class DestAccessor> 01637 void 01638 hessianOfGaussianMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 01639 DestIterator diter, DestAccessor dest, 01640 double sigma, const ConvolutionOptions<N> & opt); 01641 } 01642 \endcode 01643 01644 use argument objects in conjunction with \ref ArgumentObjectFactories : 01645 \code 01646 namespace vigra { 01647 template <class SrcIterator, class SrcShape, class SrcAccessor, 01648 class DestIterator, class DestAccessor> 01649 void 01650 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01651 pair<DestIterator, DestAccessor> const & dest, 01652 double sigma, const ConvolutionOptions<N> & opt); 01653 } 01654 \endcode 01655 01656 <b> Usage:</b> 01657 01658 <b>\#include</b> <vigra/multi_convolution.hxx> 01659 01660 \code 01661 MultiArray<3, float> source(shape); 01662 MultiArray<3, TinyVector<float, 6> > dest(shape); 01663 ... 01664 // compute Hessian at scale sigma 01665 hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma); 01666 \endcode 01667 01668 <b> Usage with anisotropic data:</b> 01669 01670 <b>\#include</b> <vigra/multi_convolution.hxx> 01671 01672 \code 01673 MultiArray<3, float> source(shape); 01674 MultiArray<3, TinyVector<float, 6> > dest(shape); 01675 TinyVector<float, 3> step_size; 01676 TinyVector<float, 3> resolution_sigmas; 01677 ... 01678 // compute Hessian at scale sigma 01679 hessianOfGaussianMultiArray(srcMultiArrayRange(source), destMultiArray(dest), sigma, 01680 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas)); 01681 \endcode 01682 01683 <b> Required Interface:</b> 01684 01685 see \ref separableConvolveMultiArray(), in addition: 01686 01687 \code 01688 int dimension = 0; 01689 VectorElementAccessor<DestAccessor> elementAccessor(0, dest); 01690 \endcode 01691 01692 \see separableConvolveMultiArray(), vectorToTensorMultiArray() 01693 */ 01694 doxygen_overloaded_function(template <...> void hessianOfGaussianMultiArray) 01695 01696 template <class SrcIterator, class SrcShape, class SrcAccessor, 01697 class DestIterator, class DestAccessor> 01698 void 01699 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01700 DestIterator di, DestAccessor dest, 01701 ConvolutionOptions<SrcShape::static_size> const & opt ) 01702 { 01703 typedef typename DestAccessor::value_type DestType; 01704 typedef typename DestType::value_type DestValueType; 01705 typedef typename NumericTraits<DestValueType>::RealPromote KernelType; 01706 01707 static const int N = SrcShape::static_size; 01708 static const int M = N*(N+1)/2; 01709 typedef typename ConvolutionOptions<N>::ScaleIterator ParamType; 01710 01711 for(int k=0; k<N; ++k) 01712 if(shape[k] <=0) 01713 return; 01714 01715 vigra_precondition(M == (int)dest.size(di), 01716 "hessianOfGaussianMultiArray(): Wrong number of channels in output array."); 01717 01718 ParamType params_init = opt.scaleParams(); 01719 01720 ArrayVector<Kernel1D<KernelType> > plain_kernels(N); 01721 ParamType params(params_init); 01722 for (int dim = 0; dim < N; ++dim, ++params) 01723 { 01724 double sigma = params.sigma_scaled("hessianOfGaussianMultiArray"); 01725 plain_kernels[dim].initGaussian(sigma, 1.0, opt.window_ratio); 01726 } 01727 01728 typedef VectorElementAccessor<DestAccessor> ElementAccessor; 01729 01730 // compute elements of the Hessian matrix 01731 ParamType params_i(params_init); 01732 for (int b=0, i=0; i<N; ++i, ++params_i) 01733 { 01734 ParamType params_j(params_i); 01735 for (int j=i; j<N; ++j, ++b, ++params_j) 01736 { 01737 ArrayVector<Kernel1D<KernelType> > kernels(plain_kernels); 01738 if(i == j) 01739 { 01740 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 2, 1.0, opt.window_ratio); 01741 } 01742 else 01743 { 01744 kernels[i].initGaussianDerivative(params_i.sigma_scaled(), 1, 1.0, opt.window_ratio); 01745 kernels[j].initGaussianDerivative(params_j.sigma_scaled(), 1, 1.0, opt.window_ratio); 01746 } 01747 detail::scaleKernel(kernels[i], 1 / params_i.step_size()); 01748 detail::scaleKernel(kernels[j], 1 / params_j.step_size()); 01749 separableConvolveMultiArray(si, shape, src, di, ElementAccessor(b, dest), 01750 kernels.begin(), opt.from_point, opt.to_point); 01751 } 01752 } 01753 } 01754 01755 template <class SrcIterator, class SrcShape, class SrcAccessor, 01756 class DestIterator, class DestAccessor> 01757 inline void 01758 hessianOfGaussianMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01759 DestIterator di, DestAccessor dest, double sigma, 01760 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01761 { 01762 ConvolutionOptions<SrcShape::static_size> par = opt; 01763 hessianOfGaussianMultiArray(si, shape, src, di, dest, par.stdDev(sigma)); 01764 } 01765 01766 template <class SrcIterator, class SrcShape, class SrcAccessor, 01767 class DestIterator, class DestAccessor> 01768 inline void 01769 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01770 pair<DestIterator, DestAccessor> const & dest, 01771 ConvolutionOptions<SrcShape::static_size> const & opt ) 01772 { 01773 hessianOfGaussianMultiArray( source.first, source.second, source.third, 01774 dest.first, dest.second, opt ); 01775 } 01776 01777 template <class SrcIterator, class SrcShape, class SrcAccessor, 01778 class DestIterator, class DestAccessor> 01779 inline void 01780 hessianOfGaussianMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01781 pair<DestIterator, DestAccessor> const & dest, 01782 double sigma, 01783 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01784 { 01785 hessianOfGaussianMultiArray( source.first, source.second, source.third, 01786 dest.first, dest.second, sigma, opt ); 01787 } 01788 01789 namespace detail { 01790 01791 template<int N, class VectorType> 01792 struct StructurTensorFunctor 01793 { 01794 typedef VectorType result_type; 01795 typedef typename VectorType::value_type ValueType; 01796 01797 template <class T> 01798 VectorType operator()(T const & in) const 01799 { 01800 VectorType res; 01801 for(int b=0, i=0; i<N; ++i) 01802 { 01803 for(int j=i; j<N; ++j, ++b) 01804 { 01805 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]); 01806 } 01807 } 01808 return res; 01809 } 01810 }; 01811 01812 } // namespace detail 01813 01814 /********************************************************/ 01815 /* */ 01816 /* structureTensorMultiArray */ 01817 /* */ 01818 /********************************************************/ 01819 01820 /** \brief Calculate th structure tensor of a multi-dimensional arrays. 01821 01822 This function computes the gradient (outer product) tensor for each element 01823 of the given N-dimensional array with first-derivative-of-Gaussian filters at 01824 the given <tt>innerScale</tt>, followed by Gaussian smoothing at <tt>outerScale</tt>. 01825 Both source and destination arrays are represented by iterators, shape objects and 01826 accessors. The destination array must have a vector valued pixel type with 01827 N*(N+1)/2 elements (it represents the upper triangular part of the symmetric 01828 structure tensor matrix). If the source array is also vector valued, the 01829 resulting structure tensor is the sum of the individual tensors for each channel. 01830 This function is implemented by calls to 01831 \ref separableConvolveMultiArray() with the appropriate kernels. 01832 01833 Anisotropic data should be passed with appropriate 01834 \ref ConvolutionOptions, the parameter <tt>opt</tt> is otherwise optional 01835 unless the parameters <tt>innerScale</tt> and <tt>outerScale</tt> are 01836 both left out. 01837 01838 <b> Declarations:</b> 01839 01840 pass arguments explicitly: 01841 \code 01842 namespace vigra { 01843 template <class SrcIterator, class SrcShape, class SrcAccessor, 01844 class DestIterator, class DestAccessor> 01845 void 01846 structureTensorMultiArray(SrcIterator siter, SrcShape const & shape, SrcAccessor src, 01847 DestIterator diter, DestAccessor dest, 01848 double innerScale, double outerScale, 01849 const ConvolutionOptions<N> & opt); 01850 } 01851 \endcode 01852 01853 use argument objects in conjunction with \ref ArgumentObjectFactories : 01854 \code 01855 namespace vigra { 01856 template <class SrcIterator, class SrcShape, class SrcAccessor, 01857 class DestIterator, class DestAccessor> 01858 void 01859 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01860 pair<DestIterator, DestAccessor> const & dest, 01861 double innerScale, double outerScale, 01862 const ConvolutionOptions<N> & opt); 01863 } 01864 \endcode 01865 01866 <b> Usage:</b> 01867 01868 <b>\#include</b> <vigra/multi_convolution.hxx> 01869 01870 \code 01871 MultiArray<3, RGBValue<float> > source(shape); 01872 MultiArray<3, TinyVector<float, 6> > dest(shape); 01873 ... 01874 // compute structure tensor at scales innerScale and outerScale 01875 structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale); 01876 \endcode 01877 01878 <b> Usage with anisotropic data:</b> 01879 01880 <b>\#include</b> <vigra/multi_convolution.hxx> 01881 01882 \code 01883 MultiArray<3, RGBValue<float> > source(shape); 01884 MultiArray<3, TinyVector<float, 6> > dest(shape); 01885 TinyVector<float, 3> step_size; 01886 TinyVector<float, 3> resolution_sigmas; 01887 ... 01888 // compute structure tensor at scales innerScale and outerScale 01889 structureTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest), innerScale, outerScale, 01890 ConvolutionOptions<3>().stepSize(step_size).resolutionStdDev(resolution_sigmas)); 01891 \endcode 01892 01893 <b> Required Interface:</b> 01894 01895 see \ref separableConvolveMultiArray(), in addition: 01896 01897 \code 01898 int dimension = 0; 01899 VectorElementAccessor<DestAccessor> elementAccessor(0, dest); 01900 \endcode 01901 01902 \see separableConvolveMultiArray(), vectorToTensorMultiArray() 01903 */ 01904 doxygen_overloaded_function(template <...> void structureTensorMultiArray) 01905 01906 template <class SrcIterator, class SrcShape, class SrcAccessor, 01907 class DestIterator, class DestAccessor> 01908 void 01909 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01910 DestIterator di, DestAccessor dest, 01911 ConvolutionOptions<SrcShape::static_size> const & opt) 01912 { 01913 static const int N = SrcShape::static_size; 01914 static const int M = N*(N+1)/2; 01915 01916 typedef typename DestAccessor::value_type DestType; 01917 typedef typename DestType::value_type DestValueType; 01918 typedef typename NumericTraits<DestValueType>::RealPromote KernelType; 01919 typedef TinyVector<KernelType, N> GradientVector; 01920 typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor; 01921 typedef typename AccessorTraits<DestType>::default_accessor GradientTensorAccessor; 01922 01923 for(int k=0; k<N; ++k) 01924 if(shape[k] <=0) 01925 return; 01926 01927 vigra_precondition(M == (int)dest.size(di), 01928 "structureTensorMultiArray(): Wrong number of channels in output array."); 01929 01930 ConvolutionOptions<N> innerOptions = opt; 01931 ConvolutionOptions<N> outerOptions = opt.outerOptions(); 01932 typename ConvolutionOptions<N>::ScaleIterator params = outerOptions.scaleParams(); 01933 01934 SrcShape gradientShape(shape); 01935 if(opt.to_point != SrcShape()) 01936 { 01937 for(int k=0; k<N; ++k, ++params) 01938 { 01939 Kernel1D<double> gauss; 01940 gauss.initGaussian(params.sigma_scaled("structureTensorMultiArray"), 1.0, opt.window_ratio); 01941 int dilation = gauss.right(); 01942 innerOptions.from_point[k] = std::max<MultiArrayIndex>(0, opt.from_point[k] - dilation); 01943 innerOptions.to_point[k] = std::min<MultiArrayIndex>(shape[k], opt.to_point[k] + dilation); 01944 } 01945 outerOptions.from_point -= innerOptions.from_point; 01946 outerOptions.to_point -= innerOptions.from_point; 01947 gradientShape = innerOptions.to_point - innerOptions.from_point; 01948 } 01949 01950 MultiArray<N, GradientVector> gradient(gradientShape); 01951 MultiArray<N, DestType> gradientTensor(gradientShape); 01952 gaussianGradientMultiArray(si, shape, src, 01953 gradient.traverser_begin(), GradientAccessor(), 01954 innerOptions, 01955 "structureTensorMultiArray"); 01956 01957 transformMultiArray(gradient.traverser_begin(), gradientShape, GradientAccessor(), 01958 gradientTensor.traverser_begin(), GradientTensorAccessor(), 01959 detail::StructurTensorFunctor<N, DestType>()); 01960 01961 gaussianSmoothMultiArray(gradientTensor.traverser_begin(), gradientShape, GradientTensorAccessor(), 01962 di, dest, outerOptions, 01963 "structureTensorMultiArray"); 01964 } 01965 01966 template <class SrcIterator, class SrcShape, class SrcAccessor, 01967 class DestIterator, class DestAccessor> 01968 inline void 01969 structureTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src, 01970 DestIterator di, DestAccessor dest, 01971 double innerScale, double outerScale, 01972 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01973 { 01974 ConvolutionOptions<SrcShape::static_size> par = opt; 01975 structureTensorMultiArray(si, shape, src, di, dest, 01976 par.stdDev(innerScale).outerScale(outerScale)); 01977 } 01978 01979 template <class SrcIterator, class SrcShape, class SrcAccessor, 01980 class DestIterator, class DestAccessor> 01981 inline void 01982 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01983 pair<DestIterator, DestAccessor> const & dest, 01984 ConvolutionOptions<SrcShape::static_size> const & opt ) 01985 { 01986 structureTensorMultiArray( source.first, source.second, source.third, 01987 dest.first, dest.second, opt ); 01988 } 01989 01990 01991 template <class SrcIterator, class SrcShape, class SrcAccessor, 01992 class DestIterator, class DestAccessor> 01993 inline void 01994 structureTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> const & source, 01995 pair<DestIterator, DestAccessor> const & dest, 01996 double innerScale, double outerScale, 01997 const ConvolutionOptions<SrcShape::static_size> & opt = ConvolutionOptions<SrcShape::static_size>()) 01998 { 01999 structureTensorMultiArray( source.first, source.second, source.third, 02000 dest.first, dest.second, 02001 innerScale, outerScale, opt); 02002 } 02003 02004 //@} 02005 02006 } //-- namespace vigra 02007 02008 02009 #endif //-- VIGRA_MULTI_CONVOLUTION_H
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|