three - dimensional ( 3d ) computer reconstruction of a target volume 
 or of a surface is an important activity in modern biomedical imaging .
the 
 accurate anatomical reconstruction in trauma or for use in image - guided 
 intervention relies on mathematical imaging technology ; and this paper develops 
 the mathematical technique of stochastic function recovery   and illustrates its use for 
 noisy boundary reconstruction .
this is an alternative approach to the standard 
 polynomial - based methods that we see as an add - on or complement to other 
 techniques in use or being developed to improve upon the reconstruction of 
 noisy boundary data to provide enhanced biomedical visualization .
the ability to distinguish features related to 
 boundaries is intrinsic to technology of visualization .
for example , in mri 
 imaging , a range of specialized methods have been developed for extracting 
 information from signals so as to reconstruct images representing internal body 
 structures .
boundary 
 recovery techniques apply to complex surgical procedures as with 
 electroanatomical mapping that tracks the position of catheters inside the body 
 with sparse signals recorded from electrodes at the tip of the catheter . 

resulting surface maps must integrate real - time measurements with preoperative 
 mr or ct images , and account for mapping data errors in registration and error 
 due to patient movement .
in general , when signals are affected by noise , it 
 must be effectively removed in order to improve the visualization and compared 
 with other medical imaging modalities , ultrasound images suffer from speckle 
 noise that often makes for weak or incomplete boundaries .
our approach is to use stochastic 
 convolution - deconvolution operators [ 1 , 57 ] , which have useful statistical properties , to smooth 
 noisy surface data in a manner that does not obliterate detail , and which 
 effectively removes gaussian noise .
the motivation for the approach is that , 
 intrinsically , stochastic interpolation using probabilistic kernels for the 
 generating function of the row space of the linear operator 
 performs well at removing noise when used to 
 approximate data .
however , the difficulty in applying these methods directly is 
 that they bias the data to a mean of zero . in approximating one - dimensional data ,
however , in 
 approximating multidimensional data , this can cause an apparent shift in the 
 approximant when working in parametric coordinates .
this is because the 
 gaussian kernels used smooth positive values to their mean value , thus 
 potentially shifting the coordinates nearer the origin .
when we approximate in 
 one dimension , the data values are shifted to the mean , however , the 
 coordinates are not touched .
if we are using parametric data , the coordinates 
 are the data values ; and this means that the data can be shifted in ( x , y ) or ( x , y , z ) space .
the less data existing or the greater the smoothing , the more 
 will be this shift . for example , when data is taken from a circle centered 
 on the origin , the approximating curve is a circular curve centered on the 
 origin , but of smaller radius , and the greater the smoothing applied in the 
 approximation ( if the data is very noisy ) , the smaller the area circumscribed 
 by the approximating curve . as the number of sample points
however , for coarse surface data that is only smoothly 
 varying , this can cause difficulties .
one approach is to make use of stochastic interpolation . creating a dense noisy data set from the sparse noisy data using 
 interpolation
is followed by approximating this fine data set to recover the 
 smoothed surface curve , thereby mitigating the shifting of the mean . while 
 workable , this approach has several disadvantages , most notably that it 
 requires a more costly interpolation step .
it also requires the application of 
 the technique more than once , and in the second or subsequent applications , the 
 approximation must be done on a fine data set , meaning that many more points 
 require approximation , again incurring a larger computational cost .
the solution is to make use of the convolution - deconvolution properties of stochastic interpolation combining the 
 densification step with the approximation step .
however , we introduce an approximate means for doing the interpolation that 
 interpolates for smooth data , but which approximates for noisy data , thus avoiding 
 the costly need to construct the inverse operator needed for interpolation .
consider the task of sampling a known function f(u ) at points f(uk ) = vk with uk  [ 0 , 1 ] so as to determine its value at xj  [ 0 , 1 ] .
the stochastic interpolant   to the data { ( uk , vk ) } , k = 1 , 
, n is given by 
 ( 1)bmnann1v , 
 where v = ( v1 , v2 ,  , vn ) is the data vector , and where ann is a row stochastic matrix whose coefficients consist of the n  n values ajk . 

choosing the generator of the row space of ann to be the bernstein functions   ( named after bernstein as the derivation of this form that can be obtained from the bernstein polynomials ) , we 
 have 
 ( 2)ajk=12[erf(uk+1xjn)+erf(xjukn ) ] 
 with  > 0 on the partition wk = ( uk+uk+1)/2 , with w0 =  and wn =  yielding a stochastic matrix . setting xj
= uj generates the entries of ann , and setting xj to any set of m consecutive values in [ 0 , 1 ] generates the coefficients of bmn , that is , the coefficients of bmn are constructed in the same manner as for ann , except that the nodes xj at which bjk is evaluated may differ from the values at which the data are given . while any probability density function ( pdf ) can be used , appropriately replacing the mean and variance of the gaussian in ( 2 ) , a pdf based on the normal distribution is consistent with the problem of 
 filtering gaussian noise . in stochastic interpolation , with the coefficient of ann generated by ( 2 ) , we can interpret ann as the discrete deconvolution of the data yielding the preimage generated by annv .
this preimage is then convolved by bmn to yield an m - vector of values that interpolates the data v at the output coordinates xj , j = 1 , 
, m ; the output coordinates are those that were used to generate the coefficients of bmn .
it is for this reason that we have elected to represent the matrix a using another symbol b since it is desirable to emphasize its role in convolution . defining wk = ( uk + uk+1)/2 with w0 = u0/2 and wn = 1 + un/2 with  constant yields a coefficient structure
in which ann is a diagonal matrix times a symmetric toeplitz matrix : ann = dtnn .
inversion or solution of these matrices can be accomplished in (n ) operations , however in the cases of interest to us , this is not necessary . it is possible to do better than this using an approximate inverse in which the row space of ann  ann is generated directly . while evaluating bmnann is still (n ) , it is significantly faster than the toeplitz matrix inversion of ann to obtain ann .
the approximate inverse is an approximation precisely because annann  inn , that is , in applying stochastic interpolation to the data v , there is an error:(3)e0=vannanniv=(innannanni)v . 

thus the interpolant to the data can be expressed using successive correction to the 
 errors using 
 ( 4)v = annanniv+e0=annanniv+aannie0+e1  =annanniv+annannie0+annannie1  + +annanniep+ep+1.substituting for ek , from   k = 0 to p , gives 
 ( 5)v = annanni(k=0p(innannanni)k)v+ep+1 , 
 and truncating the sum gives a working formula in which vv even for larger values of  , so that the method nearly interpolates .
truncating at p = 2 and applying the formula to generate m output values instead of n gives 
 ( 6)v=bmnanni(inn+gnn+gnn2)v , 
 where gnn = ( inn  annann ) .
provided that  in ( 2 ) is small , the error in constructing ann using ann is small , however if larger values of  are used , then greater smoothing is applied to 
 the data and the use of ( 6 ) becomes necessary when  is larger than 0.05 .
however , it will become apparent that it is because of this greater smoothing that it is unnecessary to apply any corrections as shown in ( 6 ) , and thus the direct computation of bmnannv is found to be convenient and efficient , requiring only a single matrix multiply . in working with stochastic data recovery , it is obvious that bmnv using the bernstein functions mollifies the data vector v , and thus provides an approximation vector of length m to the initial vector v of length n. consider the evaluation of bmnannv where the generator of the row space of ann and bmn are given by ( 2 ) .
this interpolates the data provided that b , meaning that the variance of the gaussian pdf that is used to generate bmn , is the same as the variance a that is used to generate ann . 

if instead b > a , then the preimage annv will be oversmoothed when bmn is applied to the preimage , and the result will be approximation .
similarly , if b < a , then the preimage annv will not be smoothed sufficiently , and the data will be roughened , or more appropriately it will be deconvolved . with statistical errors present in multidimensional 
 data , interpolation in parametric coordinates may yield an extremely complex 
 curve , and the errors may cause the curve to wiggle excessively , often crossing 
 over on itself .
however , as noted 
 in section 1 , simply approximating the data so as to smooth these errors may 
 introduce translation errors in the approximating surface representation , 
 particularly when the number of data points is small , or the smoothing 
 specified by  is large . since the approximation is 
 convergent , a simple work - around can be achieved by densification of the data 
 by interpolation as this will minimize this shift on subsequent smoothing .
this 
 leads to a computationally efficient approach to surface recovery that avoids 
 translation errors while smoothing the noise based on evaluating bmnannv , 
 where m is the desired number of output points 
 representing the boundary , and n are the number of input data points and b > a . to demonstrate this form , note that the intended construction is to evaluate bpnannv , 
 where p is significantly greater than n with b = a , and then applying smoothing to this densified data by multiplying by cmp , 
 where c > a . in effect
, this interpolates the data and then smooths it by the application of the approximation cmp to the densified data .
thus a two - step 
 algorithm can be described as follows : 
 compute densified interpolant bpnannv to the data vector v;compute boundry approximant cmp(bpnannv ) . 
 
compute densified interpolant bpnannv to the data vector v ; compute boundry approximant cmp(bpnannv ) .
this algorithm is equivalent to the following : 
 compute densified boundry approximate interpolant bmnannv , to the data vector v ; 
 that is , there exist c and b such that applying cmp to bpnann is the same , or nearly the same , as applying bmn to ann .
the use of ann for a wide range of a instead of ann introduces some additional smoothing , allowing 
 for less smoothing to be used on the convolution step , that is , when applying bmn .
the construction of ann is not difficult , remarkably being given by the inverse of the generator of the row space of ann .
the elements of the approximate inverse are given by the reciprocals of the coefficients of ann .
it is for this reason that the direct inversion of ann , or solution of the system using efficient toeplitz solvers , is not needed .
the only exception may be when the data is free of noise and exact interpolation without any smoothing is desired .
compute densified boundry approximate interpolant bmnannv , to the data vector v ; the reconstruction of surface data is done using a parametric representation { si}i = 1 in which si = ( x(ti ) , y(ti ) ) for two - dimensional data , and si = ( x(ti ) , y(ti ) , z(ti ) ) for three - dimensional data , where 0  ti  1 .
the recovery of the data is done using ( 1 ) or its modification using an approximate inverse ann , applied successively to the data pairs { ( ti , xi ) } , { ( ti , yi ) } in two dimensions and to the data triplets { ( ti , xi ) } , { ( ti , yi ) } , { ( ti , zi ) } in three dimensions .
for example , in the case of interpolating two - dimensional data , we apply bmnannx and bmnanny to obtain the interpolant to the position vector x and the position vector y , or we apply bmnannx and bmnanny to obtain the approximate interpolants . in reconstructing a parametrically represented surface generated from image data from pixel values , for instance , an algorithm for boundary detection and sorting is necessary . in our analysis
, it is assumed that this is available , however the errors generated in parameterizations of the surface may not be entirely random , and thus systematic errors in surface 
 representation will also be introduced . for the purpose of assessing the 
 performance of the boundary recovery algorithms
, it will be assumed that the 
 errors are random gaussian with a mean of zero , with variable variances  , 
 generated using the random variable 
 ( 7)=2log(r1 )  cos  (2r2 ) , 
 where r1 and r2 are two random numbers uniformly distributed 
 in the interval [ 0 , 1 ] .
thus the parametrically ordered data { ( x1(ti ) ,  , xd(ti ) ) } , with d = 2 or 3 , 
 and i = 1 ,  , n , 
 are perturbed to yield the data sets { ( x1(ti ) + ( 1)i ,  ,
xd(ti ) + ( d)i)}. in applying stochastic data recovery to the problem of finding the shape of a parametrically defined boundary , the problem of closed curves needs to be addressed . in the presence of large amount of noise
, the two 
 endpoints of the parametrically defined curves may not match : while ideally ( x1(t1 ) ,  , xd(t1 ) ) = ( x1(tn ) , 
, xd(tn ) ) in the case of a close loop , in the presence 
 of errors this will not be the case .
finally , in implementing the algorithm , it was found that the dependence of ( 2 ) on n in the denominator made the choice of  dependent on n , as the boundary data density increased , the recovery using any given value of  produced increasingly rougher curves as n increased , and thus it was found convenient to evaluate ajk and bjk based on a constant value of . additionally , the algorithm was applied to all of the boundary data associated with a parametric data set to construct the boundary curve , rather than 
 decomposing the data into overlapping blocks .
the merits of using all of the 
 data are a slight improvement in accuracy , while the drawbacks are that the 
 cost of evaluating the algorithm increases as the block size increases , and 
 this should be born in mind in applying the algorithm to large complex 
 three - dimensional data sets .
the values for a and b depend on the amount of noise , as well as on the presumed smoothness of the boundary data . while this seems to present 
 difficult choices , it is less complicated than it appears , as the choice for a will usually be any sufficiently small value .
for example , setting a  10 allows for representing the data as nearly piecewise linear along the boundaries .
in contrast , the choice for b requires some evaluation as this determines the smoothness of the recovered boundary .
we begin by applying the process to the recovery of the boundary of a disk defined parametrically by 12 uniformly distributed 
 points with no random errors in the data as shown in figure 1 .
the figure 
 clearly illustrates the aforementioned difficulty of attempting to approximate 
 ( smooth parametric data ) . in examining figure 1(a ) , it is also important to 
 observe that the beginning and the end of the curve are joined by a straight 
 line in this example .
the reason is that the slope of the approximant to the 
 data ( x(t1 ) , y(t1 ) ) , ( x(t2 ) , y(t2 ) ) ,  is not the same as the slope to the data  , ( x(tn1 ) , y(tn1 ) ) , ( x(tn ) , y(tn ) ) , and thus in this example , the first approximated point ( r1 , s1 ) does not agree with the last approximated point ( rm , sm ) , and are joined graphically with a straight line to close the curve .
a solution to the mismatch is to overlap the curves during reconstruction , that is , at several points away from the last point , and 
 ending the construction several points away from the first point , then only using the curve from the first to the last point . in all of the studies
presented , there is no attempt to overlap the curves in order to magnify these 
 boundary affects , and to demonstrate that they are mostly negligible , as seen in figure 1(b ) , whenever the construction is done correctly . for large data sets 
 where it may be computationally advantageous to block the data , overlapping the 
 endpoints is readily accomplished .
it is important to realize that if the number of 
 points representing the disk were to be much more than 12 , then it would have 
 been difficult to visualize the contraction of the recovered boundary curve 
 since the approximation is convergent .
thus for a sufficiently large number of 
 data points , the difference between the approximating curve and the curve 
 itself becomes arbitrarily small .
recovering the boundary of a disk becomes more 
 difficult , as shown in figure 2 , particularly if the amount of noise is quite 
 large . in this example
, the gaussian noise for a disk of radius 5 is specified as  = 3 yielding an extensively scattered data set .
while the level of noise in this illustrative example is much higher than would 
 be expected in any realistic imaging situation , the example serves a two - fold 
 purpose : ( 1 ) on the one had it shows the robustness of the method at recovering 
 a reasonable representation of the surface from the data that is more consistent 
 with noise than data , and ( 2 ) it shows that the effects of errors in 
 constructing the parameterization are less of an issue than might be presumed . 
 in the figure ,
the connectivity of the boundary data is not ordered in  moving counterclockwise around the circle , and so it is quite likely that any minor errors in parameterization , for example , 
 using even the simplest unconstrained nearest neighbor search of the data , 
 would cause unrecoverable errors in the surface representation that is 
 recovered . as expected , reducing the noise to  = 1/2 significantly improves the recovery , even for 
 half as many points , as shown in figure 3 .   
  while both of these figures are typical , the symmetry 
 and the smoothness of the boundary of the disk make the recovery somewhat less 
 challenging than for a more complex geometrical shape .
thus the performance of 
 the algorithm is examined on a star - shaped region generated from the vertex 
 data { ( 8,65 ) , ( 72,65 ) , ( 92,11 ) , ( 112,65 ) , ( 174,65 ) , 
 ( 122,100 ) , ( 142,155 ) , ( 92,121 ) , ( 42,155 ) , ( 60,100 ) , ( 8,65)}. note that the recovery makes use of  that on both steps is the same as that used in recovering the boundary of the disk . the shape of the star in figure 5 has 
 become more wiggly using the smoothing parameters the same as in the case of 
 lesser noise .
the problem is a classical one : there is no means to discern the 
 shape of the figure from the noisy data , except to note that an acceptable 
 shape is determined by the smoothness of the boundary that is intrinsic to the 
 figure . in this case , changing the smoothing can accommodate this subjective 
 assessment , as illustrated in figure 6 .
note that the recovered boundary is 
 consistent with the curve recovered from the less noisy data set : compare figures 
 4 and 6 as shown in figure 7 .
the effects of doubling the smoothing by taking b to be twice as large are clearly 
 evident .
since the noise in both cases was generated using the same seed , it is 
 only the magnitude of the excursions away from the star 's boundary that change , 
 and hence the figures are directly comparable .
this perhaps most plainly 
 illustrates that interaction between noise and smoothing , as the two curves are 
 nearly identical .
even when the noise is doubled again to  = 16 , the shape of the recovered curve using  = 0.0004 is remarkably consistent , as shown in figure 
 8 . at this level of noise
there is some loss of resolution of the limbs , 
 however the recovered boundary is recognizable as being related to the two 
 boundaries obtained in figures 4 and 6 .
while this is artificial in that the 
 amount of noise in the data in real problems is not known , it does demonstrate 
 the robustness of the algorithm at consistently recovering the boundary data 
 irrespective of the added noise .
since the approximate interpolant fails to 
 interpolate when  is large or , more appropriately , fails to 
 interpolate rapidly varying data , quite some effort was expended in developing 
 mechanisms for being able to use large  for smoothing and yet still maintain some 
 fidelity with the boundary data while constructing the preimage . as it 
 developed , this iterative correction was not necessary : the algorithm performed 
 well even without the introduction of any corrections . in part , this is due to 
 the rather simple shapes , and the relatively large amounts of noise that were 
 examined . in the case when a surface is oscillating rapidly with 
 the noise much less than this surface oscillation ,
it is clearly necessary to 
 implement the algorithm using this additional correction . indeed , given 
 sufficiently rough surface data , it may be necessary to use stochastic 
 interpolation as otherwise some high - frequency surface details will be 
 oversmoothed by the approximate inverse construction .
examination of several attempts at domain boundary 
 reconstruction in two dimensions is encouraging 
 using the stochastic data recovery techniques .
in particular , qualitative 
 assessments demonstrate the viability of utilizing stochastic function recovery 
 methods for reconstructing parametrically defined edges .
since the approach is 
 intrinsically one dimensional , its extension to three dimensions is not 
 difficult , and thus can easily be implemented and tested on more realistic 
 problems . moving to three dimensions poses no additional cost other than that 
 more data has to be processed , that is , the algorithmic costs scale directly 
 with the number of lines being evaluated , and the number of points on each line 
 at which the algorithm is applied .
it should be noted that the proposed technique does 
 not solve all problems involving noisy surface reconstruction , however it does 
 provide an additional tool for analysis .
the gaussian - based kernels ( the 
 bernstein functions ) which were used to generate the row space of the matrices , 
 were remarkably effective at cancelling the noise even under the most extreme 
 conditions where the noise essentially obliterated the image shape .
of course , 
 this noise was gaussian and so it is only reasonable that the proposed approach 
 would work well in these circumstances . for other types of noise ,
the use of alternative 
 probability density functions for generating the mollifiers is easily 
 accomplished , and thus the method has substantial design flexibility , and these 
 options need to be explored in detail to ascertain their utility at cancelling 
 these other types of noise .
the computational advantages of the technique are that 
 it requires only two matrix multiplies of the data vector , and thus the 
 approach is relatively cost effective .
moreover , the method is easily 
 implemented in parallel , and further computational gains in efficiency can be 
 achieved by blocking the data since typically only a segment of the entire data 
 vector is needed to recover the data in any region .
if fixed block sizes can be 
 implemented , then the cost of the two matrix vector multiplies can be further 
 reduced as the matrix - matrix multiply needs to be done only once , and so the 
 algorithm reduces to a single - block matrix - vector multiply .
 determining the outline or boundary contour of a two - dimensional object , or the surface of a  
 three - dimensional object poses difficulties particularly when there is substantial measurement noise or uncertainty . by adapting the mathematical approach of stochastic function recovery to this task , 
 it is possible to obtain usable estimates for these boundaries , even in the presence of large amounts of noise . 
 the technique is applied to parametric boundary data and has potential applications in biomedical imaging . 
 it should be considered as one of several techniques to improve the visualization of images . 
1. INTRODUCTION
2. DEVELOPMENT
3. RESULTS AND DISCUSSION
4. CONCLUSIONS