Keeping track of image geometry

[kmsquire]

I've been implementing the kernel hough transform (http://www2.ic.uff.br/~laffernandes/projects/kht/), and have been finding it challenging to keep the coordinates straight in my code.

Using im["x", 100, "y", 200] kind of works, but I'm not thrilled with it, and would prefer something more flexible.

What I'd like to do is write a thin Axes wrapper that takes into account the Image coordinate information and allows accessing the image in a consistent and useful manner for the task at hand--e.g., if I want to access an image like a matrix, I could create a set of axes that index by m[row, col] or m[row, col, color].

But if I wanted Cartesian axes instead, I could put an origin at, e.g., the lower left, and have my image occupy the upper-right quadrant. Or I could put the origin at the image center, flip, or scale it, as needed. It should be easily extendable to image sequences as well.

Questions:

1. Is there functionality like this already that I missed?
2. What is/should be the overlap with other packages? It would seem that Grid.jl might be useful here (and may already provide this functionality?). Plotting systems already use axes, of course, but don't interact well with Images.jl. PyPlot has this functionality, but also doesn't interact well with Images.jl

Thoughts? If this doesn't exist already, would this be a useful addition to Images.jl? Or should it be a separate package?

[timholy]

Let me get this straight: what you're hoping for is a high-performance method of writing your code in a generic way that allows you to access individual pixels?

That's pretty much the reason for the existence of Cartesian. If you can include a short snippet of what you're trying to do, I can coach you through writing it.

[kmsquire]

Thanks, Tim. I've only seen Cartesian used for creating arbitrarily dimensional loops, and I wasn't aware of it's other capabilities.

I'll give a specific example when I'm in front of a computer, but for the KHT, the struggles I've had have been more along the lines of

1. Keeping matrix (row, column) vs Cartesian (x,y) indexing straight, since I'm working with both Images and images-as-matrices.
2. Keeping my axes oriented. Lines in the hough transform are specified using (rho, theta) notation, and I find that easier to reason about in Cartesian coordinates than with the Y-axis flipped.

[kmsquire]

After looking a little closer at Cartesian, that's not really what I'm looking for.

What I'm really after is a more flexible, less costly (?) implementation for im["x", 100, "y", 200]--i.e., I want to read and write pixels and add annotations at arbitrary locations using my own (consistent) coordinate system.

[kmsquire]

"consistent" = "consistent no matter the intrinsic layout of the data"

[timholy]

I thought you were worried about storage order in the sense of color or time dimensions, and how you'd traverse the neighbors without, e.g., worrying about whether color was first or color was last.

Yes, keeping track of geometry is always a pain, especially when your images don't allow you to encode what that geometry is (much of the raison d'êtra for Images). If this is just a case of computing theta as atan2(y,x) or atan2(x,y), is that really something you need to be worried about at the pixel level? Why not figure out the geometry before or after the loop and massage the outputs accordingly?

[timholy]

That's not to say that your Axes wrapper wouldn't be interesting---I'm not fully sure I understand your intention, but it sounds like it's an interesting experiment in relabeling how arrays are indexed. Could be fun & useful.

One caution, though: for cache-efficiency reasons, IMO it's best for core algorithms to traverse arrays in their native storage order. (For many algorithms cache-order tends to be a big effect, e.g., 5-fold, so we're not talking peanuts here.) So anything that "lies" about that makes it harder to write fast algorithms. Rather than wrapping indexing, how about transforming coordinates at the time you interpret them? For your case here, that's a little like what I suggested above. In particular, flipping the x and y axes is equivalent to the transformation theta -> pi/2-theta. If you don't want to mess with the array after building it, you could define theta = a + c*atan2(y,x) and either set a=0, c = 1 or a = pi/2, c = -1 depending on the geometry of the image. A more general version of this might, for example, let you look up xreal, yreal = tform(obj, xindex, yindex). Would that do what you want?

[kmsquire]

I thought you were worried about storage order in the sense of color or time dimensions, and how you'd traverse the neighbors without, e.g., worrying about whether color was first or color was last.

Possibly that as well, although for the KHT work, I'm working with BW images. (There is some nice work on edge detection in color images that I'd like to look at eventually, though.)

Why not figure out the geometry before or after the loop and massage the outputs accordingly?

I've tried that somewhat, and haven't gotten things working yet, but that might just be because I haven't tried hard enough yet. ;-)

Thanks for the suggestion, re:traversing in native storage order (which leads back to where you started ;-). I was generally aware of that concern, but not intimately. Still, it would be nice to traverse the array in storage order, but still index in a prescribed manner.

I'll work at it a bit, see what I come up with, and post back here with something more specific. Cheers!

[timholy]

Still, it would be nice to traverse the array in storage order, but still index in a prescribed manner.

Your indexing describes how it's traversed; any time you access a particular location via indexing, that defines the traversal.

If you want cache efficiency, the closest you can come, I think, is "traverse in storage order but then tell me where I am in the terms I care about." That's the idea behind my xreal, yreal = tform(obj, xindex, yindex) suggestion. obj could be an instance of stub-type SwapMe, for example, if you effectively want to traverse in slow-axis-first order.

[timholy]

Noticing this comment again:

Plotting systems already use axes, of course, but don't interact well with Images.jl. PyPlot has this functionality, but also doesn't interact well with Images.jl

Can you elaborate on the problems? Is it that other packages define an Image type?

[kmsquire]

Can you elaborate on the problems? Is it that other packages define an Image type?

Sorry, that wasn't clear at all. It's really just that none of those packages a) take advantage of the Image type, or b) just generally work well with images.

PyPlot comes closest, with PyPlot.imshow(permutedims(img.data, (3,2,1)). (I actually just discovered this... I guess #96 has some information on how to avoid permutedims.) (Edit: Here: #96 (comment))

Winston doesn't seem only to display single channel pseudocolor images, and Gadfly is great for plotting, but doesn't work with images at all, as far as I can tell.

All of which is fine, really, since I was looking for a way to interact with the original data using an arbitrary axis, and was less concerned about plotting.

[timholy]

OK. Presumably the best approach would be for someone who needs that functionality to submit pull requests to those other packages. Nice that PyPlot already has done some of this for us.

Package interactions remain a bit of a challenge, unfortunately.

[timholy]

Update for @kmsquire, @jiahao (re JuliaImages/ImageView.jl#34). I have new implementations of ArrayViews brewing, ones that are general enough to handle any AbstractArray (I've already discussed this with Dahua, since he wrote ArrayViews). I have also done some initial work on a package called NamedAxesArrays to be built on top of that. Using the magic of stagedfunctions I think it will be possible to build all of this so that A[...] returns another NamedAxesArray for anything more than a single pixel, so you won't lose track of orientation, and yet it will be just as fast as returning an Array is now.

Moreover, I'm hoping each type of NamedAxesArray will have a custom iterator, so to calculate the intensity-weighted mean you can say

for iter in eachelement(A)
    intensity = A[iter]
    total += intensity
    xsum += iter.x*itensity
...
end
xmean = xsum/total

The key idea is to traverse the array in cache-friendly order but make sure the user can easily extract the current coordinate in physically-meaningful terms. Letting the iteration be handled in the background but giving the user a handle on the value of the :x (or :anterior, etc) coordinate at any point seems like it should provide much of what we want?

Now take a View of a NamedAxesArray, and I think you have all the functionality (shifting the origin, flipping upside down with 512:-1:1, etc) that you outlined above.

I'd be happy to have help in these endeavors! Would like to have more time to hack on these but so many other things (like HDF5 😜) are distracting me!

[kmsquire]

Sounds great, Tim!

[jiahao]

This does indeed open up a wide array of possibilities. :)

[jiahao]

Seriously, though, this could be really useful for geospatial mapping.

[timholy]

Closed by #577