-
Notifications
You must be signed in to change notification settings - Fork 12
/
astrohog3d.py
207 lines (152 loc) · 7.34 KB
/
astrohog3d.py
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
#!/usr/bin/env python
#
# This file is part of AstroHOG
#
# CONTACT: soler[AT]mpia.de
# Copyright (C) 2013-2019 Juan Diego Soler
#
#------------------------------------------------------------------------------;
import sys
import numpy as np
from astropy.convolution import convolve_fft
from astropy.convolution import Gaussian2DKernel
from scipy import ndimage
import pycircstat as circ
from nose.tools import assert_equal, assert_true
import matplotlib.pyplot as plt
import collections
import multiprocessing
from astrohog2d import *
from statests import *
from tqdm import tqdm
# ----------------------------------------------------------------------
def calculatexi(cosangles, nbins=20, thres=0.2):
hist, bin_edges = np.histogram(np.abs(cosangles), bins=nbins, range=[0.0,1.0])
bin_centres=0.5*(bin_edges[0:nbins]+bin_edges[1:nbins+1])
para=np.sum(hist[(bin_centres>=(1.0-thres)).nonzero()])
perp=np.sum(hist[(bin_centres<=thres).nonzero()])
xi=(para-perp)/(para+perp)
return xi
# --------------------------------------------------------------------------------------------------------------------------------
def HOGcorr_cubeLITE(cube1, cube2, pxsz=1., ksz=1., mode='nearest', mask1=0., mask2=0., gradthres1=0., gradthres2=0., weights=None, weightbygradnorm=False):
""" Calculates the spatial correlation between cube1 and cube2 using the HOG method
Parameters
----------
cube1 : array corresponding to the first cube to be compared
cube2 : array corresponding to the second cube to be compared
pxsz : pixel size
ksz : size of the derivative kernel in the pixel size units
Returns
-------
hogcorr :
corrframe :
Examples
--------
"""
assert cube2.shape == cube1.shape, "Dimensions of cube2 and cube1 must match"
sz1=np.shape(cube1)
if weights is None:
weights=np.ones(sz1)
if (np.size(weights)==1):
uniweights=weights
weights=uniweights*np.ones(sz1)
assert weights.shape == cube1.shape, "Dimensions of weights and ima1 must match"
pxksz=(ksz/(2*np.sqrt(2.*np.log(2.))))/pxsz #gaussian_filter takes sigma instead of FWHM as input
scube1=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[0,0,0], mode=mode)
scube2=ndimage.filters.gaussian_filter(cube2, [pxksz, pxksz, pxksz], order=[0,0,0], mode=mode)
dI1dx=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[0,0,1], mode=mode)
dI1dy=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[0,1,0], mode=mode)
dI1dz=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[1,0,0], mode=mode)
dI2dx=ndimage.filters.gaussian_filter(cube2, [pxksz, pxksz, pxksz], order=[0,0,1], mode=mode)
dI2dy=ndimage.filters.gaussian_filter(cube2, [pxksz, pxksz, pxksz], order=[0,1,0], mode=mode)
dI2dz=ndimage.filters.gaussian_filter(cube2, [pxksz, pxksz, pxksz], order=[1,0,0], mode=mode)
normGrad1=np.sqrt(dI1dx**2+dI1dy**2+dI1dz**2)
normGrad2=np.sqrt(dI2dx**2+dI2dy**2+dI2dz**2)
# Calculation of the relative orientation angles
#crossp=np.array([dI1dy*dI2dz-dI1dz*dI2dy,dI1dz*dI2dx-dI1dx*dI2dz,dI1dx*dI2dy-dI1dy*dI2dx])
#crossp=np.sqrt(np.sum(crossp**2))
dotp=dI1dx*dI2dx+dI1dy*dI2dy+dI1dz*dI2dz
cosphi=dotp/(normGrad1*normGrad2)
phi=np.arccos(cosphi)
# Excluding null gradients
bad=np.logical_or(normGrad1 <= gradthres1, normGrad2 <= gradthres2).nonzero()
phi[bad]=np.nan
# Excluding masked gradients
if np.array_equal(np.shape(cube1), np.shape(mask1)):
m1bad=(mask1 < 1.).nonzero()
phi[m1bad]=np.nan
if np.array_equal(np.shape(cube2), np.shape(mask2)):
m2bad=(mask2 < 1.).nonzero()
phi[m2bad]=np.nan
good=np.isfinite(phi).nonzero()
if (weightbygradnorm):
weights=normGrad1*normGrad2
Zx, s_Zx, meanPhi = HOG_PRS(phi[good])
rvl=circ.descriptive.resultant_vector_length(2.*phi[good], w=weights[good])
can=circ.descriptive.mean(2.*phi[good], w=weights[good])/2.
pz, Z = circ.tests.rayleigh(2.*phi[good], w=weights[good])
pv, V = circ.tests.vtest(2.*phi[good], 0., w=weights[good])
myV, s_myV, meanphi = HOG_PRS(2.*phi[good])
am=HOG_AM(phi[good])
pear, peap = stats.pearsonr(scube1[good], scube2[good])
ngood=np.size(np.isfinite(phi.ravel()).nonzero())
xi=calculatexi(cosphi[good])
#circstats=[rvl, Z, V, pz, pv, myV, s_myV, meanphi, am, pear, ngood, ssimv, msev]
circstats = {'r': rvl, 'Z': Z, 'V': V, 'meanphi': meanphi, 'xi': xi}
corrframe=phi
return circstats, corrframe, scube1, scube2
# -----------------------------------------------------------------------------------------------------------------------
def HOGcorr_cubeANDvecLITE(cube1, vec, pxsz=1., ksz=1., mode='nearest', mask1=0., mask2=0., gradthres1=0., gradthres2=0., weights=None):
""" Calculates the correlation relative orientation between cube1 and the vector field
Parameters
----------
cube1 : array corresponding to the scale field cube
vec : array corresponding to the vector field [v_0,v_1,v_2], where v_i (i=0,1,2) corresponds to the i-th index of cube1
pxsz : pixel size
ksz : size of the derivative kernel in the pixel size units
Returns
-------
hogcorr :
corrframe :
Examples
--------
"""
assert vec1[0].shape == cube1.shape, "Dimensions of vec[0] and cube1 must match"
assert vec1[1].shape == cube1.shape, "Dimensions of vec[1] and cube1 must match"
assert vec1[2].shape == cube1.shape, "Dimensions of vec[2] and cube1 must match"
scube1=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[0,0,0], mode=mode)
dcube1d0=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[1,0,0], mode=mode)
dcube1d1=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[0,1,0], mode=mode)
dcube1d2=ndimage.filters.gaussian_filter(cube1, [pxksz, pxksz, pxksz], order=[0,0,1], mode=mode)
crossp=np.array([dcube1d1*vec[2]-dcube1d2*vec[1],dcube1d2*vec[0]-dcube1d0*vec[2],dcube1d0*vec[1]-dcube1d1*vec[0]])
normcrossp=np.sqrt(np.sum(crossp**2))
dotp=dcube1d0*vec[0]+dcube1d1*vec[1]+dcube1d2*vec[2]
tempphi=np.arctan2(normcrossp, dotp)
phi=np.arctan(np.tan(tempphi))
# Excluding null gradients
normGrad1=np.sqrt(dcube1d0**2+dcube1d1**2+dcube1d2**2)
normvec=np.sqrt(vec[0]**2+vec[1]**2+vec[2]**2)
bad=np.logical_or(normGrad1 <= gradthres1, normvec <= 0.0).nonzero()
phi[bad]=np.nan
# Excluding masked gradients
if np.array_equal(np.shape(cube1), np.shape(mask1)):
m1bad=(mask1 < 1.).nonzero()
phi[m1bad]=np.nan
if np.array_equal(np.shape(cube2), np.shape(mask2)):
m2bad=(mask2 < 1.).nonzero()
phi[m2bad]=np.nan
good=np.isfinite(phi).nonzero()
Zx, s_Zx, meanPhi = HOG_PRS(phi[good])
rvl=circ.descriptive.resultant_vector_length(2.*phi[good], w=weights[good])
can=circ.descriptive.mean(2.*phi[good], w=weights[good])/2.
pz, Z = circ.tests.rayleigh(2.*phi[good], w=weights[good])
pv, V = circ.tests.vtest(2.*phi[good], 0., w=weights[good])
myV, s_myV, meanphi = HOG_PRS(2.*phi[good])
am=HOG_AM(phi[good])
pear, peap = stats.pearsonr(scube1[good], scube2[good])
ngood=np.size(good)
ssimv=np.nan #ssim(sima1[good], sima2[good])
msev =np.nan #mse(sima1[good], sima2[good])
circstats=[rvl, Z, V, pz, pv, myV, s_myV, meanphi, am, pear, ngood, ssimv, msev]
corrframe=phi
return circstats, corrframe, scube1