[Insight-developers] CovarianceCalculator::ComputeCovarianceWithoutGivenMean Bug with fix

Bradley Lowekamp blowekamp at mail.nih.gov
Thu Jun 12 16:08:15 EDT 2008


Greeting,

I was using this class on scalar histograms of a fixed range [-0.5,  
255.5], which happen to have 0 frequency in the first few bins. I was  
getting nan results. Upon looking at this function it was obvious  
where the problem was with a divide by zero.

Something to the effect of:

if (frequency == 0)
  continue;

needs to be added early on in this loop in the  
CovarianceCalculator::ComputeCovarianceWithoutGivenMean method.


  //
  // fills the lower triangle and the diagonal cells in the covariance  
matrix
   while (iter != end)
     {
     frequency = iter.GetFrequency() ;
     totalFrequency += frequency ;
     measurements = iter.GetMeasurementVector() ;
     for ( i = 0 ; i < measurementVectorSize ; ++i )
       {
       diff[i] = measurements[i] - (*m_InternalMean)[i] ;
       }

     // updates the mean vector
     double tempWeight = frequency / totalFrequency ;
     for ( i = 0 ; i < measurementVectorSize ; ++i )
       {
       (*m_InternalMean)[i] += tempWeight * diff[i] ;
       }

     // updates the covariance matrix
     tempWeight = tempWeight * ( totalFrequency - frequency ) ;
     for ( row = 0; row < measurementVectorSize ; row++ )
       {
       for ( col = 0; col < row + 1 ; col++)
         {
         m_Output(row,col) +=
           tempWeight * diff[row] * diff[col] ;
         }
       }
     ++iter ;
     }

The itkWeightedCovarianceCalculator class looks like it has this bug  
as well but in two places in similar loops.

Enjoy!

========================================================
Bradley Lowekamp
Lockheed Martin Contractor for
Office of High Performance Computing and Communications
National Library of Medicine
blowekamp at mail.nih.gov


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20080612/d1377784/attachment.htm>


More information about the Insight-developers mailing list