orientations/nm_ors.py at main · BonnenLab/orientations · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
import json
import numpy as np
import pandas as pd
from math import pi
import os
import sys
import cv2
from scipy.signal import fftconvolve
import pyrtools as pt

# Assumes folder is a folder of images and nothing else
folder = "path/to/images/"
files = sorted(os.listdir(folder))
save_folder = "path/to/save/" # the path to save the histograms to

n_ors = []

gp = np.array([0.037659, 0.249153, 0.426375, 0.249153, 0.037659])
gp = gp.reshape((1,5))
gd = np.array([-0.109604, -0.276691, 0.000000, 0.276691, 0.109604])
gd = gd.reshape((1,5))
bfilt = np.array([1, 2, 1])/4
bfilt = bfilt.reshape(3, 1)

for i in range(len(files)):
    try:
        im = cv2.imread(os.path.join(folder, files[i]))
        im = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
        H, W = im.shape[:2]

        # side length of the largest square that fits
        side = min(H, W)

        # top-left corner of centered crop
        top = (H - side) // 2
        left = (W - side) // 2

        # crop
        im = im[top:top + side, left:left + side]

        Gpyr = pt.pyramids.GaussianPyramid(im,6)

        # all hists is a list of the 6 histograms for the current frame
        all_hists = []
        # curr or n is a list of the current total orientations for each of the 6 histograms
        curr_or_n = []
        for j in range(6):

            Ix = pt.pyramids.c.wrapper.corrDn(pt.pyramids.c.wrapper.corrDn(Gpyr.pyr_coeffs[j, 0], gd), gp.T)
            Iy = pt.pyramids.c.wrapper.corrDn(pt.pyramids.c.wrapper.corrDn(Gpyr.pyr_coeffs[j, 0], gp), gd.T)

            Mxx = fftconvolve(fftconvolve(Ix*Ix, bfilt, 'valid'), bfilt.T, 'valid')
            Myy = fftconvolve(fftconvolve(Iy*Iy, bfilt, 'valid'), bfilt.T, 'valid')
            Mxy = fftconvolve(fftconvolve(Ix*Iy, bfilt, 'valid'), bfilt.T, 'valid');

            term1   = (Mxx + Myy)/2
            term2   = (term1**2 - (Mxx*Myy - Mxy**2)) ** .5


            eps = 2.2204e-16

            ev1 = term1 + term2
            ev2 = term1 - term2
            ori = np.arctan2(Mxx-ev2, Mxy)-pi/2
            energy = ev1 + ev2
            orientedness  =  ((ev1-ev2)/(ev1+ev2+eps)) **2

            height, width = orientedness.shape[:2]

            rad = .9*height/2
            x = np.linspace(-width/2, width/2-1, width)
            y = np.linspace(-height/2, height/2-1, height)
            xv, yv = np.meshgrid(x, y)
            ind_central = (xv**2 + yv**2)**.5 > rad

            ind1 = energy < max(np.nanquantile(energy[:],.68), 1e-4)
            ind2 = orientedness < .8
            ori_thresholded = ori
            ori_thresholded[ind1] = np.nan
            ori_thresholded[ind2] = np.nan
            ori_thresholded[ind_central] = np.nan

            ori_deg = ori_thresholded.flatten() * (180/pi)

            h = np.histogram(ori_deg, bins=35,range=(-87.5,87.5))

            ind_circ_lower = ori_deg < -87.5
            ind_circ_upper = ori_deg > 87.5

            fs = list(h[0])
            fs.insert(0, np.sum(ori_deg[ind_circ_lower]/ori_deg[ind_circ_lower]) + np.sum(ori_deg[ind_circ_upper]/ori_deg[ind_circ_upper]))
            fs.append(np.sum(ori_deg[ind_circ_lower]/ori_deg[ind_circ_lower]) + np.sum(ori_deg[ind_circ_upper]/ori_deg[ind_circ_upper]))
            fs = np.array(fs)
            curr_or_n.append(int(np.sum(fs)))
            all_hists.append(fs)
        # if there are nan values, then don't append
        if np.any(np.isnan(np.array(all_hists))):
          continue
        n_ors.append(curr_or_n)
        name, _ = os.path.splitext(files[i])
        out_path = os.path.join(save_folder, name + '.npy')
        np.save(out_path, np.array(all_hists))
    except:
        continue

# Serializing json
json_object = json.dumps({"N": n_ors}, indent=4)

with open(os.path.join(save_folder, 'n_ors.json'), "w") as outfile:
    outfile.write(json_object)