Function otsuThreshold

src/api/c/canny.cpp:82–144

}

Array<float> otsuThreshold(const Array<float>& in, const unsigned NUM_BINS,
                           const float maxVal) {
    Array<uint> hist = histogram<float>(in, NUM_BINS, 0, maxVal, false);

    const dim4& inDims = in.dims();
    const dim4& hDims  = hist.dims();

    const dim4 oDims(1, hDims[1], hDims[2], hDims[3]);
    vector<af_seq> seqBegin(4, af_span);
    vector<af_seq> seqRest(4, af_span);
    vector<af_seq> sliceIndex(4, af_span);

    seqBegin[0] = af_make_seq(0, static_cast<double>(hDims[0] - 1), 1);
    seqRest[0]  = af_make_seq(0, static_cast<double>(hDims[0] - 1), 1);

    Array<float> UnitP  = createValueArray<float>(oDims, 1.0f);
    Array<float> histf  = cast<float, uint>(hist);
    Array<float> totals = createValueArray<float>(hDims, inDims[0] * inDims[1]);
    Array<float> weights =
        iota<float>(dim4(NUM_BINS), oDims);  // a.k.a histogram shape

    // pixel frequency probabilities
    auto freqs        = arithOp<float, af_div_t>(histf, totals, hDims);
    auto cumFreqs     = scan<af_add_t, float, float>(freqs, 0);
    auto oneMCumFreqs = arithOp<float, af_sub_t>(UnitP, cumFreqs, hDims);
    auto qLqH         = arithOp<float, af_mul_t>(cumFreqs, oneMCumFreqs, hDims);
    auto product      = arithOp<float, af_mul_t>(weights, freqs, hDims);
    auto cumProduct   = scan<af_add_t, float, float>(product, 0);
    auto weightedSum  = reduce<af_add_t, float, float>(product, 0);

    dim4 sigmaDims(NUM_BINS - 1, hDims[1], hDims[2], hDims[3]);
    Array<float> sigmas = createEmptyArray<float>(sigmaDims);
    for (unsigned b = 0; b < (NUM_BINS - 1); ++b) {
        const dim4 fDims(b + 1, hDims[1], hDims[2], hDims[3]);
        const dim4 eDims(NUM_BINS - 1 - b, hDims[1], hDims[2], hDims[3]);

        sliceIndex[0]    = {double(b), double(b), 1};
        seqBegin[0].end  = static_cast<double>(b);
        seqRest[0].begin = static_cast<double>(b + 1);

        auto qL    = createSubArray(cumFreqs, sliceIndex, false);
        auto qH    = arithOp<float, af_sub_t>(UnitP, qL, oDims);
        auto _muL  = createSubArray(cumProduct, sliceIndex, false);
        auto _muH  = arithOp<float, af_sub_t>(weightedSum, _muL, oDims);
        auto muL   = arithOp<float, af_div_t>(_muL, qL, oDims);
        auto muH   = arithOp<float, af_div_t>(_muH, qH, oDims);
        auto diff  = arithOp<float, af_sub_t>(muL, muH, oDims);
        auto sqrd  = arithOp<float, af_mul_t>(diff, diff, oDims);
        auto op2   = createSubArray(qLqH, sliceIndex, false);
        auto sigma = arithOp<float, af_mul_t>(sqrd, op2, oDims);

        auto binRes = createSubArray<float>(sigmas, sliceIndex, false);
        copyArray(binRes, sigma);
    }

    Array<float> thresh = createEmptyArray<float>(oDims);
    Array<uint> locs    = createEmptyArray<uint>(oDims);