Unexpected output from FeatureRegister

[NeilYager]

Hello,

I can't get to the output of FeatureRegister to match my expectations. Here is an example:

import numpy as np
import scipy.ndimage as ndimage

from register.models import model
from register.samplers import sampler
from register import register

im1 = np.zeros((10, 10))
im1[2, 2] = 1
features1 = {}
features1['points'] = {}
features1['points'][1] = [2, 2]
im2 = np.zeros((10, 10))
im2[4, 5] = 1
features2 = {}
features2['points'] = {}
features2['points'][1] = [4, 5]

reg1 = register.RegisterData(im1, features=features1)
reg2 = register.RegisterData(im2, features=features2)

feature = register.FeatureRegister(model=model.Shift,
                   sampler=sampler.CubicConvolution)

p, warp, imWarped, error = feature.register(reg1, reg2)
print im1
print im2
print imWarped

The model finds the correct offset. However, imWarped matches neither im1 nor im2. Is this a bug, or am I doing/expecting something wrong?

[riaanvddool]

I think the problem might be that you don't have enough features. If you try the examples/nonlinreg_image_features.py script you should see a working example.

[NeilYager]

I don't think that's the problem here. This is a bare bones example, but I first noticed this in my data that has dozens of features.

[riaanvddool]

I can only comment on the example code you provided. Also note that it is usually a good idea to define plenty features on the borders of the two images, otherwise you will get gross distortion towards the edges.

[riaanvddool]

I made the following changes to your example:

im1 = np.zeros((10, 10)) im1[2, 2] = 1 features1 = {} features1['points'] = {} features1['points'][1] = [2, 2] features1['points'][2] = [0, 0] features1['points'][3] = [0, 9] features1['points'][4] = [9, 9] features1['points'][5] = [9, 0] im2 = np.zeros((10, 10)) im2[7, 7] = 1 features2 = {} features2['points'] = {} features2['points'][1] = [7, 7] features2['points'][2] = [0, 0] features2['points'][3] = [0, 9] features2['points'][4] = [9, 9] features2['points'][5] = [9, 0]

The answer looks plausible I think.

[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 1. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

[riaanvddool]

On the other hand, one would expect a shift model to need only one feature for a solution, I agree.

[NeilYager]

I'm a little confused. In my original example, the model correctly fits p as [ 2. 3.], which is the offset between the two features. There are 2 features and two unknown parameters, so I don't see how additional features would help for the Shift model (I can certainly see why they would help for nonlinear deformations). I still don't understand why warpedImage seems inconsistent with p.

[NeilYager]

Ah - I see you came to a similar conclusion.

[riaanvddool]

I am trying to determine now if it is a indexing issue: im2[4, 5] = 1 # ie y = 4, x = 5 vs features2['points'][1] = [4, 5] # which might be [x,y], ie should be [5,4] ?

[riaanvddool]

Nathan wrote this code, so I am a bit rusty...

[riaanvddool]

I think the indexing was the issue. When I swop the x and y around for the features I get the expected result. I then saw that for the bare-bones example, cubic-convolution was messing around the results. So I changed to nearest neighbor sampler and things look right. See my code below:

import numpy as np
import scipy.ndimage as ndimage

from register.models import model
from register.samplers import sampler
from register import register

im1 = np.zeros((10, 10)).astype(np.double)
im1[1, 1] = 1
im1[2, 2] = 1
im1[3, 3] = 1
im1[3, 1] = 1
im1[1, 3] = 1
im1[2, 3] = 1
im1[3, 2] = 1
im1[2, 1] = 1
im1[1, 2] = 1
features1 = {}
features1['points'] = {}
features1['points'][1] = [2, 2]
im2 = np.zeros((10, 10)).astype(np.double)
im2[4, 5] = 1
#im2[5, 5] = 1
features2 = {}
features2['points'] = {}
features2['points'][1] = [5, 4]

reg1 = register.RegisterData(im1, features=features1)
reg2 = register.RegisterData(im2, features=features2)

feature = register.FeatureRegister(model=model.Shift,
                   sampler=sampler.Nearest)

p, warp, imWarped, error = feature.register(reg1, reg2)

print p
print im1
print im2
print imWarped

Produces:

[[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 1. 1. 1. 0. 0. 0. 0. 0. 0.] [ 0. 1. 1. 1. 0. 0. 0. 0. 0. 0.] [ 0. 1. 1. 1. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]] [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 1. 1. 1. 0. 0. 0.] [ 0. 0. 0. 0. 1. 1. 1. 0. 0. 0.] [ 0. 0. 0. 0. 1. 1. 1. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

[nfaggian]

Thanks Riaan - Neil I will look at this soon.

[NeilYager]

Yes, that seems to fix it - thanks for your help. I guess feature points should be (x, y) instead of (y, x). In my few quick tests, Bilinear seems to be fine, but I'm not sure what's going on with CubicConvolution.

[riaanvddool]

Nathan, I dont think anything should change. The indexing should probably stay [x,y] for the features.

[riaanvddool]

I think CubicConvolution is fine for images, but narrow features like a single pixel might suffer. Unless you saw something strange when using CC on normal images?

[NeilYager]

I suppose it depends on the definition of "normal images". Working with images that have a lot of high-frequency content (like sharp edges) is common in some fields, and sub-pixel accuracy is necessary. Whole pixels disappearing can be quite serious.

[riaanvddool]

Yeah, sure. But if you are looking at 16 pixels to determine 1 pixel's value, and only 1 out of the 16 is non-zero, then you are bound to get some suppression, I think. Anyway, I think that is a separate issue.

[nfaggian]

OK - Riaan thanks for tracking this down! excellent work!!