# Interpolation into array rather than from array

**URL:** <https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538>\
**Category:** General Usage\
**Tags:** interpolations\
**Created:** [May 4, 2021, 5:50pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538 "2021-05-04T17:50:19Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [May 4, 2021, 5:50pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/1 "2021-05-04T17:50:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [May 5, 2021, 6:58am UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/2 "2021-05-05T06:58:52Z")

</div>

Note that `getindex` is deprecated for non-integer values in `Interpolations.jl`. The function call usage is preferred.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 5, 2021, 7:21am UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/3 "2021-05-05T07:21:54Z")

</div>

> [@RainerHeintzmann](#):
>
> I would like to do interpolation writing **into** an array rather than interpolation **from** an array.

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?

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [May 5, 2021, 9:29am UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/4 "2021-05-05T09:29:43Z")

</div>

Here is a dummy that maybe shows roughly what I mean. Don’t take it as a proper implementation:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 5, 2021, 10:41am UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/5 "2021-05-05T10:41:05Z")

</div>

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`)

```julia
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])

```

![Interpolations_linear_with_extrapolation](https://global.discourse-cdn.com/julialang/original/3X/0/e/0e2d1d0bda3cf8ef9c9a7a5790983f934f7429e5.png)

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [May 5, 2021, 11:33am UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/6 "2021-05-05T11:33:50Z")

</div>

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:

```julia
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")

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/3/630b52cd0c1c0a6c348c4cab33303db6de51e1dc.png)

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 5, 2021, 11:56am UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/7 "2021-05-05T11:56:26Z")

</div>

Sorry, this looks too complicated for me…  
_I’m throwing in the towel…_  
_Je jette l’éponge…_  
_Ich werfe das Handtuch…_

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [May 5, 2021, 1:04pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/8 "2021-05-05T13:04:55Z")

</div>

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](https://www.researchgate.net/publication/351347162_Resolution_Enhancement_of_Biological_Light_Microscopic_Data)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/e/fedff6dce113b229e2e49f811a707604c8232232.png)

Note that this is one step beyond of what I was asking for in the examples above.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 5, 2021, 1:31pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/9 "2021-05-05T13:31:40Z")

</div>

> [@RainerHeintzmann](#):
>
> 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:

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.

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [May 5, 2021, 3:40pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/10 "2021-05-05T15:40:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 5, 2021, 10:12pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/11 "2021-05-05T22:12:08Z")

</div>

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.”

```julia
plot(heatmap(reshape(a,1,:),c=:greys,title="a",),heatmap(reshape(b,1,:),c=:greys,title="b"),size=(1000,250), yticks=:false)

```

 ![heatmap_pixel_interpolation](https://global.discourse-cdn.com/julialang/original/3X/1/c/1c6e6f9509e79e67a336b038710bab4128d64310.png)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 9, 2021, 1:50pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/12 "2021-05-09T13:50:04Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Jake](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jake/32/46007_2.png) [@Jake](https://discourse.julialang.org/u/Jake)\
**Post date:** [May 10, 2021, 9:48pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/13 "2021-05-10T21:48:41Z")

</div>

> [@stevengj](#):
>
> 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 recently implemented the sinc interpolation scheme mentioned in the [paper](http://toc.proceedings.com/00102webtoc.pdf) on page 583. Do you have a more detailed reference of the technique you describe?

---

<div class="post-metadata">

**Author:** ![RainerHeintzmann](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rainerheintzmann/32/19726_2.png) [@RainerHeintzmann](https://discourse.julialang.org/u/RainerHeintzmann)\
**Post date:** [May 11, 2021, 12:01pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/14 "2021-05-11T12:01:15Z")

</div>

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](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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 11, 2021, 12:45pm UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/15 "2021-05-11T12:45:15Z")

</div>

> [@RainerHeintzmann](#):
>
> As for “restriction” I found surprisingly little on the internet.

I think that term is used mostly in the context of multigrid methods.

---

<div class="post-metadata">

**Author:** ![pxshen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pxshen/32/38776_2.png) [@pxshen](https://discourse.julialang.org/u/pxshen)\
**Post date:** [November 26, 2022, 4:03am UTC](https://discourse.julialang.org/t/interpolation-into-array-rather-than-from-array/60538/16 "2022-11-26T04:03:57Z")

</div>

“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!`](https://aced-differentiate.github.io/EquivariantOperators.jl/guide/#place) implements it.