worldwide , it is estimated that 17.5 million people die of cardiovascular disease ( cvd ) every year , with an approximate cost of  310.23 billion ( world heart federation . available at http://www.worldheart.org . accessed 20090514 ) . because cvd remains the main cause of death in the world , considerable amounts of resources are devoted to cardiovascular research every year . the study of cardiocyte contractility has helped unveil the fundamental processes underlying heart function in health and disease [ 1 , 2 ] . the analysis of cardiocyte mechanics has historically proven an excellent tool in providing relevant information on the excitation - contraction coupling of the heart . many inotropic factors modulate the contractile behavior of the heart , which can be conveniently studied in enzymatically dissociated ( isolated ) cardiocytes [ 25 ] . researchers commonly measure calcium transient signals , gene and protein expression , and contractility to assess the function and state of these isolated cardiocytes in all their stages [ 1 , 6 ] . isolated adult , neonatal , and embryonic cardiocytes from mammalian hearts are widely used in cardiovascular research [ 25 ] . adult cardiac ventricular myocytes have been used as analysis tool in cardiovascular research for almost thirty years , and the popularity of this approach is constantly reinforced by the numerous studies published every year . however , during the last decade , the majority of researchers performing long - term ( longer than 1 week ) studies have favored the use of embryonic and neonatal cardiocytes . yet , changes in expression of ion channels and contractile protein isoforms during the development of the cardiocytes , pose a problem when making the extrapolation to the fully developed adult cardiocyte . the most popular methods include the ones that use laser diffraction techniques   and photodiode arrays , and those that employ the monitoring of microscopic cell images [ 911 ] . among the latter of these is the edge detection method that employs a raster - line to detect changes in myocyte length by sensing left and right cell boundaries using a threshold . this edge detection method is a widespread approach in research involving adult cardiocytes [ 4 , 10 , 12 , 13 ] . geometrical and boundary characteristics of adult cardiocytes are the most commonly irregular due to gap junction ends ( jagged edge ends ) , multiple intercalated disks , and variable cell widths , which can potentially complicate the application of the edge detection system on the cardiocyte . cardiocyte motion can also occur in unexpected directions since cardiocytes will commonly rotate or move vertically  perpendicular to the raster line  depending upon the location or absence of adhesion points . these conditions can lead to complications in the implementation of the edge detection system , and consequently , can result in an inaccurate analysis [ 4 , 12 , 13 ] . the cardiac myocyte is approximately 25  m in diameter and about 100  m in length . myofibrils have distinct , repeating microanatomical units , termed sarcomeres , which are the basic contractile units that make up a myocyte . the region of myofilament structures between two z - lines is defined as a sarcomere . the z - line 's associated structures are responsible for the lateral alignment of myofibrils and the generation of active contraction in cardiac muscles . the distance between z - lines  which is equivalent to the sarcomere length  ranges from about 1.6 to 2.2  m in human hearts . the sarcomere is composed of thick and thin filaments , myosin and actin , respectively . chemical and physical interactions between the actin and myosin cause the sarcomere length to shorten , allowing the myocyte to contract during the process of excitation - contraction coupling . contractility can be defined as the intrinsic ability of the heart muscle to generate force and to shorten . at the molecular level , the contractile process originates from the change in concentrations of calcium ( ca ) ions in the myocardial cytosol . ca ions enter through the calcium channel that opens in response to the wave of depolarization that travels along the sarcolemma . these ca ions trigger the release of additional calcium from the sarcoplasmic reticulum , thereby initiating a contraction - relaxation cycle . the need for an accurate method to assess different aspects of a myocyte has led researchers to explore several techniques in order to quantify contractility . some of these methods are not very popular due to the expensive equipment required , as in the case of the scanning ion conductance microscopy method . this technique , combined with laser confocal microscopy , measures myocyte height and local calcium concentration during contractility . other methods  such as light diffraction techniques  have been applied to the study of muscle mechanics since the nineteenth century . the reliability of these studies is relatively high , although they are highly dependent upon different factors . these include the temporal resolution of the detection system , sarcomere periodicity values , and other optical artifacts . the sarcomere striation pattern detection method has also been used as a way to quantify contractility . the method can achieve high temporal resolution with the aid of ccd line array detectors and it provides a measure of average sarcomere periodicity from the entire cell or cell regions . one drawback is this method 's vulnerability to errors introduced by cell geometry and rotational and translational changes which can occur during contraction . one of the first video - based efforts to measure contraction was performed with the assistance of a device capable of capturing the extent and rate of length shortening of isolated cardiac myocytes . the video - based method uses two tracking points at each end of the myocyte to track edge displacement as the myocyte contracts . the distance between the two edges is measured using edge detection while a record of the data is stored in a separate file . the method generally produces satisfactory results and has been an approved and widely used method for measuring contractile responses of adult myocytes for over twenty years [ 4 , 12 ] . several problems have been identified with the application of the video - based edge detection method for measuring adult myocyte contractility . the first inconvenience when analyzing myocyte contractility with this method is the need to have the cell positioned parallel to the raster - line . the myocyte should be perfectly positioned in the center of the screen ( parallel to the raster - line ) , and the proper threshold conditions must be set to detect the edges and follow them through a contraction . these threshold conditions are somewhat difficult to set depending upon the characteristics of the cell . the most important source of error that can be potentially introduced during the application of this method is the result of unexpected myocyte movements . myocytes will often rotate sideways or out of the plane of focus depending upon the presence or absence of adhesion points . the changes in myocyte geometry , dynamic torquing , and rotation can lead to errors in the experiment [ 2 , 4 , 12 ] . we propose a complete computational framework - based on well - established image processing techniques for the assessment of contractility of isolated cardiac myocytes . the proposed methodology is a multi - step process that provides a comprehensive account of the cardiac myocyte contraction process . the proposed computational framework for assessing the contractility in cardiac myocytes comprises the following stages : digital video recording of the contracting cell , edge preserving total variation - based image smoothing , segmentation of the smoothed images , contour extraction from the segmented images , shape representation by fourier descriptors , and contractility assessment . previous to the application of the proposed methodology , the specimen is appropriately prepared as follows . sprage - dawley rats purchased from harlan ( ca , usa ) were used for this study . after the isolation , the cells remained in the incuvator at least 12 hours in serum - free medium 199 ( gibco 12350 , invitrogen corporation , carlsbad , ca , usa ) before any measurements were performed . once ready for measurements , the cells were washed twice using serum - free media . cardiocytes were platted without laminin to ensure that the cells were free floating ( without adhesion points ) in the wells during the digital video recording . to capture the contraction process of the isolated cardiac myocytes the following procedure was employed . cells were placed in a chamber mounted on the stage of an inverted microscope ( nikon # elwd , nikon corporation , tokyo , japan ) . only rod - shaped myocytes with clear edges were selected for recording of mechanical properties . the cells were field stimulated with a suprathreshold ( 50% ) voltage at a frequency of 0.33  hz , for a 3 millisecond duration . the stimulation was performed using a pair of platinum wires placed on opposite sides of the chamber connected to an electrical stimulator ( grass sd9 , grass technologies , west warwick , ri , usa ) . the polarity of the stimulatory electrodes was reversed automatically every 10 stimuli to prevent electrode polarization . myocyte motion was digitally recorded with a camera ( pulnix tm-1327 , jai pulnix inc . , san jose , ca , usa ) mounted on the microscope , at a rate of 30  fps . we chose a total variation-(tv-)-based method for smoothing isotropic regions while preserving the cell 's edges in order to facilitate the segmentation step of the computational framework . rudin et al .   have argued that there are a number of reasons for preferring tv - based image smoothing models over their counterparts . tv - based algorithms are relatively simple to implement and result in minimal ringing ( nonoscillatory ) while recovering sharp edges ( noninvasive ) . in other words , the tv - norm allows piecewise smooth functions with jumps and is the proper space for the analysis and recovery of discontinuous functions . also , the tv - based formulations make a priori assumptions about the noise , and therefore they can be tailored to address the specific image restoration problem at hand . furthermore , empirical evidence suggests that  the human vision favors the l - norm  . in summary , the tv - based formulations seem to be a suitable approach for restoring piecewise continuous functions from noisy and blurry signals . appendix a provides a more detailed exposition of the original tv - based formulation due to rudin et al . along with some of the improvements proposed over the years , including the ones by two of the authors of this paper [ 2325 ] . the edge preserving tv - based image smoothing model used in our experiments is given by 
 ( 1)ut|u|(u|u|)+(uu0)=0 ,  on   [  0,),u(x,0)=u0(x ) ,  on  ,gu , n=0 ,  on  (0, ) , 
							 where u and u0 are the filtered and the observed images , respectively . the dynamic parameter  is defined as 
 ( 2)=122ut(uu0 ) ,  t 
							 with the approximation to the variance of the noise  given by 
 ( 3)var(t+1)=var(ut)var(gut ) ,  t , 
							 and the dynamic time step expressed as 
 ( 4)t=5+(145)(max  ( |u|)|u|max  ( |u| ) ) ,  t , 
							 with  = 1/255 . for more details on this model the reader is referred to . the algorithm to implement the edge preserving total variation - based image smoothing model is given in appendix b. figure 1 shows an example of the application of the edge preserving tv - based image smoothing model to a frame depicting an adult myocyte . segmentation of an image produces a set of labeled partitions or segments that represent the different components or features . this simplified image allows for an easier extraction of the main contours of the image . in our application this facilitates the identification of the contours of the cell that will permit the assessment of contractility of the cardiac myocyte . , the speed of execution becomes the principal constraint  because the segmentation of several hundreds of images is required . thus , we implemented a fast and robust segmentation technique - based on the one presented in . figure 2-top shows the segmented image of the cell after applying the segmentation procedure . for the contour extraction step of the framework we employ the built - in matlab ( the mathworks , inc . , natick , ma , usa ) function that creates a contour plot of image data . we convolve the extracted contour with a gaussian kernel to eliminate the typical noise produced in the segmentation procedure . we also resample the contour points in all the frames so that they will have the same number of contour points and for this to be equal to 2 , m  . this facilitates the implementation of the discrete fourier transform algorithm of the next section . fourier descriptors have been extensively proposed for the purpose of shape recognition , retrieval , classification , and analysis [ 2735 ] . among the contour - based shape representation methods the ones that have proven more promising for our application are nonetheless , the lower frequency fourier descriptors are the ones that contain the main information about the overall features of the shape . the higher frequency fourier descriptors , in turn , contain information relative to the finer details of the shape . therefore , only a relative small number of fourier descriptors are usually employed to capture the overall features of the shape . for completeness , we present both methods in the following subsections . let p(n ) = ( x(n ) , y(n ) ) , for 0  n  n  1 , be a discrete function representing the coordinates of a ( closed ) contour of an image 's shape in the cartesian space , such as the one in figure 2(b ) . in the complex plane , we can define this contour as a complex coordinates function q , such that 
 ( 5)q(n)=x(n)+jy(n ) . then , the fourier descriptors for the contour of the shape  described by q  can be computed using the discrete fourier transform ( dft ) . these fourier descriptors are the normalized fourier coefficients 
 ( 6)q(k)=1nn=0n1q(n)ej2kn / n ,  k=0,  ,n1 , 
							 which represent the discrete contour of the shape in the frequency domain [ 38 , 39 ] . in order to use the fourier descriptors as abstract representation of image features in each frame , it is customary to make them invariant to translation , scale , rotation , and their starting point . ( in some cases , retaining the step information can be advantageous . ) for our particular application , we want the fourier descriptors to change covariantly with the shape of the cell . in other words , we want the fourier coefficients to capture the contractions of the cell , disregarding only the translation and rotation of the cell . ( we also want to make the fourier descriptors independent of their starting point . ) translation of the contour function by    results in a change only to the first fourier descriptor , q(0 ) . therefore , by setting q(0 ) = 0 , we move the centroid of the contour onto 0 , and make the fourier descriptors invariant to translations . invariance with respect to the starting point can be achieved by subtracting the step of the second fourier descriptor , 1 = tan(im  q(1)/req(1 ) ) , weighted by k , from the step of all the fourier descriptors , q(k)e . rotation of the contour function by an angle  corresponds to a constant step shift of  in the fourier descriptors . we can make the fourier descriptors rotation invariant by computing the magnitude of the fourier descriptors , |q(k)| . the performance of this method is almost identical that of the centroid distance function method explained in the next subsection . a shape signature  a one - dimensional function derived from the shape boundary coordinates p(n ) = ( x(n ) , y(n ) ) , for 0  n  n  1can be used to generate fourier descriptors of the shape . fourier descriptors derived from centroid distance function generally outperform other shape signatures [ 36 , 40 ] . the centroid distance function of a shape is expressed by the distance of the points on the shape 's boundary from the centroid ( xc , yc ) of the shape 
 ( 7)r(n)=((x(n)xc)2+(y(n)yc)2)1/2 , 
							 where 
 ( 8)xc=1nn=0n1x(n ) ,  yc=1nn=0n1y(n ) . figure 3 shows the centroid distance function of the cell shape used in the description of these methods . since function ( 7 ) is real - valued , there are only n/2 distinct frequencies in the fourier transform . thus , only half of the fourier descriptors will be necessary to describe the shape . also , by construction therefore , we only need to make the fourier descriptors invariant to rotation and the starting point by identical procedures as in the case of the complex coordinates function method . figure 4 shows two identical cell shapes , one of which has been translated and rotated with respect to the other . along with the two shapes , figure 4 shows both of their first 30 fourier descriptors superimposed . we observe that both sets of fourier descriptors match almost perfectly for the case of translation , rotation , and starting point invariance . figure 5 shows two cell shapes of which one is slightly smaller and has been translated and rotated with respect to the other . along with the two shapes , figure 5 shows both of their first 30 fourier descriptors superimposed . we observe that their fourier descriptors are able to capture this change in shape size by making the fourier descriptors variant to scale but invariant to translation , rotation , and starting point . we tested the proposed approach by assessing the contractile responses in isolated adult rat cardiocytes , and compared them against the classic raster - line approach [ 911 ] . we used a sequence of digitized images obtained as previously described for both the proposed method and the raster - line technique . our results show good qualitative agreement between both methods as far as frequency , pacing , and overall behavior of the contractions are concerned ( see figure 6 ) . nonetheless , the raster - line method  being a one - dimensional technique  is unable to capture the contraction processes occurring outside its domain of influence . the proposed method , on the other hand , captures the contraction of the cell as a two dimensional event over the entire boundary of the cell . the proposed methodology was also able to capture a slower recovery period than the raster - line method ( see figure 7 ) , which can be attributed to the dimensionality characteristics of both methods . this means that the proposed method is capable of not only assessing the myocyte 's length , but also its overall changes in shape and geometry . in other words , it is capable of assessing the myocyte 's dimensional changes during contraction while remaining invariant to rotation , translation and starting point . we presented a complete computational framework for the comprehensive assessment of contractile responses of isolated adult cardiac myocytes . the proposed methodology comprises the following stages : digital video recording of the contracting cell , edge preserving total variation - based image smoothing , segmentation of the smoothed images , contour extraction from the segmented images , shape representation by fourier descriptors , and contractility measurements . these stages are - based on mathematically sound and computational robust algorithms that are very well established in the image processing community . being a two - dimensional technique  is capable of capturing the contractile processes that are otherwise missed by the one - dimensional techniques . this capability makes the method suitable for determining myocyte contraction in cells that usually bend or move during contraction , for example , atrial myocytes and isolated smooth muscle cells , or in cardiac myocytes which develop spatially nonuniform oscillatory contractile activity induced by intracellular calcium fluctuations [ 10 , 41 ] . our future work entails the application of the proposed method to analyzing the contractility of myocytes that have been exposed to a drug over a given period of time . we have been investigating one category of the mechanisms that may be responsible for the observed effects on heart cells from a synthetic antidiabetic drug , rosiglitazone ( avandia , glaxosmithkline , brentford , uk ) of the thiazolidinedione ( tzd ) family of insulin - sensitizing compounds used in the treatment of type ii diabetes . we are anticipating that the proposed method will be an essential tool in that it will complement the analysis of our drug studies , which have been also performed using microarray , ca transient , gene and protein expression measurements . furthermore , we are in the process of deploying a more sophisticated image acquisition technology that includes a high - speed camera . this will allow for a more in - depth analysis of the contraction process undergone by the cardiac myocyte . rudin et al .   proposed removing noise from images by minimizing the tv norm of the estimated solution . they derived a constrained minimization algorithm as a time - dependent nonlinear pde , where the constraints are determined by the noise statistics . they stated that the space of bounded ( total ) variation ( bv ) is the proper class for many basic image processing tasks . thus , given a noisy image u0 = f +  , where the true image f has been perturbed by additive white noise  , the restored image u  f is the solution of 
 ( a.1)min  ubv()tv(u)=min  ubv()|u|dx , 
							 subject to the following constraints involving the noise : 
								 ( a.2)12  (uu0)2dx=12||2,1||u0  dx1||u  dx , 
							 where || represents the area of the image . the first constraint uses a priori information that the standard deviation of the noise is  , while the second constraint assumes that the noise has zero mean . ( it can also be shown that ( a.1 ) and ( a.2 ) imply that the noise is normally distributed . ) the tv - norm does not penalize discontinuities in u , and ; therefore , it allows the recovery of the edges in the observed image u0 . to solve this minimization problem we would usually solve its associated euler - lagrange equation , namely , 
								 ( a.3)(u|u|)+(uu0)=0 , 
							 on a closed domain  , and subject to homogeneous neumann boundary conditions on the boundary . the solution procedure proposed in   employs a parabolic equation with time as an evolution ( scale ) parameter , or equivalently , the gradient descent method , that is , 
								 ( a.4)ut(u|u|)+(uu0)=0 , 
							 for t > 0 , on a closed domain  , with the observed image as initial condition , u(x , 0 ) = u0(x ) , and homogeneous neumann boundary conditions , g  u , n = 0 , on the boundary . for the parameter  , they suggested a dynamic value (t ) estimated by rosen 's gradient - projection method , which as t converges to 
								 ( a.5)=12||2    [ |u|u0tu|u|]dx . 
							 existence and other successful implementations of this minimization problem include the second order cone programming , convex programming , duality , and a fast and exact minimization method - based on graph cuts [ 47 , 48 ] . nonetheless , the rudin - osher - fatemi model , in its original form , presents several practical challenges . this evolution scheme is not trivial to implement since it is highly nonlinear and not well - posed . it can also run into trouble when |u|  0 beyond machine accuracy . in practice , it is very common to use a slightly modified version of the tv - norm : 
								 ( a.6)  (|u|2+)1/2  dx , 
							 where  is a small positive number which  smoothes out the corner  at |u| = 0 . the two other practical ( observable ) limitations presented by the rudin - osher - fatemi original model are the loss of contrast [ 51 , 52 ] and the  staircase  effect , that is , a strong preference for piecewise constant patches [ 53 , 54 ] . the rudin - osher - fatemi model has been extensively studied and improved upon by many scientists [ 20 , 4952 , 5465 ] . marquina and osher   proposed a different version of the transient parabolic equation that helps speed up the convergence of the time - marching scheme . the new evolution equation is 
								 ( a.7)ut|u|(u|u|)+|u|g(guu0)=0 , 
							 for t > 0 , on a closed domain  , with the observed image as initial condition , u(x , 0 ) = u0(x ) , and homogeneous neumann boundary conditions , g  u , n = 0 , on the boundary  , and where g is a blurring operator ( gaussian kernel ) . this approach fixes the staircase problem of the original scheme and is used for the removal of both blur and noise . version of the rudin   euler - lagrange equation as modified by marquina and osher . they referred to this approach as parameter - free adaptive total variation - based noise removal and edge strengthening model . for our current application , we will use this method for estimating the unknown noise level and their pixel - wise definition for the parameter . given the assumption that the image has been perturbed by additive white noise , u0 = f +  , and that this noise is independent from the signal , the variance of the noisy image must be equal to the sum of the variance of the true image and the variance of the noise , that is , var(u0 ) = var(gu0 ) + var( ) . here , the variance of the ( unknown ) true image is approximated by the variance of the convolved noisy image with a gaussian kernel of width  = 1 . this parameter will be updated at every iteration which provides a positive effect . for the parameter they proposed a variation of the method suggested in . instead of integrating ( or summing ) over the domain  , they assumed a pixel - wise parameter as 
								 ( a.8)|u|=122[ux(ux(u0)x)+uy(uy(u0)y ) ] . for a detailed explanation of the attributes of the dynamic parameter  , the reader can consult . due to the high nonlinearity of the tv - based models , to ensure stability song   has shown that the cfl condition for the rudin - osher - fatemi model is t/x he has also shown that the cfl condition for the marquina - osher model is t/x  c , with c > 0 . as a rule of thumb , gilboa has suggested ( assuming  x = 1 ) setting the value of t = /5 , where  is the regularization constant used in ( a.6 ) . have shown that for explicit discretization schemes , the stability condition for the perona - malik - type models ( assuming  x = 1 and for  all  s : g(s )  1 ) is t < 1/2d , with d being the number of dimensions of the data . in his dissertation , bazn   argued that since the explicit discretization schemes used in the tv - based models produce updates of the following form : 
								 ( a.9)ut+1=ut+tf(ut,ut , u0,u0, ) , 
							 then , in practice , the smaller the time - step , the slower the restoration process . t(x , t ) , which does not only make the tv - based schemes more stable ( smooth ) , but also speeds up the restoration process . the proposed adaptive time - step is 
								 ( a.10)t(x)=5+(12d5)(max  ( |u|)|u|max  ( |u| ) ) ,  t , 
							 where , as before ,  is the regularization constant used in ( a.6 ) , and d is the number of dimensions of the data . for a detailed description of the characteristics of the adaptive time - step the reader is referred to . the   algorithm to implement the parameter - free adaptive tv - based noise removal and edge strengthening model is in table 1 .
we describe a computational framework for the comprehensive assessment 
 of contractile responses of enzymatically dissociated adult cardiac myocytes . 
 the proposed methodology comprises the following stages : digital video recording of the contracting cell , edge preserving total variation - based image 
 smoothing , segmentation of the smoothed images , contour extraction from the segmented images , shape representation by fourier descriptors , and contractility assessment . 
 the different stages are variants of mathematically 
 sound and computationally robust algorithms very well established in the image processing community . 
 the physiologic application of the methodology is evaluated by assessing overall contraction in enzymatically dissociated adult rat cardiocytes . 
 our results demonstrate the effectiveness of the proposed approach in characterizing the true , two - dimensional ,  shortening  in the contraction process of adult cardiocytes . 
 we compare the performance of the proposed method to that of a popular edge detection system in the literature . 
 the proposed method not only provides a more comprehensive assessment of the myocyte contraction process but also can potentially eliminate historical concerns and sources of errors caused by myocyte rotation or translation during contraction . 
 furthermore , the versatility of the image processing techniques makes the method suitable for determining myocyte shortening in cells that usually bend or move during contraction . 
 the proposed method can be utilized to evaluate changes in contractile behavior resulting from drug intervention , disease modeling , transgeneity , or other common applications to mammalian cardiocytes .