Splat features (weights image and adding)

[Jacob Hinkle]

Commonly we do something like the following:

onesIm = Image3D(grid,mType)
splatJ = Image3D(grid,mType)
splatOnes = Image3D(grid,mType)
splatSum = Image3D(grid,mType)
splatOnesSum = Image3D(grid,mType)

SetMem(onesIm,1.0)
SetMem(splatSum,0.0)
SetMem(splatOnesSum,0.0)

for i in xrange(N):
    Splat(splatJ, phi[i], J[i])
    Splat(splatOnes, phi[i], onesIm)

    Add(splatSum,splatJ)
    Add(splatOnesSum,splatOnes)

Div(I0, splatSum, splatOnesSum)

This could be made way more efficient in both space and time.

It's done this way because ordinary splatting just takes an image and a deformation and does a single splat. But if we're doing splatting it's because we need a Jacobian-weighted sum. I propose we add the following options to Splat:

• weightIm: if present, write the weights image while we have trilerp weights computed.
• replace: just means we replace the existing images instead of adding to them. This could easily be implemented by making addition the only actual implementation of splat, and if replace is present we do a SetMem (or better yet, memset or cuda_memset) if replace is provided to zero out the original images.

With that behavior, we could replace the above chunk with the following:

splatSum = Image3D(grid,mType)
splatOnesSum = Image3D(grid,mType)

SetMem(splatSum,0.0)
SetMem(splatOnesSum,0.0)

for i in xrange(N):
    Splat(splatSum, phi[i], J[i], weightIm=splatOnesSum, replace=False)

Div(I0, splatSum, splatOnesSum)

I propose keeping the current interface to Splat but adding the above as optional arguments that get routed to templated versions of Splat that can do combinations of the above stuff. I think this just amounts to needing four new template instantiations in some of the *FieldOpers.{h,cxx} files.

[Jacob Hinkle]

Along with this, we should include a unit test that checks that this is actually the adjoint of ApplyH. I will probably have time to do all this after Feb 1

[Jacob Hinkle]

By the way, I think we could also use this to handle the normalize=True case easily, even though I'm not aware of a single time we've used this.

[Jacob Hinkle]

Initial attempt at splatting with weights (CPU only)

See #15. This needs testing. Splat() looks much more complicated on the GPU and I'm not sure exactly all that's going on there, so until I do I don't wanna touch it. The CPU approach is straight-forward so I've updated that part. If you have a MEM_HOST image and corresponding weights image allocated you can pass those in as the first two parameters to Splat and not send normalize option. Since this might touch quite a bit of existing work, and since there is a workaround, we should wait to merge this until we have a unit test verifying all the possible scenarios.

→ <<cset fa432ec982b3>>

[Sam Preston]

Yeah, splatting is done on the GPU with fixed-point precision. Not that much different but you have to be careful what operations you use on the data. I need to look at this again too to make sure I understand it.

[Jacob Hinkle]

For reference, see my comments in Shooting.h

[Jacob Hinkle]

In this commit message, I noted that Splat is also not currently the adjoint of ApplyH. They should be for every choice of boundary condition and interp method. Also, the default boundary condition and interp method should be the same for both functions. I'm not sure how to augment that test so that it checks all these conditions, but I will look into it when I get time to fix this problem.

[Jacob Hinkle]

See #32. Whenever we overhaul splat, we should make the interface match the interpolation interface precisely, in terms of weights computation and background strategies.