Document: distance calculation for morphological erosion

[graeme-winter]

dials/src/dials/algorithms/image/filter/distance.h

Line 67 in 9af1ec1

void chebyshev_distance(const af::const_ref<InputType, af::c_grid<2> > &src,

This is supposed to compute the chess board king move distance between a pixel A and the nearest pixel B of a given class, and checks North / East then South / West. As far as I can tell this ignores pixels in the North East or the South / West quadrant?

This feels like a bug to me. But I don't think it is a significant bug, as the erosion probably achieved by other pixel combinations. I only stumbled across this as appearing asymmetric when I was trying to understand how the dispersion extended spot finding actually worked.

[graeme-winter]

OK, comments also mention other quadrants, so more a case of "I don't understand this code" than "this is a bug"

[graeme-winter]

It turns out it is not a bug, but the code could do with some annotation

#include <iostream>
#include <vector>

#define NX 64
#define NY 64

std::vector<uint16_t> dumb(std::vector<uint16_t> image) {
  std::vector<uint16_t> result(NY * NX);
  for (int i = 0; i < NY; i++) {
    for (int j = 0; j < NX; j++) {
      uint16_t mask = image[i * NX + j];
      for (int _i = -3; _i <= 3; _i++) {
        if ((i + _i) < 0 || (i + _i) >= NY)
          continue;
        for (int _j = -3; _j <= 3; _j++) {
          if ((j + _j) < 0 || (j + _j) >= NX)
            continue;
          mask = std::min(mask, image[(i + _i) * NX + j + _j]);
        }
      }
      result[i * NX + j] = mask;
    }
  }
  return result;
}

std::vector<uint16_t> erode(std::vector<uint16_t> dst, uint16_t distance) {
  std::vector<uint16_t> mask(NY * NX);
  for (int i = 0, k = 0; i < NY; i++) {
    for (int j = 0; j < NX; j++, k++) {
      if (dst[k] <= distance) {
        mask[k] = 0;
      } else {
        mask[k] = 1;
      }
    }
  }
  return mask;
}

std::vector<uint16_t> distance(std::vector<uint16_t> image) {
  std::vector<uint16_t> dst(NY * NX);

  uint16_t max_distance = NY + NX;

  for (int i = 0, k = 0; i < NY; i++) {
    for (int j = 0; j < NX; j++, k++) {
      uint16_t N = (i > 0) ? dst[(i - 1) * NX + j] : max_distance;
      uint16_t W = (j > 0) ? dst[i * NX + j - 1] : max_distance;
      uint16_t NW = (i > 0 && j > 0) ? dst[(i - 1) * NX + j - 1] : max_distance;
      uint16_t NE =
          (i > 0 && j < NX - 1) ? dst[(i - 1) * NX + j + 1] : max_distance;
      if (image[k] == 0) {
        dst[k] = 0;
      } else {
        dst[k] = 1 + std::min(std::min(N, W), std::min(NE, NW));
      }
    }
  }

  for (int i = NY - 1, k = NY * NX - 1; i >= 0; i--) {
    for (int j = NX - 1; j >= 0; j--, k--) {
      uint16_t S = (i < NY - 1) ? dst[(i + 1) * NX + j] : max_distance;
      uint16_t E = (j < NX - 1) ? dst[i * NX + j + 1] : max_distance;
      uint16_t SW =
          (i < NY - 1 && j < NX - 1) ? dst[(i + 1) * NX + j + 1] : max_distance;
      uint16_t SE =
          (i < NY - 1 && j > 0) ? dst[(i + 1) * NX + j - 1] : max_distance;
      uint16_t X = 1 + std::min(std::min(S, E), std::min(SW, SE));
      if (image[k] == 0) {
        dst[k] = 0;
      } else {
        dst[k] = std::min(X, dst[k]);
      }
    }
  }

  return dst;
}

int main(int argc, char **argv) {
  std::vector<uint16_t> image(NY * NX);
  for (int i = 0; i < NY; i++) {
    for (int j = 0; j < NX; j++) {
      if ((j - 32) * (j - 32) + (i - 32) * (i - 32) < 256) {
        image[i * NX + j] = 1;
      } else {
        image[i * NX + j] = 0;
      }
    }
  }

  std::vector<uint16_t> dst = distance(image);
  std::vector<uint16_t> sqr = dumb(image);
  std::vector<uint16_t> erd = erode(dst, 3);

  for (int i = 0, k = 0; i < NY; i++) {
    for (int j = 0; j < NX; j++, k++) {
      if (erd[k] != sqr[k])
        std::cout << k << " " << erd[k] << " " << sqr[k] << std::endl;
    }
  }

  return 0;
}

Reimplementation demonstrates it does what it says on the tin.

Documentation needed to explain this => reopen and relabel.

[graeme-winter]

https://github.com/dials/dials/wiki/How-Does-Spot-Finding-Work