[Insight-users] Expensive Gaussian calls

[Peter Cech]

On Mon, Sep 14, 2009 at 09:26:18PM +0200, motes motes wrote:
> I iterate a 512*512 image.
> 
> For each pixel I compute the gaussian of a distance measure (computed
> in constant time). But since the gaussian is rather expensive to
> evaluate I first do:
> 
>     //GaussianPoint gp;
>     //gp[0] = dist;
>     if (m_Gaussian->IsInside(gp)) {
> 
>         double weight = 0.0;
>         m_Gaussian->ValueAt(gp, weight);
>     }
> 
> 
> I have tried just to do the:
> 
> m_Gaussian->IsInside(gp)
> 
> call and it adds around 10 seconds to the execution time (I am using a
> release build). And the m_Gaussian->ValueAt(gp, weight); adds even
> more!
> 
> What makes:
> 
> m_Gaussian->IsInside(gp)
> m_Gaussian->ValueAt(gp, weight);
> 
> so expensive??

The above description lacks some information. I'll give it a stab, but
in general, you'll have better chances to get answers if you provide
complete information:

 - What are the types of variables and classes here? (m_Gaussian,
   GaussianPoint)
 - A frame of reference for the additional runtime. What is the runtime
   without the call? How many times is the call made?
 - Some more details what you want to do. There might be a simpler/more
   efficient way to do it.

Now to my guess. I'll assume that m_Gaussian is of type
GaussianSpatialObject. The ValueAt method calls IsInside, so you can
change your code to:

  double weight;
  if (!m_Gaussian->ValueAt(gp, weight)) {
    // gp is not inside, weight was set to a special value marking
    // outside, but we want 0.0
    weight = 0.0;
  }

The methods are very expensive, because the spatial object can be
transformed in space and the point has to be transformed with inverse
transform before actual computations on it are done. If you do not need
the flexibility of spatial transforms (m_Gaussian is always isotrpoic),
use Statistics::GaussianDistribution,
Statistics::GaussianTransferFunction or compute the gaussian directly.

Peter