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 (ImageTransformations.jl/warp.jl at master Â· JuliaImages/ImageTransformations.jl Â· GitHub) 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.