Commit 3934ca0c authored by Martin Reinecke's avatar Martin Reinecke

recursive algorithm

parent a141d66e
......@@ -21,7 +21,7 @@ if __name__ == "__main__":
ends = np.random.random((2, n))
return starts, ends
nlos = 1000
nlos = 100
starts, ends = make_random_los(nlos)
R = ift.library.LOSResponse(s_space, starts=starts, ends=ends)
......
......@@ -9,94 +9,61 @@ from ..field import Field
from .. import dobj
def _dist_from_LOS(start, end, pos):
start = start.reshape((-1,1))
end = end.reshape((-1,1))
dir = end-start
dlen = np.sqrt(np.sum(dir*dir))
dir *= 1./dlen
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<dlen, d1, d3))
def _comp_traverse(start, end, shp, dist):
ndim = start.shape[0]
class Pixellister(object):
def __init__(self, start, end):
self._start = start
self._end = end
self._dir = end - start
self._dirlen = np.sqrt(np.sum(self._dir*self._dir))
self._dir *= 1./self._dirlen
def _losdist(self, pos):
t1 = pos - self._start
dotpr = np.sum(self._dir*t1)
if dotpr < 0:
return np.sqrt(np.sum(t1*t1))
elif dotpr <= self._dirlen:
t2 = self._start + self._dir*dotpr
d1 = pos-t2
return np.sqrt(np.sum(d1*d1))
else:
d3 = pos-self._end
return np.sqrt(np.sum(d3*d3))
def build_pixellist(self, lo, hi, pixels):
Rmax=5.
diff = hi-lo
if np.all(diff==1):
dist = self._losdist(lo)
if dist < Rmax:
pixels.append((lo, dist))
return
mid = (lo+hi-1.)/2.
t1 = mid-lo
if self._losdist(mid) < np.sqrt(np.sum(t1*t1)) + Rmax:
dim = np.argmax(diff)
imid = lo[dim]+diff[dim]//2
lmid = lo.copy()
lmid[dim] = imid
hmid = hi.copy()
hmid[dim] = imid
self.build_pixellist(lo, hmid, pixels)
self.build_pixellist(lmid, hi, pixels)
def _comp_traverse(start, end, shp):
nlos = start.shape[1]
inc = np.full(len(shp), 1)
for i in range(-2, -len(shp)-1, -1):
inc[i] = inc[i+1]*shp[i+1]
combo = [0]
for i in range(len(shp)):
tmp = [x+inc[i] for x in combo]
tmp += [x-inc[i] for x in combo]
combo += tmp
pmax = np.array(shp)
out = [None]*nlos
for i in range(nlos):
dir = end[:, i]-start[:, i]
dirx = np.where(dir == 0., 1e-12, dir)
d0 = np.where(dir == 0., ((start[:, i] > 0)-0.5)*1e12,
-start[:, i]/dirx)
d1 = np.where(dir == 0., ((start[:, i] < pmax)-0.5)*-1e12,
(pmax-start[:, i])/dirx)
(dmin, dmax) = (np.minimum(d0, d1), np.maximum(d0, d1))
dmin = dmin.max()
dmax = dmax.min()
dmin = np.maximum(0., dmin)
dmax = np.minimum(1., dmax)
dmax = np.maximum(dmin, dmax)
# hack: move away from potential grid crossings
dmin += 1e-7
dmax -= 1e-7
if dmin > dmax: # no intersection
out[i] = (np.full(0, 0), np.full(0, 0.))
continue
# determine coordinates of first cell crossing
c_first = np.ceil(start[:, i]+dir*dmin)
c_first = np.where(dir > 0., c_first, c_first-1.)
c_first = (c_first-start[:, i])/dirx
pos1 = np.asarray((start[:, i]+dmin*dir), dtype=np.int)
pos1 = np.sum(pos1*inc)
cdist = np.empty(0, dtype=np.float64)
add = np.empty(0, dtype=np.int)
for j in range(ndim):
if dir[j] != 0:
step = inc[j] if dir[j] > 0 else -inc[j]
tmp = np.arange(start=c_first[j], stop=dmax,
step=abs(1./dir[j]))
cdist = np.append(cdist, tmp)
add = np.append(add, np.full(len(tmp), step))
idx = np.argsort(cdist)
cdist = cdist[idx]
add = add[idx]
cdist = np.append(np.full(1, dmin), cdist)
cdist = np.append(cdist, np.full(1, dmax))
corfac = np.linalg.norm(dir*dist)
cdist *= corfac
wgt = np.diff(cdist)
mdist = 0.5*(cdist[:-1]+cdist[1:])
add = np.append(pos1, add)
add = np.cumsum(add)
zone = set()
for x in combo:
zone.update(add+x)
npix = np.prod(shp)
zone = [z for z in zone if z>=0 and z<npix]
fullidx = np.unravel_index(zone, shp)
fullidx=np.array(fullidx)
xdist = _dist_from_LOS(start[:,i], end[:,i], fullidx)
xxwgt = np.exp(-9*xdist*xdist)
twgt = np.sum(xxwgt)
xxwgt /= twgt
out[i] = (np.array(zone), np.array(xxwgt))
print i
pixels = []
pl = Pixellister(start[:,i], end[:,i])
pl.build_pixellist(np.array(shp)*0, np.array(shp), pixels)
pix1 = [np.sum(inc*p[0]) for p in pixels]
wgt = [np.exp(-p[1]*p[1]) for p in pixels]
out[i] = (pix1, wgt)
return out
......@@ -125,7 +92,6 @@ class LOSResponse(LinearOperator):
every task).
"""
def __init__(self, domain, starts, ends):
super(LOSResponse, self).__init__()
self._domain = DomainTuple.make(domain)
......@@ -159,8 +125,7 @@ class LOSResponse(LinearOperator):
# 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[0]))
self._local_shape)
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