doc update for hansenlaw, vmi

[stggh]

Updated docs for doc/transform_methods/hansenlaw.rst and abel.tools.vmi.py.

NB I renamed calculate_angular_distributions() to radial_integration() to better align with the partner method angular_integration(). Modified examples to suit.

Hmm, some of the listed changes were just a 80 column reformat, sorry.

@DanHickstein livehtml is neat, and certainly speeds up the editing process.

Edit: NB changed name of parameter method to center in abel.tools.center.find_center(image, center="<method>") to match center_image(image, center="<method>")

[stggh]

Issue: abel.tools.symmetry.put_image_quadrants() does not know about the original image size.

get_image_quadrants() correctly handles any shaped image, adding one row or column if the respective axis size is odd. The reassembly, put_image_quadrants() depends on the odd-evenness of both axes. The single passed parameter, odd_size=True assumes, incorrectly, a square original image.

Ideas:

1. tuple odd_axis=None, 0, 1, or (0,1)
2. We encapsulate get/put_image_quadrants() into a class, keeping the knowledge of the original axis sizes

Edit:
3. pass the original image shape to put_image_quadrants.

I (now) prefer 3.

[stggh]

Point 3. has been implemented.

[Steve:.+Abel/abel/tests]$ python test_tools_symmetry.py 
test image of shape  (4, 5)
 odd-size axes tagged with '-1'
[[ 1.  1. -1.  0.  0.]
 [ 1.  1. -1.  0.  0.]
 [ 2.  2. -1.  3.  3.]
 [ 2.  2. -1.  3.  3.]]

reoriented quadrant Q0
[[-1.  0.  0.]
 [-1.  0.  0.]]

reoriented quadrant Q1
[[-1.  1.  1.]
 [-1.  1.  1.]]

reoriented quadrant Q2
[[-1.  2.  2.]
 [-1.  2.  2.]]

reoriented quadrant Q3
[[-1.  3.  3.]
 [-1.  3.  3.]]

reassembled image, shape =  (4, 5)
[[ 1.  1. -1.  0.  0.]
 [ 1.  1. -1.  0.  0.]
 [ 2.  2. -1.  3.  3.]
 [ 2.  2. -1.  3.  3.]]

another test:

[Steve:.+Abel/abel/tests]$ python test_tools_symmetry.py 
test image of shape  (5, 5)
 odd-size axes tagged with '-1'
[[ 1.  1. -1.  0.  0.]
 [ 1.  1. -1.  0.  0.]
 [-1. -1. -1. -1. -1.]
 [ 2.  2. -1.  3.  3.]
 [ 2.  2. -1.  3.  3.]]

reoriented quadrant Q0
[[-1.  0.  0.]
 [-1.  0.  0.]
 [-1. -1. -1.]]

reoriented quadrant Q1
[[-1.  1.  1.]
 [-1.  1.  1.]
 [-1. -1. -1.]]

reoriented quadrant Q2
[[-1.  2.  2.]
 [-1.  2.  2.]
 [-1. -1. -1.]]

reoriented quadrant Q3
[[-1.  3.  3.]
 [-1.  3.  3.]
 [-1. -1. -1.]]

reassembled image, shape =  (5, 5)
[[ 1.  1. -1.  0.  0.]
 [ 1.  1. -1.  0.  0.]
 [-1. -1. -1. -1. -1.]
 [ 2.  2. -1.  3.  3.]
 [ 2.  2. -1.  3.  3.]]

[stggh]

@DhrubajyotiDas and @DanHickstein
Given the timing issues with the basis set methods (basex and three_point) it may be a good idea to standard all the basis set examples to one size image. The (random) examples that I have run produced:

basex_basis_151_301_151_151.npy
basex_basis_1_501_1_251.npy
basex_basis_251_501_251_251.npy
basex_basis_256_511_256_256.npy
basex_basis_260_693_260_347.npy
basex_basis_512_1023_512_512.npy

eg. example_hansenlaw_basex.py @DhrubajyotiDas zoom factor may be zoom(im, 0.4892578125) to produce a 501x image size, which then matches the sample image sizes. I have made this change to this PR.

[DhrubajyotiDas]

@stggh that's a great point. I'll look into standardizing the basis set examples.

[stggh]

This PR is ready to be merged. The main changes:

• Updated the text in doc/transform_methods/hansenlaw.rst
• abel.tools.center.center_image() variable name change of method= to center=, to match the other methods.
• abel.tools.symmetry.put_image_quadrants() now reassembles non-square images, original_image_shape variable replaced the single odd_size that applied to square images
• NB abel.tools.symmetry.get_image_quadrants() I removed averaging. The averaging denominator complicates the use of use_quadrants=(boolean tuple), to such an extent that it is virtually useless. e.g. it is impossible to map a single quadrant, Qi, to the whole image. I have reached the conclusion that it is better not to average. Removed averaging/average from docs.
• updated the unit test test_tools_symmetry and examples to reflect the above changes

[DanHickstein]

Ah, yes, I see that you were saying about \infty before - sorry for the confusion.

[DanHickstein]

Thanks for making these changes @stggh!

That is an interesting issue with put_image_quadrants() and I see no easy solution. Passing the original image size as you have done, works, and is robust, but I think it might be annoying to users. Would it be a good idea for put_image_quadrants() to have a default of original_image_shape='odd', whereby both the rows and the columns would be assumed to have odd values? (Not that the image is square, just that both rows and column are odd.) Or, your current implementation is fine.

I think that the averaging the quadrants is the way to go, as it avoids weird situations where the adding quadrants would lead to situations where half of the image is twice as intense as the other half. The existing code seems like it would throws a divide-by-zero errors if the user attempted to map a single quadrant to the whole image (i.e., symmetry_axis=(0,1)). We could likely fix this by first including:

is symmetry_axis==(0, 1):
    Q = (Q0*use_quadrants[0]+
         Q1*use_quadrants[2]+
         Q2*use_quadrants[2]+
         Q3*use_quadrants[3])/(np.sum(use_quadrants))

    return Q, Q, Q, Q 

Whatever solution we adopt, we should probably also check for situations where the user selects options for symmetry_axis and use_quadrants that don't make sense. For example, anything that does not produce a full image should be caught an an error raised:

symmetry_axis == None, any quadrant == False    -> at least one quadrant is empty
symmetry_axis == 0, top quadrants both False    -> top of image is empty
symmetry_axis == 0, bottom quadrants both False -> bottom of image is empty
symmetry_axis == 1, left quadrants both False   -> left of image is empty
symmetry_axis == 1, right quadrants both False  -> right of image is empty

With these changes, I don't think that averaging should pose a problem, and would produce something closer to what the user expects.

And I agree about the basis sets. Paring down to just one (or two) basis sets should make the examples a lot faster to run.

[stggh]

@DanHickstein, sorry for the delay I was trapped in a gridlocked road, due to a car accident further along.

I gave up with averaging because it causes issues about which quadrants should be averaged first. eg. with the old pre-PR code try selecting only Q0, ..., Q3 and see the issue. No doubt @DanHickstein you will come up with a clever solution.

I will revisit. Just heading off to the river for a few hours to escape from the 33C autumn heat.

[DanHickstein]

it causes issues about which quadrants should be averaged first

Yes, you are correct that this is a glaring flaw with the current code.

But, the order of averaging should only be an issue with symmetry_axis=(0,1), in the other cases there is only one (or zero) averaging steps for each pair of quadrants.

If we simply check for symmetry_axis=(0,1) first, then we should avoid this problem.

[DanHickstein]

Okay, I think that the following code (included before the if 0 in symmetry_axis line in get_image_quadrants) should allow for the averaging to work:

if ((symmetry_axis == [None] and (use_quadrants[0]==False
                                   or use_quadrants[1]==False 
                                   or use_quadrants[2]==False
                                   or use_quadrants[3]==False)) or        # at least one empty
        (symmetry_axis == [0] and use_quadrants[0]==False and use_quadrants[1]==False) or # top empty
        (symmetry_axis == [0] and use_quadrants[2]==False and use_quadrants[3]==False) or # bot empty
        (symmetry_axis == [1] and use_quadrants[1]==False and use_quadrants[2]==False) or # left empty
        (symmetry_axis == [1] and use_quadrants[0]==False and use_quadrants[3]==False)): # right empty
        raise ValueError('At least one quadrant would be empty. Please check symmetry_axis '
                          'and use_quadrant values to ensure that all quadrants will have a '
                          'defined value.')

    if symmetry_axis==(0, 1):
        Q = (Q0*use_quadrants[0]+
             Q1*use_quadrants[2]+
             Q2*use_quadrants[2]+
             Q3*use_quadrants[3])/(np.sum(use_quadrants))

        return Q, Q, Q, Q

And here is code to test all situations, it should probably be incorporated into a unit test:

import abel
import numpy as np

z = np.zeros((5, 5))

for symmetry_axis in (None, 0,1,(0,1)):
    for Q0 in (True, False):
        for Q1 in (True, False):
            for Q2 in (True, False):
                for Q3 in (True, False):

                    print symmetry_axis, Q0, Q1, Q2, Q3

                    try:
                        abel.tools.symmetry.get_image_quadrants(z,  
                                    symmetry_axis=symmetry_axis, 
                                    use_quadrants=(Q0,Q1,Q2,Q3))
                    except ValueError as e:
                        print e

[stggh]

@DanHickstein you have a very logical brain. The top section of your code is quite a mouthful. I have implemented it + one additional test case symmetry_axis=(0,1) and (Edit later corrected) not np.any((use_quadrants)).

I have returned (what I could find) the doc references to average.

Test code run, not yet implemented as a unit test.

Edit: I am not 100% get_image_quadrants() is working correctly. Whole and single quadrant speed profiles look too similar, but may be that's the case for dribinski sample image. It's late. I will take another look tomorrow.

[DanHickstein]

Each row corresponds to a definite theta, that the reconstruction does not care about the other rows in an image was hard to get my head around.

Yeah, with BASEX three-point, this important consequence of assuming only cylindrical symmetry (independent rows) gets hidden in the matrix multiplications, while the other algorithms (HL and Direct) make it clear that each row is independent.

I think that the reason that this seemed unintuitive to you (and me) is that photolelectron/photoion images have a higher degree of symmetry than just cylindrical (not quite sure what word to use to describe it. Basically they have intensity that is "slowly varying" as a function of theta). To get good transforms with low signal-to-noise photoelectron data, it's probably a good idea to use a method that makes more assumptions than just cylindrical symmetry, like pBASEX or POP (#30).

---

As far as the crazy if statement in the above. It could be split into several if statements each with their own ValueError. This would take up more lines of code, but would likely be more informative for the user.

[DanHickstein]

And good point about symmetry_axis=(0,1) and not np.any((use_quadrants)).

Actually, shouldn't we just generally raise a ValueError if not np.any((use_quadrants))? (i.e., all symmetry cases require at least one quadrant.)

[DhrubajyotiDas]

@DanHickstein Just to add a clarification to the independent-row discussion, with three-point, the independence is pretty explicit, as seen in this line.

    for i, P in enumerate(IM):
        inv_IM[i] = np.dot(D, P)
    return inv_IM

We are, quite literally, cycling through every row in IM one-at-a-time and doing stuff to it, and rebuilding the transformed inv_IM row-by-row.

---

I think that the reason that this seemed unintuitive to you (and me) is that photolelectron/photoion images have a higher degree of symmetry than just cylindrical (not quite sure what word to use to describe it. Basically they have intensity that is "slowly varying" as a function of theta).

This example might help visualize this. The left is a line-of-sight image of a flame, the the right is the BASEX inversion of that image. You can think of this image as composed of only quadrants 1 and 0.

[stggh]

@DhrubajyotiDas neat, and nice colours.

get/put_image_quadrants() is almost doing the correct thing. If I select one quadrant, there is a blank cross in the image center

Edit: This was just the input image. I think all may be well with get/put

Edit2: Not quite, there remains a vertical blank space ...

/Edit2

Edit corrected code added % 2 to size calculation

Apologies for all the re-edits.

# -*- coding: utf-8 -*-
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

import abel

def gauss(r, r0, sigma):
    return np.exp(-(r-r0)**2/sigma**2)

def test_quad(image_shape=(201, 201)):
    rows, cols = image_shape
    r2 = rows//2 + rows % 2
    c2 = cols//2 + cols % 2
    x = np.linspace(-c2, c2, cols)*201/c2
    y = np.linspace(-r2, r2, rows)*201/r2

    X, Y = np.meshgrid(x, y)
    R, THETA = abel.tools.polar.cart2polar(X, Y)

    IM = np.zeros(image_shape)

    R0, dummy = abel.tools.polar.cart2polar(X[:r2, -c2:], Y[:r2, -c2:])
    IM[:r2, -c2:] = gauss(R0, 50, 4)
    R1, dummy = abel.tools.polar.cart2polar(X[:r2, :c2], Y[:r2, :c2])
    IM[:r2, :c2] =  gauss(R1, 100, 4)
    R2, dummy = abel.tools.polar.cart2polar(X[-r2:, :c2], Y[-r2:, :c2])
    IM[-r2:, :c2] = gauss(R2, 150, 4)
    R3, dummy = abel.tools.polar.cart2polar(X[-r2:, -c2:], Y[-r2:, -c2:])
    IM[-r2:, -c2:] = gauss(R3, 200, 4)

    print("Enter use_quadrants, symmetry_axis eg (1,0,0,0),(0,1)")
    plt.ion()
    while True:
        uq, sa = input("use_quadrants, symmetry_axis? ")

        try:
            Q = abel.tools.symmetry.get_image_quadrants(IM, symmetry_axis=sa,
                         use_quadrants=uq)
        except ValueError, e:
            print(str(e))
            continue

        rIM = abel.tools.symmetry.put_image_quadrants(Q,
                         original_image_shape=IM.shape, symmetry_axis=sa)

        lbl = "{:s} {:s}".format(str(uq), str(sa))
        flbl = lbl.translate(None, "(),' '")
        plt.subplot(211)
        plt.title("input: "+lbl)
        plt.imshow(IM)
        plt.subplot(212)
        plt.title("reassembled")
        plt.imshow(rIM)
        plt.savefig(flbl+".png",dpi=100)
        plt.draw()

if __name__ == "__main__":
    test_quad()

[DanHickstein]

@DhrubajyotiDas, that is beautiful! Is that real data?

@stggh, that looks good with the put/get methods. Should this be merged now?

[stggh]

Not quite ready, sorry for all the re-edits. There remains a vertical center blank line, which is no doubt a bug in my test code, above, but it will need to be found.

Edit (I have not learnt): Part of the issue is the way I construct the test image. The whole odd-size image results in individual quadrants are even-size, and some have some common boundaries.

I now think that the code is ready to be merged.

pushed the raise a ValueError if not np.any((use_quadrants)) correction mentioned above

[DhrubajyotiDas]

@DanHickstein Yes, pretty much straight out of the instrument without any processing.

I was considering putting up the image and its transform as a flame example on this repo, but they are fairly bulky (> 5 MB). I know in the past we've taken steps to keep the repo down to a reasonable size, which is one of the reasons why I haven't done so yet.

@stggh get/put_image_quadrants() looks pretty cool. I will definitely find a use for it in some of my other projects.

[stggh]

@DhrubajyotiDas get/put_quadrants() is my nemesis. It opens up a can of worms concerning how an image may be divided, further compounded by the odd-width requirements of the Abel transform methods.

Note, the more clever Python code is due to @rth and @DanHickstein. I just implemented the concept of quadrant combining (from my own VMI analysis), and it has evolved in PyAbel.

[stggh]

@DanHickstein this PR merge has done something weird to the doc links, e.g. README.rst

:doc:Hansen–Law <transform_methods/hansenlaw>

[DanHickstein]

I think that this was my previous PR. I made the links look good when the documentation is built (see http://pyabel.readthedocs.org/en/latest/readme_link.html) but they don't work on GitHub.

I don't know of an easy fix to make them work in both situations... Maybe I will just shorten the readme and move some of the information on the transform methods to a different doc?

[DanHickstein]

A solution to this problem with the links is discussed here: http://comments.gmane.org/gmane.comp.python.sphinx.devel/7539 but it looks a little complicated to implement.