natelust/least_asymmetry

"RuntimeError: Arrays must be square"

Closed this issue · 3 comments

I am attempting to cycle through a stack of ~60k images, each of size 32x32 and find the center of the Spitzer PSF that is expected to be at the dead center of the image (pixel 15x15).

I use this code:

def fit_least_asymmetry_centering(imageCube, yguess, xguess, asym_rad):
    '''
        imageCube : nD-array (assume: (60k, 32, 32))
        yguess    : initial y center guess (assume: 15)
        xguess    : initial x center guess (assume: 15)
        asym_rad  : width of asymmetry subframe -- attempt to control corner cases (assume: 5)
    '''
    y,x = 0,1
    
    nFrames = imageCube.shape[0]
    
    yinds0, xinds0 = indices(imageCube[0].shape)
    
    nAsymParams = 2 # Xc, Yc
    centering_LeastAsym  = np.zeros((nFrames, nAsymParams))

    for kframe in tqdm(range(nFrames), total=nFrames):
        center_asym = least_asymmetry.actr(imageCube[kframe], [yguess, xguess], \
                           asym_rad=asym_rad, asym_size=5, maxcounts=2, method='gaus', \
                           half_pix=False, resize=False, weights=False)[0]

        try:
            centering_LeastAsym[kframe]  = center_asym[::-1]
        except:
            print('Least Asymmetry FAILED: Setting centering_LeastAsym[%s] to Initial Guess: [%s,%s]' % (kframe, yguess, xguess))
            centering_LeastAsym[kframe]  = np.arange([yguess, xguess])

    return centering_LeastAsym

from sklearn.externals import joblib
from pylab import *;ion()
import least_asymmetry

xgf_imageCube = joblib.load('xgf_image_cube_array.pickle') # input 60k x 32 x 32 nD-array

yguess   = 15
xguess   = 15
asym_rad = 5 # tried 8 before too

fit_least_asymmetry_centering(xgf_imageCube, yguess, xguess, asym_rad)

and receive this error RuntimeError: Arrays must be square (see Traceback below)

This error is generated from the make_asym.cc code on line 35.

I can process about 1% of the images before this error is triggered. So I assume that the first 600 images worked fine, but then they 601st image breaks.

I added std::cout << data.shape[0] << '\t' << data.shape[1]; to make_asym.cc on line 26 just to see what shape the python generator was sending to cpp code. The error is generated when the std::cout reads 17 0, which means that data.shape[1] == 0

I tracked the error this far, but cannot figure out why it's happening.

Any idea?

Right now, I bypass the whole problem with a simple "try/except" in the python script.

Traceback (most recent call last):
  File "least_asym_test.py", line 68, in <module>
    fit_least_asymmetry_centering(xgf_imageCube, 15,15,8)
  File "least_asym_test.py", line 51, in fit_least_asymmetry_centering
    half_pix=False, resize=False, weights=False)[0]
  File "/Users/jonathan/anaconda/lib/python3.6/site-packages/least_asymmetry-0.2-py3.6-macosx-10.7-x86_64.egg/least_asymmetry/asym.py", line 248, in actr
    asym = np.array(list(map(make_asym, views, lb_views, w_truth_dup)))
RuntimeError: Arrays must be square

By using the try/except bypass, I can say that this error occurs about 25 times in 60k images.

This issue occurs when there is no minimum in the asymmetry space. Basically the routine is following the gradient in asymmetry space trying to get the minimum value to be at the center of a square array. If in the course of following the gradient, the array which is used to calculate an asymmetry value extends beyond the data array it is attempting to measure this error will be thrown. In other words, it is trying to calculate a radial profile around a certain x,y position but there is no data in the input array from which to generate the profile. The code used a pre-calculated square array do do this calculation (it is important it always be symmetric so a rectangle would not work) for the sake of efficiency. The routine is checking the pre-computed array matches the dimensions of a view into the data array and throws an error if they do not match.

There are a few quick options you can go with.

  1. give more area (a larger array) to the routine and hope that this will be enough such that the gradient in asym space flows down hill to a minimum.

  2. decrease the asym_rad to something like 2 or 3 pixels. If the gradient in asym space is such that you would find a minimim near an edge this may give you enough pixels to work with. Remember you need enough data to correspond to the asym transformation radius plus the asym image radius. in other words if you want to transform a 5x5 area of an image using a 5 pixel asym radius you will need a data array of at least 10x10 (assuming the center is already in the center of the data image) if is not, you will need additional pixels in the dimensions the routine needs to walk to find the center. If it must go up one pixel and right one pixel, you will need a data array of at least 11x11. You may be able to just decrease the asym_rad for the images that fail the try block (decreasing the rad will also decrease the accuracy.

  3. Investigate the failing images. It may be possible that something weird is going on (bluring large noise spikes, etc.) and the images are bad anyway. At least this may give you an idea of what the algorithm is seeing.

Does this help? I think I probably can put in a more helpful error message to give a hint to this for future uses.

Yes. That definitely helped. All 27 images had a cosmic ray hit (strong single-pixel spike) nearby the PSF. So the algorithm was definitely being led astray.

I think that it would be fruitful for the algorithm to print an error message that flags the frames and says "Frame ##### FAILED; examine more closely; assigning pathological value for centering". Then I would have looked at those frames and said "yep. these are cosmic ray hits". Instead, I spent about 2 hours digging through the python and cpp code looking for the size of the arrays.

The error message "Array must be square" implied to me that the input array (my images) were not square, which made no sense to me.

Thank you for the help!