Generate profiles through a 2D array at an angle without altering pixels

Question

I'd like to plot two profiles through the highest intensity point in a 2D numpy array, which is an image of a blob (i.e. a line through the semi-major axis, and another line through the semi-minor axis). The blob is rotated at an angle theta counterclockwise from the standard x-axis and is asymmetric.

It is a 600x600 array with a max intensity of 1 (at only one pixel) that is located right at the center at (300, 300). The angle rotation from the x-axis (which then gives the location of the semi-major axis when rotated by that angle) is theta = 89.54 degrees. I do not want to use scipy.ndimage.rotate because it uses spline interpolation, and I do not want to change any of my pixel values. But I suppose a nearest-neighbor interpolation method would be okay.

I tried generating lines corresponding to the major and minor axes across the image, but the result was not right at all (the peak was far less than 1), so maybe I did something wrong. The code for this is below:

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage

def profiles_at_angle(image, axis, theta):
    theta = np.deg2rad(theta)
    if axis == 'major':
        x_0, y_0 = 0, 300-300*np.tan(theta)
        x_1, y_1 = 599, 300+300*np.tan(theta)
    elif axis=='minor':
        x_0, y_0 = 300-300*np.tan(theta), 599
        x_1, y_1 = 300+300*np.tan(theta), -599
    num = 600
    x, y = np.linspace(x_0, x_1, num), np.linspace(y_0, y_1, num)
    z = ndimage.map_coordinates(image, np.vstack((x,y)))

    fig, axes = plt.subplots(nrows=2)
    axes[0].imshow(image, cmap='gray')
    axes[0].axis('image')
    axes[1].plot(z)
    plt.xlim(250,350)
    plt.show()

profiles_at_angle(image, 'major', theta)

Did I do something obviously wrong in my code above? Or how else can I accomplish this? Thank you.

Edit: Here are some example images. Sorry for the bad quality; my browser crashed every time I tried uploading them anywhere so I had to take photos of the screen.

Figure 1: This is the result of my code above, which is clearly wrong since the peak should be at 1. I'm not sure what I did wrong though.

Figure 2: I made this plot below by just taking the profiles through the standard x and y axes, ignoring any rotation (this only looks good coincidentally because the real angle of rotation is so close to 90 degrees, so I was able to just switch the labels and get this). I want my result to look something like this, but taking the correction rotation angle into account.

Edit: It could be useful to run tests on this method using data very much like my own (it's a 2D Gaussian with nearly the same parameters):

image =  np.random.random((600,600))

def generate(data_set):
    xvec = np.arange(0, np.shape(data_set)[1], 1)
    yvec = np.arange(0, np.shape(data_set)[0], 1)
    X, Y = np.meshgrid(xvec, yvec)
    return X, Y

def gaussian_func(xy, x0, y0, sigma_x, sigma_y, amp, theta, offset):
    x, y = xy

    a = (np.cos(theta))**2/(2*sigma_x**2) + (np.sin(theta))**2/(2*sigma_y**2)
    b = -np.sin(2*theta)/(4*sigma_x**2) + np.sin(2*theta)/(4*sigma_y**2)
    c = (np.sin(theta))**2/(2*sigma_x**2) + (np.cos(theta))**2/(2*sigma_y**2)

    inner = a * (x-x0)**2
    inner += 2*b*(x-x0)*(y-y0)
    inner += c * (y-y0)**2
    return (offset + amp * np.exp(-inner)).ravel()

xx, yy = generate(image)
image = gaussian_func((xx.ravel(), yy.ravel()), 300, 300, 5, 4, 1, 1.56, 0)
image = np.reshape(image, (600, 600))

I mean the 2D array you'd like to take cross-sections of: can you save that as an image and post it ? — kevinkayaks, Sep 08 '18 at 01:20
that's not adequate to reproduce your bug and fix it. it needs to be an image of exactly the 2d array you're working with — kevinkayaks, Sep 08 '18 at 01:22
@kevinkayaks if you mean I should post the actual full data, I'm not sure how to do that (it's 360,000 values so I can't post it here, and from what I understand SO doesn't have a feature to attach a file to questions). However, the graphical image of the exact data is shown above. — curious_cosmo, Sep 08 '18 at 02:03
save the data to an image and post the image. This is very different than a photo of your computer monitor. It's very unlikely anyone will help you without verifiable test data — kevinkayaks, Sep 08 '18 at 02:05
As I said in my edit, my browser on my work computer where I have the image saved would not let me upload any images (or even email them), so I was forced to take a photo of the screen and upload those from my personal laptop. I will have to try to save it with a flashdrive next time I'm at my work computer; in the meantime, I don't think that's a very good reason to downvote, as maybe someone else can provide insight on the method of slicing images at an angle (it could be with random example data). — curious_cosmo, Sep 08 '18 at 02:10
Possible duplicate of [How to extract an arbitrary line of values from a numpy array?](https://stackoverflow.com/questions/7878398/how-to-extract-an-arbitrary-line-of-values-from-a-numpy-array) — kevinkayaks, Sep 08 '18 at 02:17
Nope, I read through that and actually based my code above on the answer, but my issue is a bit different. I've now posted a very accurate example that is representative of my data -- simply generated from a 2D Gaussian. — curious_cosmo, Sep 08 '18 at 02:28
Your `y_0` and `y_1` are not centered on 300 in the minor-axis case, but you said the plot was of the major-axis... — Davis Herring, Sep 08 '18 at 02:40

kevinkayaks · Answer 1 · 2018-09-08T21:44:00.450

0

This should do it for you. You just did not properly compute your lines.

theta = 65
peak = np.argwhere(image==1)[0]
x = np.linspace(peak[0]-100,peak[0]+100,1000)
y = lambda x: (x-peak[1])*np.tan(np.deg2rad(theta))+peak[0]
y_maj = np.linspace(y(peak[1]-100),y(peak[1]+100),1000)
y = lambda x: -(x-peak[1])/np.tan(np.deg2rad(theta))+peak[0]
y_min = np.linspace(y(peak[1]-100),y(peak[1]+100),1000)
del y

z_min = scipy.ndimage.map_coordinates(image, np.vstack((x,y_min)))
z_maj = scipy.ndimage.map_coordinates(image, np.vstack((x,y_maj)))

fig, axes = plt.subplots(nrows=2)
axes[0].imshow(image)
axes[0].plot(x,y_maj)
axes[0].plot(x,y_min)
axes[0].axis('image')

axes[1].plot(z_min)
axes[1].plot(z_maj)

plt.show()

edited Sep 08 '18 at 21:44

answered Sep 08 '18 at 03:01

kevinkayaks

2,636
1
14
30

I don't think this does what I want, since `map_coordinates` uses spline interpolation, which will alter some of the pixel values shown in the profiles and make it appear smoother than it actually is. – curious_cosmo Sep 10 '18 at 17:08
this answer shows you how to define the lines like you need. The stackoverflow link I sent you shows how to interpolate (or not) however you want. I'd round the lines to integers, take the unique values, then use those to index into your image array. it just takes a tiny bit of effort and another 3 lines of code – kevinkayaks Sep 10 '18 at 18:29

Generate profiles through a 2D array at an angle without altering pixels

1 Answers1