void detectDots()

in IsometricPatternMatcher/DotExtractor.h [138:221]
66 lines of code
18 McCabe index (conditional complexity)

  void detectDots(const Eigen::AlignedBox2i roi) {
    std::atomic<INDEX_TYPE> atomicIndex(0);
    std::atomic<INDEX_TYPE> rejectedAtomicIndex(0);
    for (int y = roi.min().y(); y < roi.max().y(); ++y) {
      for (int x = roi.min().x(); x < roi.max().x(); ++x) {
        const auto* rm = blurredImage_.RowPtr(y - 1) + x;
        const auto* r0 = blurredImage_.RowPtr(y) + x;
        const auto* rp = blurredImage_.RowPtr(y + 1) + x;

        const float ixy = *r0;

        /* Ignore non-maximal pixels.

          Calculate and(ixy >= X) where X is all 8 other pixels surrounding the
         center pixel.
        */

        bool isMaximalPixel = ixy > rm[-1] && ixy > rm[0] && ixy > rm[1] &&
                              ixy > r0[-1] && ixy > r0[1] && ixy > rp[-1] &&
                              ixy > rp[0] && ixy > rp[1];

        if (isMaximalPixel) {
          const float dyy = (*rp + *rm - 2.0f * ixy);
          const float dxx = (r0[1] + r0[-1] - 2.0f * ixy);
          const float dxy = (rm[-1] - rm[1] - rp[-1] + rp[1]) / 4.0f;

          const float det = (dxx * dyy - dxy * dxy);

          if (det != 0.0f) {
            const float dx = (r0[1] - r0[-1]) / 2.0f;
            const float dy = (*rp - *rm) / 2.0f;
            const float invdet = 1.0f / det;

            // x,y update (newton step)
            Eigen::Vector2f fp((dxy * dy - dyy * dx) * invdet,
                               (dxy * dx - dxx * dy) * invdet);

            if (fp.lpNorm<Eigen::Infinity>() < maxDelta_) {
              if (dxx * dxx + dyy * dyy > hessThresh_) {
                INDEX_TYPE memoryIndex = atomicIndex++;
                if (memoryIndex < maxNumDots_) {
                  if (clampDelta_ > 0.f) {
                    fp[0] = std::clamp(fp[0], -clampDelta_, clampDelta_);
                    fp[1] = std::clamp(fp[1], -clampDelta_, clampDelta_);
                  }

                  dots_[memoryIndex] =
                      Eigen::Matrix<float, 2, 1>(x + fp[0], y + fp[1]);

                } else {  // maxdots
                  INDEX_TYPE idx = rejectedAtomicIndex++;
                  rejectedDots_[idx] = Eigen::Vector2f(x + fp[0], y + fp[1]);
                  rejectedStatus_[idx] = RejectedDotStatus::MAXNUM;
                }
              } else {
                INDEX_TYPE idx = rejectedAtomicIndex++;
                rejectedDots_[idx] = Eigen::Vector2f(x + fp[0], y + fp[1]);
                rejectedStatus_[idx] = RejectedDotStatus::HESSIAN;
              }
            } else {
              INDEX_TYPE idx = rejectedAtomicIndex++;
              rejectedDots_[idx] = Eigen::Vector2f(x + fp[0], y + fp[1]);
              rejectedStatus_[idx] = RejectedDotStatus::MAXDELTA;
            }
          } else {
            INDEX_TYPE idx = rejectedAtomicIndex++;
            rejectedDots_[idx] = Eigen::Vector2f(x, y);
            rejectedStatus_[idx] = RejectedDotStatus::DETERMINANT;
          }
        }
      }
    }

    INDEX_TYPE ai = atomicIndex == 0
                        ? 0
                        : atomicIndex - 1;  // to compensate for the last ++
    numDots_ = std::clamp(ai, static_cast<INDEX_TYPE>(0), maxNumDots_);

    INDEX_TYPE rai =
        rejectedAtomicIndex == 0
            ? 0
            : rejectedAtomicIndex - 1;  // to compensate for the last ++
    numRejectedDots_ = std::clamp(rai, static_cast<INDEX_TYPE>(0), maxNumDots_);
  }