Median difference

[ns-rse]

As noted in #52 which revises the filtering process which differs somewhat from processing DNA images we need to first
apply a Median difference which...

shifts the lines so that the median of height differences (between vertical neighbour pixels) becomes zero. Therefore
it better preserves large features while it is more sensitive to completely bogus lines.

Options

Scikit-image has
skimage.filters.median()
function which "Return[s] the local median of an image." and it has two "behaviour"/methods, rank which uses the
native skimge.filters.rank.median() and ndimage which calls SciPy
scipy.ndimage.median_filter()
which "Calculate[s] a multidimensional median filter.". This functions has the option axes which...

If None, input is filtered along all axes. Otherwise, input is filtered along the specified axes. When axes is
specified, any tuples used for size, origin, and/or mode must match the length of axes. The ith entry in any of these
tuples corresponds to the ith entry in axes.

The function also notes...

For 2-dimensional images with uint8, float32 or float64 dtypes the specialised function
scipy.signal.medfilt2d
may be faster. It is however limited to constant mode with cval=0.

cval is a value to fill past edges if the mode is constant which adds a fixed value beyond the boundaries, it
defaults to 0.0.

Checking scipy.signal.medfilt2d it
takes a 2-d array as input and then requires the kernel-size which is "the size of the median filter window in eeach
dimension. Elements of kernel_size should be odd". Whether this permits 0 as a parameter so that the kernel only
extends in the vertical plane will require investigating.

Regardless I think we might be able to use these functions to perform the required median difference filter to the
image but I will have to construct some simple examples (i.e. tests!) to check that they are doing what we require. If
not we'll have to develop the function from scratch.

[ns-rse]

Started working on this and went to check how to setout factory methods in TopoStats and discovered we already have a
Filters.median_flatten()
method in TopoStats! 🙌

We should be able to use this, although I will implement tests as the tests in TopoStats are slightly occluded using
random numbers in a 20x20 array.

[ns-rse]

Had a useful conversation with @SylviaWhittle (thanks Sylvia!) and this isn't quite the same as what Gwyddion does so
I'm going to check the source code of Gwyddion to see exactly how the Median difference is implemented and compare it
to that in TopoStats and the methods available in SciPy.

Gwyddion

Not very conversant in C code but I think I've found where this is done.

modules/process/linematch.c

This has a set of enum

 47 │ typedef enum {
 48 │     LINE_MATCH_POLY         = 0,
 49 │     LINE_MATCH_MEDIAN       = 1,
 50 │     LINE_MATCH_MEDIAN_DIFF  = 2,
 51 │     LINE_MATCH_MODUS        = 3,
 52 │     LINE_MATCH_MATCH        = 4,
 53 │     LINE_MATCH_TRIMMED_MEAN = 5,
 54 │     LINE_MATCH_TMEAN_DIFF   = 6,
 55 │     LINE_MATCH_FACET_TILT   = 7,
 56 │ } LineMatchMethod;

...which map to the options when clicking on the "Align Rows"...

121 │ static const GwyEnum methods[] = {
122 │     { N_("linematch|Polynomial"),        LINE_MATCH_POLY,         },
123 │     { N_("Median"),                      LINE_MATCH_MEDIAN,       },
124 │     { N_("Median of differences"),       LINE_MATCH_MEDIAN_DIFF,  },
125 │     { N_("Modus"),                       LINE_MATCH_MODUS,        },
126 │     { N_("linematch|Matching"),          LINE_MATCH_MATCH,        },
127 │     { N_("Trimmed mean"),                LINE_MATCH_TRIMMED_MEAN, },
128 │     { N_("Trimmed mean of differences"), LINE_MATCH_TMEAN_DIFF,   },
129 │     { N_("Facet-level tilt"),            LINE_MATCH_FACET_TILT,   },
130 │ };

Further down in the execute() function we find...

 374 │     /* Perform the correction. */
 375 │     if (method == LINE_MATCH_POLY) {
 376 │         if (max_degree == 0)
 377 │             linematch_do_trimmed_mean(myfield, mymask, shifts, masking, 0.0);
 378 │         else
 379 │             gwy_data_field_row_level_poly(myfield, mymask, masking, max_degree, shifts);
 380 │     }
 381 │     else if (method == LINE_MATCH_MEDIAN)
 382 │         linematch_do_trimmed_mean(myfield, mymask, shifts, masking, 0.5);
 383 │     else if (method == LINE_MATCH_MEDIAN_DIFF)
 384 │         linematch_do_trimmed_diff(myfield, mymask, shifts, masking, 0.5);
 385 │     else if (method == LINE_MATCH_MODUS)
 386 │         linematch_do_modus(myfield, mymask, shifts, masking);
 387 │     else if (method == LINE_MATCH_MATCH)
 388 │         linematch_do_match(myfield, mymask, shifts, masking);
 389 │     else if (method == LINE_MATCH_FACET_TILT)
 390 │         linematch_do_facet_tilt(myfield, mymask, shifts, masking);
 391 │     else if (method == LINE_MATCH_TRIMMED_MEAN)
 392 │         linematch_do_trimmed_mean(myfield, mymask, shifts, masking, trim_fraction);
 393 │     else if (method == LINE_MATCH_TMEAN_DIFF)
 394 │         linematch_do_trimmed_diff(myfield, mymask, shifts, masking, trim_fraction);
 395 │     else {
 396 │         g_assert_not_reached();
 397 │     }

Which points us to linematch_do_trimmed_mean() for Median (and linematch_do_trimmed_diff() for Median Differences).

 439 │ static void
 440 │ linematch_do_trimmed_mean(GwyDataField *field, GwyDataField *mask, GwyDataLine *shifts,
 441 │                           GwyMaskingType masking, gdouble trimfrac)
 442 │ {
 443 │     GwyDataLine *myshifts;
 444 │
 445 │     myshifts = gwy_data_field_find_row_shifts_trimmed_mean(field, mask, masking, trimfrac, 0);
 446 │     gwy_data_field_subtract_row_shifts(field, myshifts);
 447 │     gwy_data_line_assign(shifts, myshifts);
 448 │     g_object_unref(myshifts);
 449 │ }

So its calling gwy_data_field_find_row_shifts_trimmed_mean() which is defined in...

gwyddion/libprocess/correct.c

1579 │ /**
1580 │  * gwy_data_field_find_row_shifts_trimmed_mean:
1581 │  * @data_field: A data field.
1582 │  * @mask: Mask of values to take values into account/exclude, or %NULL for full @data_field.
1583 │  * @masking: Masking mode to use.  See the introduction for description of masking modes.
1584 │  * @trimfrac: Fraction of lowest values and highest values to discard when trimming.
1585 │  * @mincount: Minimum number of values in a row necessary for per-row calculation.  Rows which are essent
1586 │  *            completely masked are not shifted with respect to a global value.  Pass a non-positive numb
     │ se an
1587 │  *            automatic minimum count.
1588 │  *
1589 │  * Finds row shifts to misaligned row correction using trimmed row means.
1590 │  *
1591 │  * For zero @trimfrac the function calculates row means.  For @trimfrac of 1/2 or larger it calculates ro
     │ ns.
1592 │  * Values between correspond to trimmed means.
1593 │  *
1594 │  * Returns: A newly created data line containing the row shifts, for instance row means, medians or trimm
     │ s.
1595 │  *
1596 │  * Since: 2.52
1597 │  **/
1598 │ GwyDataLine*
1599 │ gwy_data_field_find_row_shifts_trimmed_mean(GwyDataField *data_field,
1600 │                                             GwyDataField *mask,
1601 │                                             GwyMaskingType masking,
1602 │                                             gdouble trimfrac,
1603 │                                             gint mincount)
1604 │ {
1605 │     GwyDataLine *shifts;
1606 │     gint xres, yres;
1607 │     const gdouble *d, *m;
1608 │     gdouble total_median;
1609 │     gdouble *sdata;
1610 │
1611 │     if (!_gwy_data_field_check_mask(data_field, &mask, &masking))
1612 │         return NULL;
1613 │     xres = data_field->xres;
1614 │     yres = data_field->yres;
1615 │
1616 │     if (mincount <= 0)
1617 │         mincount = GWY_ROUND(log(xres) + 1);
1618 │
1619 │     shifts = gwy_data_line_new(yres, data_field->yreal, FALSE);
1620 │     shifts->off = data_field->yoff;
1621 │     gwy_data_field_copy_units_to_data_line(data_field, shifts);
1622 │     total_median = gwy_data_field_area_get_median_mask(data_field, mask, masking, 0, 0, xres, yres);
1623 │
1624 │     d = data_field->data;
1625 │     m = mask ? mask->data : NULL;
1626 │     sdata = shifts->data;
1627 │
1628 │ #ifdef _OPENMP
1629 │ #pragma omp parallel if(gwy_threads_are_enabled()) default(none) \
1630 │             shared(d,m,sdata,xres,yres,masking,shifts,total_median,trimfrac,mincount)
1631 │ #endif
1632 │     {
1633 │         gdouble *buf = g_new(gdouble, xres);
1634 │         gint ifrom = gwy_omp_chunk_start(yres), ito = gwy_omp_chunk_end(yres);
1635 │         gint i, j, count;
1636 │         const gdouble *row, *mrow;
1637 │
1638 │         for (i = ifrom; i < ito; i++) {
1639 │             if (m) {
1640 │                 row = d + i*xres;
1641 │                 mrow = m + i*xres;
1642 │                 count = 0;
1643 │                 if (masking == GWY_MASK_INCLUDE) {
1644 │                     for (j = 0; j < xres; j++) {
1645 │                         if (mrow[j] > 0.0)
1646 │                             buf[count++] = row[j];
1647 │                     }
1648 │                 }
1649 │                 else {
1650 │                     for (j = 0; j < xres; j++) {
1651 │                         if (mrow[j] < 1.0)
1652 │                             buf[count++] = row[j];
1653 │                     }
1654 │                 }
1655 │                 if (count >= mincount)
1656 │                     sdata[i] = trimmed_mean_or_median(buf, count, trimfrac);
1657 │                 else
1658 │                     sdata[i] = total_median;
1659 │             }
1660 │             else {
1661 │                 gwy_assign(buf, d + i*xres, xres);
1662 │                 sdata[i] = trimmed_mean_or_median(buf, xres, trimfrac);
1663 │             }
1664 │         }
1665 │
1666 │         g_free(buf);
1667 │     }
1668 │
1669 │     zero_level_row_shifts(shifts);
1670 │     return shifts;
1671 │ }

There is some removal of masked values here but ultiamtely sdata[i] is assigned the result of either
trimmed_mean_or_median() or total_median.

We see early on in the function that total_median = gwy_data_field_area_get_median_mask(data_field, mask, masking, 0, 0, xres, yres); which is defined in libprocess/stats.[c|h] whilst trimmed_mean_or_median() is defined in
libprocess/correct.c

trimmed_mean_or_median()

1553 │ static gdouble
1554 │ trimmed_mean_or_median(gdouble *array, guint n, gdouble p)
1555 │ {
1556 │     guint nlowest = GWY_ROUND(p*n);
1557 │     guint nhighest = GWY_ROUND(p*n);
1558 │
1559 │     if (nlowest + nhighest + 1 >= n)
1560 │         return gwy_math_median(n, array);
1561 │     return gwy_math_trimmed_mean(n, array, nlowest, nhighest);
1562 │ }

gwy_data_field_area_get_median_mask()

2542 │ /**
2543 │  * gwy_data_field_area_get_median_mask:
2544 │  * @data_field: A data field.
2545 │  * @mask: Mask specifying which values to take into account/exclude, or %NULL.
2546 │  * @mode: Masking mode to use.  See the introduction for description of masking modes.
2547 │  * @col: Upper-left column coordinate.
2548 │  * @row: Upper-left row coordinate.
2549 │  * @width: Area width (number of columns).
2550 │  * @height: Area height (number of rows).
2551 │  *
2552 │  * Computes median value of a data field area.
2553 │  *
2554 │  * Returns: The median value.
2555 │  *
2556 │  * Since: 2.18
2557 │  **/
2558 │ gdouble
2559 │ gwy_data_field_area_get_median_mask(GwyDataField *dfield,
2560 │                                     GwyDataField *mask,
2561 │                                     GwyMaskingType mode,
2562 │                                     gint col, gint row,
2563 │                                     gint width, gint height)
2564 │ {
2565 │     gdouble median = 0.0;
2566 │     const gdouble *datapos, *mpos;
2567 │     gdouble *buffer;
2568 │     gint i, j, xres;
2569 │     guint n;
2570 │
2571 │     /* empty_is_noop=TRUE for compatibility */
2572 │     if (!_gwy_data_field_check_area(dfield, col, row, width, height, TRUE)
2573 │         || !_gwy_data_field_check_mask(dfield, &mask, &mode))
2574 │         return median;
2575 │
2576 │     xres = dfield->xres;
2577 │     if (mask) {
2578 │         buffer = g_new(gdouble, width*height);
2579 │         datapos = dfield->data + row*xres + col;
2580 │         mpos = mask->data + row*xres + col;
2581 │         n = 0;
2582 │         for (i = 0; i < height; i++) {
2583 │             const gdouble *drow = datapos + i*xres;
2584 │             const gdouble *mrow = mpos + i*xres;
2585 │
2586 │             if (mode == GWY_MASK_INCLUDE) {
2587 │                 for (j = 0; j < width; j++) {
2588 │                     if (mrow[j] > 0.0)
2589 │                         buffer[n++] = drow[j];
2590 │                 }
2591 │             }
2592 │             else {
2593 │                 for (j = 0; j < width; j++) {
2594 │                     if (mrow[j] <= 0.0)
2595 │                         buffer[n++] = drow[j];
2596 │                 }
2597 │             }
2598 │         }
2599 │
2600 │         if (n)
2601 │             median = gwy_math_median(n, buffer);
2602 │
2603 │         g_free(buffer);
2604 │
2605 │         return median;
2606 │     }
2607 │
2608 │     if (col == 0 && width == xres
2609 │         && row == 0 && height == dfield->yres)
2610 │         return gwy_data_field_get_median(dfield);
2611 │
2612 │     buffer = g_new(gdouble, width*height);
2613 │     datapos = dfield->data + row*xres + col;
2614 │     if (height == 1 || (col == 0 && width == xres))
2615 │         gwy_assign(buffer, datapos, width*height);
2616 │     else {
2617 │         for (i = 0; i < height; i++)
2618 │             gwy_assign(buffer + i*width, datapos + i*xres, width);
2619 │     }
2620 │     median = gwy_math_median(width*height, buffer);
2621 │     g_free(buffer);
2622 │
2623 │     return median;
2624 │ }

gwy_math_median() is defined in libgwyddion/gwymath-rank.c

265 │ /**
266 │  * gwy_math_median:
267 │  * @n: Number of items in @array.
268 │  * @array: Array of doubles.  It is shuffled by this function.
269 │  *
270 │  * Finds median of an array of values using Quick select algorithm.
271 │  *
272 │  * See gwy_math_kth_rank() for details of how the values are shuffled.
273 │  *
274 │  * Returns: The median value of @array.
275 │  **/
276 │ gdouble
277 │ gwy_math_median(gsize n, gdouble *array)
278 │ {
279 │     return gwy_math_kth_rank(n, array, n/2);
280 │ }

...which is, as noted, just the median as k is half the number of observations (did look at gwy_math_kth_rank() but
no need to dig deeper).

Conclusion

I think we can therefore use the Filters.median_flatten() function from TopoStats! 🎉