Skip to content

Commit

Permalink
no-op rearrangement, variable renaming, comments
Browse files Browse the repository at this point in the history
  • Loading branch information
Dima Kogan authored and Dima Kogan committed Aug 23, 2020
1 parent ac45347 commit 6039abe
Showing 1 changed file with 12 additions and 17 deletions.
29 changes: 12 additions & 17 deletions mrcal/utils.py
Expand Up @@ -3312,11 +3312,12 @@ def intrinsics_implied_Rt10(q0, v0, v1,
weights = nps.clump(weights, n = len(weights.shape)-2)

i_nan_v0 = ~np.isfinite(v0)
i_nan_v1 = ~np.isfinite(v1)
v0[i_nan_v0] = 0.
weights[i_nan_v0[...,0]] = 0.0
weights[i_nan_v0[...,1]] = 0.0
weights[i_nan_v0[...,2]] = 0.0

i_nan_v1 = ~np.isfinite(v1)
v1[i_nan_v1] = 0.
weights[..., i_nan_v1[...,0]] = 0.0
weights[..., i_nan_v1[...,1]] = 0.0
Expand All @@ -3328,30 +3329,26 @@ def intrinsics_implied_Rt10(q0, v0, v1,
if np.count_nonzero(i)<3:
raise Exception("Focus region contained too few points")

v0_cut = v0 [...,i, :]
v1_cut = v1 [ i, :]
wcut = weights[...,i ]

if not atinfinity:
p0 = v0
p0_cut = v0_cut
v0_cut = v0 [...,i, :]
v1_cut = v1 [ i, :]
weights = weights[...,i ]

def residual_jacobian_rt(rt):

# rtp0 has shape (...,N,3)
rtp0, drtp0_dr, drtp0_dt, _ = \
mrcal.transform_point_rt(rt, p0_cut,
mrcal.transform_point_rt(rt, v0_cut,
get_gradients = True)

# shape (...,N,3,6)
drtp0_drt = nps.glue( drtp0_dr, drtp0_dt, axis=-1)

# inner(a,b) ~ cos(x) ~ 1 - x^2/2
# inner(a,b)/(mag(a)*mag(b)) ~ cos(x) ~ 1 - x^2/2
# Each of these has shape (...,N)
mag_rtp0 = nps.mag(rtp0)
inner = nps.inner(rtp0, v1_cut)
th2 = 2.* (1.0 - inner / mag_rtp0)
x = th2 * wcut
x = th2 * weights

# shape (...,N,6)
dmag_rtp0_drt = nps.matmult( nps.dummy(rtp0, -2), # shape (...,N,1,3)
Expand All @@ -3371,7 +3368,7 @@ def residual_jacobian_rt(rt):
(nps.dummy(inner, -1) * dmag_rtp0_drt - \
nps.dummy(mag_rtp0, -1) * dinner_drt) / \
nps.dummy(mag_rtp0*mag_rtp0, -1) * \
nps.dummy(wcut,-1)
nps.dummy(weights,-1)
return x.ravel(), nps.clump(J, n=len(J.shape)-1)


Expand All @@ -3383,19 +3380,19 @@ def residual_jacobian_r(r):
mrcal.rotate_point_r(r, v0_cut,
get_gradients = True)

# inner(a,b) ~ cos(x) ~ 1 - x^2/2
# inner(a,b)/(mag(a)*mag(b)) ~ cos(x) ~ 1 - x^2/2
# Each of these has shape (N)
inner = nps.inner(rv0, v1_cut)
th2 = 2.* (1.0 - inner)
x = th2 * wcut
x = th2 * weights

# shape (N,3)
dinner_dr = nps.matmult( nps.dummy(v1_cut, -2), # shape (N,1,3)
drv0_dr # shape (N,3,3)
# matmult has shape (N,1,3)
)[:,0,:]

J = -2. * dinner_dr * nps.dummy(wcut,-1)
J = -2. * dinner_dr * nps.dummy(weights,-1)
return x, J


Expand All @@ -3414,8 +3411,6 @@ def jacobian(rt, f):

# # gradient check
# import gnuplotlib as gp
# p0 = v0
# p0_cut = v0_cut
# rt0 = np.random.random(6)*1e-3
# x0,J0 = residual_jacobian_rt(rt0)
# drt = np.random.random(6)*1e-7
Expand Down

0 comments on commit 6039abe

Please sign in to comment.