Areas and centroids

moments.cpp

array areas(const array& in)
{
    // Must be 2d image
    dim4 dims = in.dims();

    array flt_in = flat(in);

    // Corner case: no components
    array non_zero_ind = where(flt_in != 0);
    if (non_zero_ind.isempty())
        return non_zero_ind;

    // Otherwise sort the array
    array vals;
    array idx;
    sort(vals, idx, flt_in);
    array locs = where(diff1(vals) != 0);

    // Corner case: 1 giant component
    if (locs.isempty()) {
        return join(1,
                    vals(0),
                    constant(static_cast<float>(flt_in.dims(0)), 1, 1));
    }

    locs = join(0,
                locs,
                constant(static_cast<float>(flt_in.dims(0)-1), 1));

    array areas  = diff1(locs);
    array labels = flt_in((idx(locs)));

    // Remove the entry for zero
    labels = labels(seq(1, end));
    return join(1, labels, areas);
}

array centroids(const array& in)
{
    dim4 dims = in.dims();

    array flt_in = flat(in);

    // Corner case: no components, return empty array
    array non_zero_ind = where(flt_in != 0);
    if (non_zero_ind.isempty())
        return non_zero_ind;

    array vals;
    array idx;
    sort(vals, idx, flt_in);

    // Compute x/y indices
    const unsigned m = in.dims(0);
    array y = floor(idx/ m);
    array x = idx - (y * m);

    // Scan indices for the same values
    array scanned_x = segsum(x, vals);
    array scanned_y = segsum(y, vals);
    array locs = where(diff1(vals) != 0);

    // Corner case: 1 giant component, return center of img
    if (locs.isempty()) {
        const float x_half = static_cast<float>(in.dims(0)) / 2.0f;
        const float y_half = static_cast<float>(in.dims(1)) / 2.0f;
        array xa = constant(x_half, 1, 1);
        array ya = constant(y_half, 1, 1);
        return join(1, join(1, vals(0), xa), ya);
    }

    locs = join(0,
                locs,
                constant(static_cast<float>(flt_in.dims(0)-1), 1));

    // Get areas and labels, trim zeros if not using diff1
    array areas  = diff1(locs);
    locs = locs(seq(1, end));
    array labels = flt_in((idx(locs)));

    // Build answer array
    array ans = join(1, labels, scanned_x(locs) / areas);
    ans       = join(1, ans,    scanned_y(locs) / areas);
    return ans;
}