Extracting the indices of the contour

[oashour]

Let's take the sample code from the documentation

for cl in levels(contours(x,y,z))
    lvl = level(cl) # the z-value of this contour level
    for line in lines(cl)
        xs, ys = coordinates(line) # coordinates of this line segment
    end
end

I need the indices of xs and ys in x and y, i.e.in a simple world this would be:

ix = indexin(xs, x)
iy = indexin(ys, y)
this_contour = z[ix, iy]

However, xs and ys do not belong exactly to x and y, and I can't figure out how to extract the indices. The reason I need to do this is that I have an auxiliary array w that has the same shape as z, and I would like to extract the part corresponding to the contours.

Is there some way to do this? Does contour.jl return the indices somehow?

[sjkelly]

In short, there are not provisions for this in the top level API. However, the data structures used for the extraction preserve the indices as a hash:

Contour.jl/src/Contour.jl

Lines 224 to 232 in 32602e7

	# Process ambiguous cells (case 5 and 10) using
	# a bilinear interpolation of the cell-center value.
	if case == 0x05
	cells[(xi, yi)] = 0.25*sum(elts) >= h ? NWSE : NESW
	elseif case == 0x0a
	cells[(xi, yi)] = 0.25*sum(elts) >= h ? NESW : NWSE
	else
	cells[(xi, yi)] = edge_LUT[case]
	end

The key would be to modify this function here:

Contour.jl/src/Contour.jl

Lines 242 to 264 in 32602e7

	function chase!(cells, curve, x, y, z, h, start, entry_edge, xi_range, yi_range, ::Type{VT}) where VT
	ind = start
	# When the contour loops back to the starting cell, it is possible
	# for it to not intersect with itself. This happens if the starting
	# cell contains a saddle-point. So a loop is only closed if the
	# contour returns to the starting edge of the starting cell
	loopback_edge = entry_edge
	@inbounds while true
	exit_edge = get_next_edge!(cells, ind, entry_edge)
	push!(curve, interpolate(x, y, z, h, ind, exit_edge, VT))
	ind, entry_edge = advance_edge(ind, exit_edge)
	!((ind[1], ind[2], entry_edge) != (start[1], start[2], loopback_edge) &&
	ind[2] ∈ yi_range && ind[1] ∈ xi_range) && break
	end
	return ind
	end

The ind parameter, I believe is exactly what you would want. You can just push that to the array as a tuple or whatever is suitable.

It should be possible to add as a feature here. Unfortunately I am quite busy with other projects so I cannot do it myself. I am happy to advise where I can and review any PRs. Though we should probably create a non-API breaking way to do this. There is an example of adding a new algorithm in #56 that is non-breaking (though that algorithm is broke, I think, it has been a while :).

[oashour]

Thank you for the detailed response! Is it reasonable to modify ContourLevel{T} to add another field indices, i.e.

struct ContourLevel{T}
    level::T
    lines::Vector{Curve2{T}}
    indices::Vector{Int64}
end

with a function

indices(cl::ContourLevel) = cl.indices

To mimic how lines works? Or is this considered breaking? I'll have some time to work on this over the weekend so hopefully I can come up with something that doesn't break the API, it seems straightforward enough I hope.

P.S. I'm very new to collaborating with other people on stuff so bear with me.

[sjkelly]

I would consider modifying ContourLevel to be breaking (this is user facing). We could add a new ContourLevelIndexed with what you proposed. However since chase! is internal we can modify that with extra parameters to track if we want to push the indices or not.

[oashour]

Sorry for the trivial questions but this is a learning moment for me. Why is adding a field to a user facing structure considered breaking? I can understand renaming or removing fields of course, but what's wrong with adding fields?

As for chase!, are you suggesting we modify the definition of chase from

Contour.jl/src/Contour.jl

Line 242 in 32602e7

function chase!(cells, curve, x, y, z, h, start, entry_edge, xi_range, yi_range, ::Type{VT}) where VT

to add two new parameters

function chase!(cells, curve, x, y, z, h, start, entry_edge, xi_range, yi_range, save_ind, ind_arr, ::Type{VT}) where VT 

where save_ind is a boolean flag and ind_arr is an array to store the indices. Then we can simply modify this part:

Contour.jl/src/Contour.jl

Line 257 in 32602e7

ind, entry_edge = advance_edge(ind, exit_edge)

to something like

        ind, entry_edge = advance_edge(ind, exit_edge)
	if save_ind
		push!(ind_arr, ind)
	end

Is this what you're suggesting? I'm not sure how to make ind_arr user accessible, and defining a second type, ContourLevelIndexed seems cumbersome.

I don't have access to a computer with Julia at the moment but I'll try to play around with this as soon as I do. I'm just stuck theory crafting till then.

[sjkelly]

I appreciate the questions, so no worries. Adding a field may seem non-breaking in the Julia ecosystem. There may be macros or memory assumptions about the layout of data that may make it breaking. Similarly we are going to use more memory in this approach and the algo will likely be slower, so it is a good idea to keep it as a different algorithm. This is me being over conservative to avoid any regressions or upstream issues.

The approach seems about right. If we setup a different type assuming the entry point will be something like contour(x,y,x,ContourIndexed()). We can duck type and save_ind is not needed. E.g. check the type of ind_arr if it is not Nothing means we can save. This means that branch should be elided for the default algorithm during compilation.

I may have some time next weekend to work on this. Also feel free to ping me on the Julia slack if you are there.