I would like to do interpolation writing into an array rather than interpolation from an array. I looked around a bit and found Interpolations.jl. Yet it seems like it only supports to interpolate on read (getindex)instead of on write (setindex). One could argue that you can invert the projection function and then obtain the destination array by interpolation on read. Yet it is not always easy or even possible to invert the mapping function. Does anyone know of a package that also supports interpolation on write?
This means, that when writing to a subpixel coordinate mymap(CartesianIndex) it would write (i.e. add) to a couple of pixels with appropriate weights. I understand, that it may be necessary to construct and divide by a sum-of-weights image to compensate for the weights not summing to one.
Any help is appreciated.
Note that getindex is deprecated for non-integer values in Interpolations.jl. The function call usage is preferred.
I am not sure that’s a well-defined problem in the context of interpolation. Consider gridpoints x_1, x_2, x_3, with values y_1 etc, and a linear interpolation.
Now supposes you want to “write” (x, y), with x_1 < x < x_2. What would be the new values of y_1, y_2, y_3? Or do you mean to insert a new point?
Here is a dummy that maybe shows roughly what I mean. Don’t take it as a proper implementation:
julia> a=[0.0,0.0,1.0,1.0,0.0,0.0]; b=a.*0;
julia> dest_pos = [1.2,2.3,3.3,4.6,5.5,6.4];
julia> w_l = dest_pos .- floor.(dest_pos); w_r = 1 .- w_l;
julia> d_l = floor.(Int,dest_pos); d_r = d_l[1:end-1] .+1;
julia> b[d_l] .= w_l[d_l] .* a[1:end];
julia> b[d_r] .+= w_r[d_r] .* a[1:end-1];
julia> b
6-element Vector{Float64}:
0.0
0.0
0.2999999999999998
1.0
0.5
0.0
Could not make sense of your example…
Is below what you want? (i.e., given a at integer positions 1:6, compute c at positions dest_pos)
using Interpolations, Plots; gr()
a = [0., 0., 1., 1., 0., 0.]; n = 1:length(a)
intlin = LinearInterpolation(n, a, extrapolation_bc=Line());
dest_pos = [1.2, 2.3, 3.3, 4.6, 5.5, 6.4];
c = intlin(dest_pos)
plot(n, a); scatter!([n dest_pos],[a c], label=["a" "c"], c=[:red :green])
Thanks for this nice graphical visualization. You show a nice example of what I would call Interpolation on the source data. What you call dest_pos in your script really are “source positions”, not destination positions. The difference becomes very clear, if we zoom a lot. Let me demonstrate what I mean:
using Interpolations, Plots; gr()
a=[0.0,0.0,1.0,1.0,0.0,0.0];
dest_pos = [10.2,20.3,30.3,40.6,50.5,60.4];
b = zeros(70);
w_r = dest_pos .- floor.(dest_pos); w_l = 1 .- w_r;
d_l = floor.(Int,dest_pos); d_r = d_l .+1;
b[d_l] .= w_l .* a;
b[d_r] .+= w_r .* a;
plot(dest_pos, a, label="on source")
scatter!(b, label="on destination")
The result of your routine would all reside on the “linear source” line, whereas I want the results to be as shown. Of course one can argue how much sense such a zoom makes but the same arguments can be made if you zoom out a lot and use your routine. The you would be “missing” a lot of data rather than averaging it, which is what my routine would do.
Sorry, this looks too complicated for me…
I’m throwing in the towel…
Je jette l’éponge…
Ich werfe das Handtuch…
Its not that complicated. The description on page 123 of the following document has a nice picture which shows how to interpolate considering BOTH, source and destination voxels:
https://www.researchgate.net/publication/351347162_Resolution_Enhancement_of_Biological_Light_Microscopic_Data
Note that this is one step beyond of what I was asking for in the examples above.
I think know what you mean, though in practice I would call this operation a type of “smoothing” rather than “interpolation” per se. (Effectively you are taking a source function, convolving it with a smoothing kernel, and then sampling the smoothed function on a grid.) A common way to do this is to multiply by the transpose of the interpolation matrix (that takes grid points to source points) — this is also called a restriction operation.
I don’t recall offhand having seen a package for this, however.
I am not sure that I understand what we are talking about. The “both” or the “destination” interpolation/restriction? In continuous space, I can see the convolution argument, but in practice this continuous space is never calculated. If we talk about the “destination-only” version as indicated at the top, the thinking would be to have in addition to getindex also a setindex routine which adds onto the destination grid. Now that getindex with float arguments is deprecated, this would be called different and packaged into a routine, but it can still be done.
A combine algorithm for the shown example would be an interpolation with a trapezoidal shaped interpolation kernel for regular grids. For irregular ordered grids this kernel would ideally adapt to the local spacing. The big advantage: Not loosing SNR when zooming out (e.g with noisy data), approaching binning for very small zooms and approaching (e.g. linear) interpolation when zooming in.
The problem seems to be about interpolating pixels or voxels, not points. Maybe it would have been better to have started from the start.
Nevertheless, the Julia code provided still looks confusing because we do not know where the “spatial” coordinates nor the pixel intensities to interpolate are?
It would be better if you could provide a MWE that is easier to understand.
For instance displaying the pixels as heatmaps may allow a better understanding, including of your interpolation constraints (“the integral intensity of an image should be kept constant, etc.”
plot(heatmap(reshape(a,1,:),c=:greys,title="a",),heatmap(reshape(b,1,:),c=:greys,title="b"),size=(1000,250), yticks=:false)
In particular, suppose you have a vector d of samples and you form some interpolation (e.g. billinear, spline) to obtain a vector s of interpolated points, represented by a linear operation s = Md. It sounds like what you want is the transposed “restriction” operation d = M^T s.
I recently implemented the sinc interpolation scheme mentioned in the paper on page 583. Do you have a more detailed reference of the technique you describe?
Is there a link to the paper describing the interpolation method, rather that the link to the index of the journal?
As for “restriction” I found surprisingly little on the internet. But I saw that in the warp function of the Interpolations.jl toolbox (https://github.com/JuliaImages/ImageTransformations.jl/blob/master/src/warp.jl) it is said: “This approach is known as backward mode warping.”, which may imply that “forward mode warping” also exists? But it does not seem to be implemented yet.
I think that term is used mostly in the context of multigrid methods.
“Interpolation into array” is known as the particle mesh method, meaning sum of the values written onto the array scaled by the discretization volume is equal to the discrete value. (Though unsure if that’s what the OP wants.) EquivariantOperators.place! implements it.