Add stereographic projection

[DSimonne]

See https://pubs.acs.org/doi/abs/10.1021/acsnano.3c11534

[ccechatelier]

Hi, I will clean my code, and suggest a version that you can work with asap!

[ccechatelier]

In order to provide this projection, one need to use the interpolation parameters provided by cdiutils. Here is an example for the initial particle in the aforementionned paper (scan 163) :

import numpy as np
import matplotlib.pyplot as plt
import colorcet
from cdiutils.converter import Interpolator3D # maybe this has changed in cdiutils, need to check
from scipy.interpolate import griddata

scan = 163
filename_orthoparam = f'S{scan}/cdiutils_S{scan}_interpolation_parameters.npz'
filename_pynxdata = f'S{scan}/pynx_phasing/S{scan}_pynx_input_data.npz'

M = np.load(filename_orthoparam,'r')['q_space_linear_transformation_matrix']
Qspacetransitions = np.load(filename_orthoparam,'r')['q_space_transitions']
Qspaceshift = np.load(filename_orthoparam,'r')['q_space_shift']

I = np.load(filename_pynxdata, 'r')['data']

q_lab_interpolator = Interpolator3D(original_shape=I.shape,
                                    original_to_target_matrix=M,
                                    target_voxel_size=[3e-4, 3e-4, 3e-4],
                                    verbose=True)

I_ortho = q_lab_interpolator(I)

X, Y, Z = Interpolator3D.zero_centered_meshgrid(I_ortho.shape, [3e-4, 3e-4, 3e-4])
X += Qspaceshift[0]
Y += Qspaceshift[1]
Z += Qspaceshift[2]
Xmean = Qspaceshift[0]*0.9999
Ymean = Qspaceshift[1]*0.9999
Zmean = Qspaceshift[2]*0.9999

nxI, nyI, nzI = I_ortho.shape

#take only the upper part of the sphere
I_partial = np.copy(I_ortho)
I_partial[:, :, :nzI//2+1] = 0

#define matrix of radii
volume_R = np.sqrt((X - Xmean)**2 + (Y - Ymean)**2 + (Z - Zmean)**2)

#### THOSE VARIABLES ARE CRITICAL
radius_mean = 0.015 # where to look at, distance from the Bragg (center of the box)
dr = 0.0002 # thickness of the volume to look at

# Following mask : True where we keep the data
mask = np.logical_and( (radius_mean - dr) < volume_R, volume_R < (radius_mean + dr) )

# make up some randomly distributed data
x = radius_mean*-(X - Xmean)/(radius_mean - (Z - Zmean))
y = radius_mean*-(Y - Ymean)/(radius_mean - (Z - Zmean))

x1 = x[mask] / radius_mean * 90
y1 = y[mask] / radius_mean * 90

# Some plottwists
I_partial = I_partial[:, ::-1, ::-1]
mask = mask[:, ::-1, ::-1]

I_interest = I_partial[mask]

grid_x, grid_y = np.mgrid[-90:90:250j, -90:90:250j]
I_polar = griddata( (x1, y1), I_interest, (grid_x, grid_y), method='linear')
I_polar[np.isnan(I_polar)] = 0

# to clean the outskirt of the plot
I_log_threshold = -0.5
I_polar_log10[np.where(I_polar_log10 < I_log_threshold )] = np.nan

######### PLOT

fig, ax = plt.subsplots()

plt.contourf(grid_x, grid_y, I_polar_log10, 20, cmap='turbo', vmin=-0.5, vmax=1.1)

plt.axis('equal')
plt.grid()
fig.patch.set_visible(False)
ax.axis('off')

# the different angle circles
circle=plt.Circle((0,0),10,color='0.1',fill=False,linestyle='dotted', linewidth=2)
fig.gca().add_artist(circle)
circle=plt.Circle((0,0),30,color='0.1',fill=False,linestyle='dotted', linewidth=2)
fig.gca().add_artist(circle)
circle=plt.Circle((0,0),50,color='0.1',fill=False,linestyle='dotted', linewidth=2)
fig.gca().add_artist(circle)
circle=plt.Circle((0,0),70,color='0.1',fill=False,linestyle='dotted', linewidth=2)
fig.gca().add_artist(circle)
circle=plt.Circle((0,0),90,color='k',fill=False, linewidth=2)
fig.gca().add_artist(circle)

for ii in range(0,365,45):
    ax.plot([0, 90*np.cos(ii*np.pi/180)], 
                [0, 90*np.sin(ii*np.pi/180)], 
                color='0.1',
                linestyle='dotted', 
                linewidth=2.5)

#if one wants to save the polar data
file_data= open(f'pole_stereo_Iortho.dat', "w")  
string=''
for ii in range(len(grid_x)):
    for jj in range(len(grid_y)):
        string+=str(grid_x[ii,0])+" "+str(grid_y[0,jj])+" "+str(I_polar[ii,jj])+'\n'
file_data.write(string)
file_data.close()

Initial Bragg peak :

Half of it :

Selected region :

Result :

Hope it's clear enough! Cheers, C.