Method GetPointGradient

Common/DataModel/vtkImageData.cxx:570–669

------------------------------------------------------------------------------ Given structured coordinates (i,j,k) for a point in a structured point dataset, compute the gradient vector from the scalar data at that point. The scalars s are the scalars from which the gradient is to be computed. This method will treat structured point datasets of any dimension.

// The scalars s are the scalars from which the gradient is to be computed.
// This method will treat structured point datasets of any dimension.
void vtkImageData::GetPointGradient(int i, int j, int k, vtkDataArray* s, double g[3])
{
  const double* ar = this->Spacing;
  double sp, sm;
  const int* extent = this->GetExtent();

  vtkIdType dims[3];
  this->GetDimensions(dims);
  vtkIdType ijsize = dims[0] * dims[1];

  // Adjust i,j,k to the start of the extent
  i -= extent[0];
  j -= extent[2];
  k -= extent[4];

  // Check for out-of-bounds
  if (i < 0 || i >= dims[0] || j < 0 || j >= dims[1] || k < 0 || k >= dims[2])
  {
    g[0] = g[1] = g[2] = 0.0;
    return;
  }

  // i-axis
  if (dims[0] == 1)
  {
    g[0] = 0.0;
  }
  else if (i == 0)
  {
    sp = s->GetComponent(i + 1 + j * dims[0] + k * ijsize, 0);
    sm = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
    g[0] = (sm - sp) / ar[0];
  }
  else if (i == (dims[0] - 1))
  {
    sp = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
    sm = s->GetComponent(i - 1 + j * dims[0] + k * ijsize, 0);
    g[0] = (sm - sp) / ar[0];
  }
  else
  {
    sp = s->GetComponent(i + 1 + j * dims[0] + k * ijsize, 0);
    sm = s->GetComponent(i - 1 + j * dims[0] + k * ijsize, 0);
    g[0] = 0.5 * (sm - sp) / ar[0];
  }

  // j-axis
  if (dims[1] == 1)
  {
    g[1] = 0.0;
  }
  else if (j == 0)
  {
    sp = s->GetComponent(i + (j + 1) * dims[0] + k * ijsize, 0);
    sm = s->GetComponent(i + j * dims[0] + k * ijsize, 0);
    g[1] = (sm - sp) / ar[1];
  }
  else if (j == (dims[1] - 1))