# Help implement Kalman filter

**URL:** <https://discourse.julialang.org/t/help-implement-kalman-filter/8718>\
**Category:** General Usage\
**Tags:** question, kalman\
**Created:** [January 31, 2018, 2:45pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718 "2018-01-31T14:45:36Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [January 31, 2018, 2:45pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/1 "2018-01-31T14:45:36Z")

</div>

Hi all!  
I have a specific problem that I think lends itself to a Kalman filter. Since I’ve never used a Kalman filter before and I’m not sure how exactly I could apply it on my data, I thought I’d reach-out here.

## My data

I have time-lapse images of roots growing. As they grow they emit a bit of light from their tip. I want to track the trajectory of the root tips.

 ![1](https://global.discourse-cdn.com/julialang/original/3X/4/b/4b33f13555ff5a44e3f110e79ca8d34130bed626.gif)

The image of the root tip looks like a blurry point source – the distribution of the light intensities sometimes looks almost like a gaussian (albeit not in the below example).

![1](https://global.discourse-cdn.com/julialang/original/3X/1/5/15ddb19d508f349c2349e51a97f48f8b28f5c0ac.jpeg)

I can get an ok estimation of the location of the beginning of a root-tip so I know where to start in the image. I also have good guesses about the tip’s:

1. speed
2. angular speed (how quick it can change directions)
3. range of possible directions (the root doesn’t grow “up”, it only grows with gravitation)
4. intensity change (how the tip brightness degrades with time)

## How

I read [this](http://www.bzarg.com/p/how-a-kalman-filter-works-in-pictures/) excellent explanation about Kalman filters, I checked out the `Kalman`, `StateSpace`, `DataAssim`, and `StateSpaceRoutines` packages and I thought that I should be able to:

1. Start with an initial tip location
2. Guess the location of the tip in the next frame
3. Extract a large enough region of interest (ROI) from that next frame centered on the guess
4. Use the light distribution in that ROI to feed the Kalman filter with an estimation of the new location

## Implementation

I’m a bit lost… Does anyone here have a good understanding of how I could use one of the available tools in Julia to accomplish this?

Thanks!

## Example data

This will generate some data to play with:

```julia-auto
sz = (100,100)
nframe = 10
imgs = rand(sz..., nframe);
t = 1:nframe
x = 1.1*t + sz[2]/2 + rand(nframe)
y = 7.8*t + 1 + rand(nframe)
i = [CartesianIndex(round(Int, y[frame]), round(Int, x[frame]), frame) for frame in 1:nframe]
imgs[i] = 5e3
using ImageFiltering
for frame in 1:nframe
    imgs[:,:,frame] = imfilter(imgs[:,:,frame], Kernel.gaussian(2))
end

```

---

<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:** [January 31, 2018, 3:11pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/2 "2018-01-31T15:11:41Z")

</div>

Your data is not in state space form (that I could easily imagine as applicable for the Kalman filter). Eg your state could be (x, y) coordinates of the root tips, and then perhaps you could use a random walk or an I(1) process (“velocity” changes randomly), but those coordinates would need to be extracted first from the pixels. AFAICT that is orthogonal to the application of the Kalman filter. (BTW, amazing stuff).

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [January 31, 2018, 3:27pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/3 "2018-01-31T15:27:39Z")

</div>

I’ve successfully followed the tip by just finding the closest and brightest pixel below the previous one. So I could generate an approximation for where the tips might be.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [January 31, 2018, 7:44pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/4 "2018-01-31T19:44:14Z")

</div>

Yes, amazing pictures. To have something more robust, instead of taking the brightest pixel, you can do the following: Given the current position, you can apply a window or kernel and then compute the average of the light location (mean of the light distribution) within the window.

```julia
ij(x,y) = round(Int, y), round(Int, x)

# initial guess x y
x = ...
y = ...
M = imgs[:,:,1] 

kernel(x, y, σ) = 1/(2*π*σ^2)*exp(-hypot(x,y)^2/(2*σ^2))

function track(M, x, y, σ = 2, n = 7)
    si = 0.0
    sj = 0.0
    s = 0.0
    i0, j0 = ij(x, y)
    for i in i0-n:i0+n
        for j in j0-n:j0+n
            im = mod1(i, size(M,1))
            jm = mod1(j, size(M,2))
            
            si += i*M[im,jm]*kernel(y - i, x - j, σ)
            sj += j*M[im,jm]*kernel(y - i, x - j, σ)
            s += M[im,jm]*kernel(y - i, x - j, σ)
        end
    end
    sj/s, si/s
end

xx = zeros(nframe)
yy = zeros(nframe)

for frame in 1:nframe
    x, y = track(imgs[:, :, frame], x, y, 10, 16)
    xx[frame] = x
    yy[frame] = y
end

[xx yy]

```

In the example after playing a bit with the window parameters, say sigma=10 and n=16, I get

```julia
 x xtrue y ytrue
51.9744 51.8124 9.98118 9.76172
 52.9121 53.0537 16.3462 17.1948 
 53.8929 53.8186 24.1064 24.813  
 54.88 54.799 32.098 33.1901 
 55.8867 56.2031 40.0934 40.5572 
 57.7943 57.5816 47.2083 48.399  
 58.8921 58.6591 55.0959 55.9413 
 58.9889 58.9721 63.0704 64.3214 
 59.8979 60.299 71.0971 71.7049 
 61.7865 61.6383 79.0698 79.9845

```

That seems to work nicely. You can also use `track` to obtain a velocity vector.

---

<div class="post-metadata">

**Author:** ![oatlzzvztd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oatlzzvztd/32/3285_2.png) [@oatlzzvztd](https://discourse.julialang.org/u/oatlzzvztd)\
**Post date:** [January 31, 2018, 8:32pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/5 "2018-01-31T20:32:27Z")

</div>

This video looks amazing! Is there any hope that you could publish the data set somewhere in the open so that others could play around with it?

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [February 1, 2018, 9:51am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/6 "2018-02-01T09:51:10Z")

</div>

Thanks!

Cool, but could you please explain what the function `ij` is or does (my only guess was `round(Int, x)`, but that doesn’t seem to work)?

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 1, 2018, 9:53am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/7 "2018-02-01T09:53:58Z")

</div>

Ah, sorry, You had it almost.

```julia
ij(x,y) = round(Int, y), round(Int, x)

```

I’ll edit.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [February 1, 2018, 9:58am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/8 "2018-02-01T09:58:37Z")

</div>

Wow! Awesome!  
So what do you mean by `track` in:

> [@mschauer](#):
>
> You can also use track to obtain a velocity vector.

?

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 1, 2018, 10:01am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/9 "2018-02-01T10:01:50Z")

</div>

`track` is just the function name. If you take

```julia
xnew, ynew = track(imgs[:, :, frame], xpred, ypred, 10, 16)
dx = xcur - xnew
dy = ycur - ynew

```

you can use that in filtering.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [February 1, 2018, 10:34am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/10 "2018-02-01T10:34:58Z")

</div>

Right, right.

So, I played around with it, and it seems like if I increase the noise by one order of magnitude:

```julia
imgs[i] = 5e2

```

then it really breaks down. I might need to play with the window size and kernel variance…  
A similar solution I’ve used before is to calculate the weighted mean of the window, using the pixel brightness as the weights. Not sure how your solution really differs from that…

While this is cool, I can’t help thinking if there is a way to use online Kalman filtering here…

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 1, 2018, 10:39am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/11 "2018-02-01T10:39:59Z")

</div>

As @Tamas_Papp said: Kalman-filtering is mostly orthogonal, if you use a function like `track` to obtain the noisy position signal, you can use Kalman filtering to update the position `x`, `y` taking a linear model for the dynamics of the root location and the current position into account.

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [February 1, 2018, 10:41am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/12 "2018-02-01T10:41:19Z")

</div>

Ok, so I understand I need to have a step between the timelapse images and the kalman filter, namely your track function or something similar that spits out an estimation of the coordinates. I can then use similar estimators that spit out the speed, angular speed, direction, and intensity change, but how do I (if at all) feed all that into a kalman filter?

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [February 1, 2018, 11:08am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/13 "2018-02-01T11:08:43Z")

</div>

Please have at it, and please let me know if you find some robust way of getting the tip’s trajectories:  
[dropbox link](https://www.dropbox.com/s/97482fgvov7nrbh/example.zip?dl=0) to a zip file with 25 TIF UInt16 images (2 MB each), temporally sorted (named 1, 2, 3, …, 25). As a starter, there is a tip in image `1.TIF` at column 798 and row 636. Not to be all negative and all, but keep in mind that while you could easily follow it by detecting the brightest point below it in the consecutive image, it would be better if we could find some more robust and accurate way of doing this.  
Thanks in advance!

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [February 1, 2018, 4:28pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/14 "2018-02-01T16:28:17Z")

</div>

> [@yakir12](#):
>
> find some robust way of getting the tip’s trajectories

I think the following is a pretty decent way of getting an approximation of the tip’s locations:

```julia
W = 10
cols = [col for row in -W:W for col in -W:W]
rows = [row for row in -W:W for col in -W:W]
win = CartesianIndex.(rows, cols)
frame = img[p .+ win]
w = Float64.(frame)
w -= mean(w)
clamp!(w, 0, 1)
yx = Float64.([rows cols])
μ = mean(yx, StatsBase.Weights(w), 1)
p = CartesianIndex(round.(Int, μ)...) + p

```

So now I have a list of estimated locations. I should therefore be able to use a Kalman filter on that list in order to improve on these estimates. Since this next step is separate from acquiring the tip estimates, I’ll start a new topic just for _that_.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [February 1, 2018, 4:46pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/15 "2018-02-01T16:46:34Z")

</div>

That is about what I also found working on your example data. I applied my solution to the normalized frame, but normalizing and thresholding the window works even better I suppose.

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

---

<div class="post-metadata">

**Author:** ![yakir12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yakir12/32/297_2.png) [@yakir12](https://discourse.julialang.org/u/yakir12)\
**Post date:** [February 1, 2018, 6:21pm UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/16 "2018-02-01T18:21:12Z")

</div>

Cool. Well, now all that remains, unless I misunderstood something, is to refine that with the Kalman filter (see [here](https://discourse.julialang.org/t/how-to-kalman-filter/8745)).

---

<div class="post-metadata">

**Author:** ![mzaffalon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mzaffalon/32/214168_2.png) [@mzaffalon](https://discourse.julialang.org/u/mzaffalon)\
**Post date:** [February 3, 2018, 7:12am UTC](https://discourse.julialang.org/t/help-implement-kalman-filter/8718/17 "2018-02-03T07:12:27Z")

</div>

[Offtopic]: This looks impressive! Can you please post a link (or two) to the science behind the light emission?