Contour Smoothing Option: Dual Marching Squares

[quiqueg]

Summary

A new "dual linear contour smoothing" option offers improved contour quality with minimal additional performance cost.

Background

Dual Marching Squares is a simple extension of the original (primal) Marching Squares algorithm.

Primal Marching Squares contour smoothing entails linear interpolation of isoline points along grid cell edges.

Dual Marching Squares still performs linear interpolation along grid cell edges, but it uses the dual of the isocontour that was generated from that linear interpolation step. In other words, compared to primal Marching Squares, Dual Marching Squares adds an additional smoothing step in which the mid-points of all segments of a given isoline are joined together; the result is used as the isoline.

Summary of the benefits of Dual Marching Squares

Compared to primal Marching Squares, Dual Marching Squares:

• Produces isocontours that better approximate curved isolines, particularly in areas of high curvature, i.e. "sharp features".
• Resolves ambiguity at saddle points, i.e. grid cells in which all four edges are bipolar edges.
  • In primal Marching Squares, this ambiguity occurs when two opposite corners of a grid cell are both above the isovalue, while the remaining 2 opposite corners are below the isovalue. The orientation of the saddle is ambiguous:
  

API changes

In this Pull Request, I have added a new "linearDual" option for the smooth argument to the contour smooth function:

# contours.smooth([smooth])	Smoothing method
false	No smoothing
true or "linear" (Default)	Linear smoothing (Primal Marching Squares)
"linearDual"*	Dual linear smoothing (Dual Marching Squares)

* New option

Implementation details

The smoothLinearDual smoothing function has been added to contours.js alongside the existing smoothLinear function. It can be enabled by invoking contours.smooth("linearDual").

First it calls smoothLinear and then modifies the returned isoring by shifting each vertex to the midpoint of the next isoring segment. Effectively, it creates the final isoring by joining the isoring midpoints iteratively, in a single pass.

According to the current API behavior, any truthy value for smooth will enable linear smoothing, whereas passing false will disable smoothing. In order to maintain backward compatibility, my changes simply check if smooth is equal to "linearDual". If yes, then the dual linear smoothing method is used; if no, then fallback to the original behavior. I have enhanced the documentation to suggest that "linear" can also be passed and will be treated the same as true.

We could consider adding more smoothing options in the future. Some examples could be logarithmic interpolation (both primal and dual), or bilinear interpolation (note that dual linear interpolation is not the same as bilinear interpolation). However, I have not explored potential use cases nor the effectiveness of the contour smoothing that those options might produce.

Go to demo

The demo page linked above demonstrates the effects of dual linear contour smoothing, side-by-side with no smoothing and linear smoothing.

Controls are provided to simulate different data densities, as well as different numbers of contour lines. Finally, the demo can be animated to show isolines at all possible isovalues in the domain.

Select snapshots

Medium-density data

With medium-density data, both linear and dual linear contour smoothing noticeably improve the contour quality, but the "no smoothing" option can still be suitable for some applications which do not require precise contours. Dual linear contour smoothing shows a subtle improvement over linear smoothing.

Low/Medium-density data

With low/medium-density data, both linear and dual linear contour smoothing significantly improve the contour quality. Dual linear contour smoothing fixes some specific issues (see Notable improvements below) that are otherwise present with linear smoothing.

Low-density data

With low-density data, both linear and dual linear contour smoothing significantly improve the contour quality. Dual linear contour smoothing fixes some specific issues (see Notable improvements below) that are otherwise present with linear smoothing.

Notable improvements

Reduced grid artifacts

In the original Primal Marching Squares algorithm with linear smoothing, all contour line points fall along grid lines, i.e. they cannot fall within the grid cells. This can result in the grid lines surfacing as artifacts in the final contour drawing: see the middle image below where several contour lines appear to "fold" over an invisible line.

Dual Marching Squares fixes this issue by moving contour line points to inside the grid cells. See the right-most image below, where the "invisible line" grid artifacts are no longer present:

Improved handling of saddles

A "saddle" occurs when two isorings with the same value intersect the same grid square. With Primal Marching Squares, this situation can result in two isorings converging at a point, as demonstrated in the middle image below.

Dual Marching Squares fixes this issue by connecting each pair of isoring segments at their mid-points, which has the effect of bringing the isoring closer to its center and thus away from the former saddle point: see the right-most image below.

Performance implications

The performance impact of both linear and dual linear smoothing is minimal. See the table below for the added runtime cost associated with each smoothing option:

		Low-density data		High-density data
Smoothing method		Contour quality	Performance cost	Contour quality	Performance cost
None	false	Poor	+0%	Good	+0%
Linear (Default)	true or "linear"	Good	+2.40%	Best	+1.55%
Dual linear	"linearDual"	Best	+2.70%	Best	+2.00%

In all cases, an overwhelming majority of the algorithm's time is spent on the stitching, isoring processing, and contour processing logic (CPU self time, without smoothing):

	Stitching	Isoring processing	Contour processing	Other
Low-density data	64.36%	30.96%	0.70%	3.98%
High-density data	29.15%	33.77%	33.53%	3.55%

Low-density data

The data source used for these performance tests is the topography of Maungawhau (the classic volcano dataset from R), which is used for several d3-contour examples such as https://observablehq.com/@d3/volcano-contours

No smoothing

Linear smoothing

Dual linear smoothing

High-density data

The data source used for these performance tests is the same GeoTIFF data used here: https://observablehq.com/@d3/geotiff-contours-ii

No smoothing

Linear smoothing

Dual linear smoothing

[Fil]

Nice work! 👏

As a minor remark I don't think we need "VS Code settings and build task". Also, I would simplify the option's name as "dual" (contours.smooth("dual")).

[Fil]

Suggested change

	- `true` or `"linear"` (Default): Linear interpolation smoothing, as described in the original [Marching Squares algorithm](https://en.wikipedia.org/wiki/Marching_squares).
	- `true` or `"linear"` (default): Linear interpolation smoothing, as described in the original [Marching Squares algorithm](https://en.wikipedia.org/wiki/Marching_squares).

[Fil]

Suggested change

	Additionally, the density of the source data impacts contour smoothness: the differences between the smoothing methods are more noticeable with low-density data than they are with high-density data. With _very_ high-density data, there is no discernible quality difference between the three methods.
	Additionally, the density of the source data impacts contour smoothness: the differences between the smoothing methods are more noticeable with low-density data than they are with medium-density data. With high-density data, there is no discernible quality difference between the three methods.

[Fil]

remove (this is interesting in the PR but too overwhelming for the README)

[Fil]

Suggested change

	return arguments.length ? (smooth = _ === 'linearDual' ? smoothLinearDual : _ ? smoothLinear : noop, contours) : smooth === smoothLinear || smooth === smoothLinearDual;
	return arguments.length ? (smooth = _ === 'dual' ? smoothLinearDual : _ ? smoothLinear : noop, contours) : smooth === smoothLinear || smooth === smoothLinearDual;

[Fil]

Suggested change

	tape("contours.smooth('linearDual')(values) returns the expected result for an empty polygon", function(test) {
	tape("contours.smooth('dual')(values) returns the expected result for an empty polygon", function(test) {

[Fil]

Suggested change

	tape("contours.smooth('linearDual')(values) returns the expected result for a simple polygon", function(test) {
	tape("contours.smooth('dual')(values) returns the expected result for a simple polygon", function(test) {

[Fil]

Suggested change

	var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
	var contours = d3.contours().smooth('dual').size([10, 10]).thresholds([0.5]);

[Fil]

Suggested change

	tape("contours.smooth('linearDual')(values) returns the expected result for polygon in the corner", function(test) {
	var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
	tape("contours.smooth('dual')(values) returns the expected result for polygon in the corner", function(test) {
	var contours = d3.contours().smooth('dual').size([10, 10]).thresholds([0.5]);

[Fil]

Suggested change

	tape("contours.smooth('linearDual')(values) returns the expected result for a complex polygon", function(test) {
	var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
	tape("contours.smooth('dual')(values) returns the expected result for a complex polygon", function(test) {
	var contours = d3.contours().smooth('dual').size([10, 10]).thresholds([0.5]);

[Fil]

Suggested change

	tape("contours.smooth('linearDual')(values) returns the expected result for a multipolygon with holes", function(test) {
	var contours = d3.contours().smooth('linearDual').size([10, 10]).thresholds([0.5]);
	tape("contours.smooth('dual')(values) returns the expected result for a multipolygon with holes", function(test) {
	var contours = d3.contours().smooth('dual').size([10, 10]).thresholds([0.5]);

[Fil]

The algorithm in this PR (apply marching squares then smooth) doesn't seem to match the description in https://www.boristhebrave.com/2018/04/15/dual-contouring-tutorial/ ? Is this expected?

Bilinear interpolation is addressed in #46

[BorisTheBrave]

Yes, this is not dual contouring, it is a different technique. Dual contouring is vastly more common, so I would recommend against labelling this option "dual" as it may cause confusion.

I would also disagree with your analysis.

Produces isocontours that better approximate curved isolines, particularly in areas of high curvature, i.e. "sharp features".
Resolves ambiguity at saddle points, i.e. grid cells in which all four edges are bipolar edges.

Is this actually true? These are facts about dual contouring, but not what you've done here. As you are smoothing the existing contours, you cannot be preserving sharp features better than it, and you are stuck with the existing implementations handling at saddle points.