Commit b76e0847 authored by Martin Reinecke's avatar Martin Reinecke

preparations

parent ca2ceecf
......@@ -24,9 +24,7 @@ if __name__ == "__main__":
nlos = 1000
starts, ends = make_random_los(nlos)
R = ift.library.LOSResponse(
s_space, starts=starts, ends=ends,
sigmas_up=np.full(nlos, 0.1), sigmas_low=np.full(nlos, 0.1))
R = ift.library.LOSResponse(s_space, starts=starts, ends=ends)
Rh = R*ht
......
import numpy as np
from scipy.special import erfc
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import aslinearoperator
from ..operators.linear_operator import LinearOperator
......@@ -10,11 +9,22 @@ from ..field import Field
from .. import dobj
def _gaussian_error_function(x):
return 0.5*erfc(x*np.sqrt(2.))
def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
def _dist_from_LOS(start, end, pos):
start = start.reshape((-1,1))
end = end.reshape((-1,1))
dir = end-start
dir *= 1./np.sqrt(np.sum(dir*dir))
t1 = pos - start
dotpr = np.sum(dir*t1,axis=0)
t2 = start + dir*dotpr.reshape((1,-1))
d1 = pos-t2
d1 = np.sqrt(np.sum(d1*d1,axis=0))
d2 = np.sqrt(np.sum(t1*t1,axis=0))
d3 = pos-end
d3 = np.sqrt(np.sum(d3*d3,axis=0))
return np.where(dotpr<0, d2, np.where(dotpr<1, d1,d3))
def _comp_traverse(start, end, shp, dist):
ndim = start.shape[0]
nlos = start.shape[1]
inc = np.full(len(shp), 1)
......@@ -67,24 +77,12 @@ def _comp_traverse(start, end, shp, dist, lo, mid, hi, erf):
cdist *= corfac
wgt = np.diff(cdist)
mdist = 0.5*(cdist[:-1]+cdist[1:])
wgt = apply_erf(wgt, mdist, lo[i], mid[i], hi[i], erf)
add = np.append(pos1, add)
add = np.cumsum(add)
out[i] = (add, wgt)
return out
def apply_erf(wgt, dist, lo, mid, hi, erf):
wgt = wgt.copy()
mask = dist > hi
wgt[mask] = 0.
mask = (dist > mid) & (dist <= hi)
wgt[mask] *= erf((dist[mask]-mid)/(hi-mid))
mask = (dist <= mid) & (dist > lo)
wgt[mask] *= erf((dist[mask]-mid)/(mid-lo))
return wgt
class LOSResponse(LinearOperator):
"""Line-of-sight response operator
......@@ -95,22 +93,21 @@ class LOSResponse(LinearOperator):
Parameters
----------
domain : RGSpace or DomainTuple
The operator's input domain. This must be a single RGSpace.
The operator's input domain. This must be a single RGSpace, with
all distances the same
starts, ends : numpy.ndarray(float) with two dimensions
Arrays containing the start and end points of the individual lines
of sight. The first dimension must have as many entries as `domain`
has dimensions. The second dimensions must be identical for both arrays
and indicated the total number of lines of sight.
sigmas_low, sigmas_up : numpy.ndarray(float) (optional)
For expert use. If unsure, leave blank.
Notes
-----
`starts, `ends`, `sigmas_low`, and `sigmas_up` have to be identical on
`starts and `ends` have to be identical on
every calling MPI task (i.e. the full LOS information has to be provided on
every task).
"""
def __init__(self, domain, starts, ends, sigmas_low=None, sigmas_up=None):
def __init__(self, domain, starts, ends):
super(LOSResponse, self).__init__()
self._domain = DomainTuple.make(domain)
......@@ -119,20 +116,15 @@ class LOSResponse(LinearOperator):
(len(self._domain) != 1)):
raise TypeError("The domain must be exactly one RGSpace instance.")
if not np.all(self.domain[0].distances == self.domain[0].distances[0]):
raise ValueError("domain distances must all be equal")
ndim = len(self.domain[0].shape)
starts = np.array(starts)
nlos = starts.shape[1]
ends = np.array(ends)
if sigmas_low is None:
sigmas_low = np.zeros(nlos, dtype=np.float32)
if sigmas_up is None:
sigmas_up = np.zeros(nlos, dtype=np.float32)
sigmas_low = np.array(sigmas_low)
sigmas_up = np.array(sigmas_up)
assert starts.shape[0] == ndim, "dimension mismatch"
assert nlos == sigmas_low.shape[0], "dimension mismatch"
assert starts.shape == ends.shape, "dimension mismatch"
assert sigmas_low.shape == sigmas_up.shape, "dimension mismatch"
self._local_shape = dobj.local_shape(self.domain[0].shape)
local_zero_point = (np.array(
......@@ -142,19 +134,16 @@ class LOSResponse(LinearOperator):
diffs = ends-starts
difflen = np.linalg.norm(diffs, axis=0)
diffs /= difflen
real_ends = ends + sigmas_up*diffs
lzp = local_zero_point.reshape((-1, 1))
dist = np.array(self.domain[0].distances).reshape((-1, 1))
localized_pixel_starts = (starts-lzp)/dist + 0.5
localized_pixel_ends = (real_ends-lzp)/dist + 0.5
localized_pixel_ends = (ends-lzp)/dist + 0.5
# get the shape of the local data slice
w_i = _comp_traverse(localized_pixel_starts,
localized_pixel_ends,
self._local_shape,
np.array(self.domain[0].distances),
difflen-sigmas_low, difflen, difflen+sigmas_up,
_gaussian_error_function)
np.array(self.domain[0].distances))
boxsz = 16
nlos = len(w_i)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment