py::tuple hahog()

in opensfm/src/features/src/hahog.cc [123:202]

• 62 lines of code
• 3 McCabe index (conditional complexity)

py::tuple hahog(foundation::pyarray_f image, float peak_threshold,
                float edge_threshold, int target_num_features) {
  if (!image.size()) {
    return py::none();
  }

  std::vector<float> points;
  std::vector<float> desc;
  vl_size numFeatures;
  vl_size dimension = 128;

  {
    py::gil_scoped_release release;

    // create a detector object
    VlCovDet *covdet = vl_covdet_new(VL_COVDET_METHOD_HESSIAN);
    // set various parameters (optional)
    vl_covdet_set_first_octave(covdet, 0);
    vl_covdet_set_peak_threshold(covdet, peak_threshold);
    vl_covdet_set_edge_threshold(covdet, edge_threshold);

    // process the image and run the detector
    vl_covdet_put_image(covdet, image.data(), image.shape(1), image.shape(0));
    vl_covdet_set_non_extrema_suppression_threshold(covdet, 0);
    vl_covdet_detect(covdet, std::numeric_limits<vl_size>::max());

    // select the best features to keep
    numFeatures = run_features_selection(covdet, target_num_features);

    // compute the orientation of the features (optional)
    std::vector<VlCovDetFeature> vecFeatures =
        vlfeat_covdet_extract_orientations(covdet, numFeatures);
    numFeatures = vecFeatures.size();

    // get feature descriptors
    VlSiftFilt *sift = vl_sift_new(16, 16, 1, 3, 0);
    vl_index i;
    vl_index patchResolution = 15;
    double patchRelativeExtent = 7.5;
    double patchRelativeSmoothing = 1;
    vl_size patchSide = 2 * patchResolution + 1;
    double patchStep = (double)patchRelativeExtent / patchResolution;
    points.resize(4 * numFeatures);
    desc.resize(dimension * numFeatures);
    std::vector<float> patch(patchSide * patchSide);
    std::vector<float> patchXY(2 * patchSide * patchSide);

    vl_sift_set_magnif(sift, 3.0);
    for (i = 0; i < (signed)numFeatures; ++i) {
      const VlFrameOrientedEllipse &frame = vecFeatures.at(i).frame;
      float det = frame.a11 * frame.a22 - frame.a12 * frame.a21;
      float size = sqrt(fabs(det));
      float angle = atan2(frame.a21, frame.a11) * 180.0f / M_PI;
      points[4 * i + 0] = frame.x;
      points[4 * i + 1] = frame.y;
      points[4 * i + 2] = size;
      points[4 * i + 3] = angle;

      vl_covdet_extract_patch_for_frame(covdet, &patch[0], patchResolution,
                                        patchRelativeExtent,
                                        patchRelativeSmoothing, frame);

      vl_imgradient_polar_f(&patchXY[0], &patchXY[1], 2, 2 * patchSide,
                            &patch[0], patchSide, patchSide, patchSide);

      vl_sift_calc_raw_descriptor(
          sift, &patchXY[0], &desc[dimension * i], (int)patchSide,
          (int)patchSide, (double)(patchSide - 1) / 2,
          (double)(patchSide - 1) / 2,
          (double)patchRelativeExtent / (3.0 * (4 + 1) / 2) / patchStep,
          VL_PI / 2);
    }
    vl_sift_delete(sift);
    vl_covdet_delete(covdet);
  }

  return py::make_tuple(
      foundation::py_array_from_data(points.data(), numFeatures, 4),
      foundation::py_array_from_data(desc.data(), numFeatures, dimension));
}