<jupyter_start><jupyter_text>Napari unwarper - auto alignment<jupyter_code>import napari
from tifffile import imread
import numpy as np

from pointpats import PointPattern
import matplotlib.pyplot as plt
original_tif = '/Users/hotte/Downloads/AVG_FOV_256_zoom11test_AfterFixPCB_00001.tif'
original_tif_ = imread(original_tif)
print(original_tif_.shape)
viewer = napari.view_image(original_tif_, rgb=False, name='Grid image')
grid_image_layer = viewer.layers['Grid image']
grid_image_layer.gamma = 1
viewer.add_points(name='Original', 
                  edge_width=1, 
                  edge_color='orangered',  
                  face_color = 'white',
                  opacity = .5,    
                  size=8,
                  blending='translucent'
              )
viewer.layers['Original'].mode ='add'
viewer.layers['Original'].symbol ='x'
layer_original = viewer.layers['Original']
usr_points = layer_original.data
#usr_points = usr_points[:,::-1]
len(usr_points)<jupyter_output><empty_output><jupyter_text>viewer.add_points(name='Firsts',                   edge_width=1,                   edge_color='yellowgreen',                    face_color = 'white',                  opacity = 1,                      size=8,                  blending='translucent'              )viewer.layers['Firsts'].mode ='add'viewer.layers['Firsts'].symbol ='o' layer_firsts = viewer.layers['Firsts']first_points = layer_firsts.data len(first_points)<jupyter_code>first_points
usr_point_pat = PointPattern(usr_points)
nn4, nnd4 = usr_point_pat.knn(4)
nnd_1 = nnd4[:,0]
med_dist = np.median(nnd_1)

print(f'Median spacing points: {med_dist:.2f}')
#plt.hist(nnd_1, bins=20);
#plt.axvline(np.median(nnd_1), color='k', ls=':')
angles = []
for no, (nn, nnd) in enumerate(zip(nn4, nnd4)): 
    filt = nn[nnd < med_dist * 1.5]
    

    for filt_idx in filt: 
        angle = np.arctan2(np.diff([usr_points[no][0], usr_points[filt_idx][0]]),
                           np.diff([usr_points[no][1], usr_points[filt_idx][1]]))[0]
        angles.append(angle)
angles = np.degrees(angles)
plt.hist(angles, bins=80);
from scipy import stats
from scipy.signal import find_peaks
x = angles
bins = np.arange(-180,180,1)
kde = stats.gaussian_kde(x , bw_method=.08)
xx = np.linspace(-250, 250, 1000)
figure = plt.figure(figsize=(10,5))
ax = figure.add_subplot(111)
ax.hist(x, density=True, bins=bins, alpha=0.5)
ax.plot(xx, kde(xx))

peaks = find_peaks(kde(xx), width=20, prominence=.001)
ax.scatter(xx[peaks[0]], np.zeros_like(peaks[0])+0.03, s=200, marker='|', color='red')
filt_peaks = [] 
for peak in np.mod(xx[peaks[0]],90):
    if np.abs(90-np.abs(peak)) < np.abs(peak):
        filt_peaks.append(90-np.abs(peak))
    else:
        filt_peaks.append(np.abs(peak))
print(f'Median orientation in degrees: {np.median(filt_peaks):.2}')
# adapted from 
# https://stackoverflow.com/questions/29630052/ordering-coordinates-from-top-left-to-bottom-right/64659569#64659569

points = []
line_locations= []
keypoints_to_search = usr_points[:]
while len(keypoints_to_search) > 0:
    a = sorted(keypoints_to_search, key=lambda p: p[1] + p[0])[0]   # find upper left point
    b = sorted(keypoints_to_search, key=lambda p: p[0] - p[1])[-1]  # find bottom left point
    line_locations.append([[a[0], a[1]], [b[0], b[1]]])

    # convert to numpy 3d point
    a = np.array([a[1], a[0]])
    b = np.array([b[1], b[0]])

    row_points = []
    remaining_points = []
    dists = []
    for k in keypoints_to_search:
        p = np.array([k[1], k[0]])
        dist = np.linalg.norm(np.cross(np.subtract(p, a), np.subtract(b, a))) / np.linalg.norm(b)   # distance between keypoint and line a->b
        dists.append(dist)
        if dist < med_dist/2:
            row_points.append(k)
        else:
            remaining_points.append(k)
    print((len(row_points)))
    points.extend(sorted(row_points, key=lambda h: h[0]))
    keypoints_to_search = remaining_points
    break
viewer.add_shapes(np.asarray(line_locations), 
                  name ='Row lines',
                  shape_type='line', 
                  face_color = [0]*4, 
                  edge_color = "mediumturquoise", 
                  edge_width = .5)

def dist(pointA, pointB):
    dist = np.sqrt(np.square(pointB[0]-pointA[0]) + np.square(pointB[1]-pointA[1]))
    return dist
line_locations
points
start_point_line = points[0]
end_point_line   = points[-1]
med_dist
dist = np.sqrt(np.square(x2 - x1) + np.square(y2 - y1) )

properties = {
    'point_idx': np.arange(len(points)),
}

viewer.add_points(name='Sorted points',
                  data=points,
                  text= {
                                'text': '{point_idx}',
                                'size': 12,
                                'weight':'bold',
                                'color': 'mediumseagreen',
                                'translation': np.array([-5, -5]),
                            },
                  properties=properties,
                  edge_width=1, 
                  edge_color='#000000',  
                  face_color = 'white',
                  opacity = .8,    
                  size=3,
                  blending='translucent'
              )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1441976254308694, 'fraction_numerical': 0.030639601685178094, 'mean_word_length': 2.5603271983640083, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 11, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '600627', 'license': 'MIT', 'n_tokens_mistral': 1968, 'n_tokens_neox': 1815, 'n_words': 309, 'path': 'Unwarp procedures.ipynb', 'repo_name': 'horsto/napari-mini-unwarp', 'size': '5222'}