Commit 2c7b4895 authored by Martin Reinecke's avatar Martin Reinecke

tweaking

parent 2dbf90a7
Pipeline #28344 passed with stages
in 2 minutes and 37 seconds
......@@ -17,8 +17,8 @@ if __name__ == "__main__":
ss = ht(sh)
def make_random_los(n):
starts = np.random.random((2, n))
ends = np.random.random((2, n))
starts = np.random.random((2, n))*2-0.5
ends = np.random.random((2, n))*2-0.5
return starts, ends
nlos = 1000
......@@ -28,7 +28,7 @@ if __name__ == "__main__":
Rh = R*ht
signal_to_noise = 10
signal_to_noise = 100
noise_amplitude = Rh(sh).val.std()/signal_to_noise
N = ift.ScalingOperator(noise_amplitude**2, Rh.target)
n = N.draw_sample()
......
......@@ -7,30 +7,26 @@ from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from .. import dobj
from array import array
class Pixellister(object):
Rmax = np.sqrt(3.)*0.5
def __init__(self, start, end):
def __init__(self, start, end, inc):
self._start = start
self._end = end
self._dir = end - start
self._dirlen = np.linalg.norm(self._dir, axis=0)
self._dir *= 1./self._dirlen.reshape((1, -1))
self._inc = inc
def _losdist(self, idx, pos):
pos2 = pos.reshape((-1, 1))
sidx = self._start[:, idx]
eidx = self._end[:, idx]
didx = self._dir[:, idx]
dlidx = self._dirlen[idx]
t1 = pos2 - sidx
t1 = pos2 - self._start[:, idx]
dotpr = np.sum(didx*t1, axis=0)
return np.where(
dotpr < 0, np.linalg.norm(t1, axis=0),
np.where(dotpr > dlidx, np.linalg.norm(pos2-eidx, axis=0),
np.linalg.norm(pos2 - sidx - didx*dotpr, axis=0)))
dotpr = np.maximum(0., np.minimum (self._dirlen[idx], dotpr))
return np.linalg.norm(t1 - didx*dotpr, axis=0)
def build_pixellist(self, idx, lo, hi, pixels):
diff = hi-lo
......@@ -39,7 +35,8 @@ class Pixellister(object):
if np.all(diff == 1):
for i in range(len(dist)):
if dist[i] < self.Rmax:
pixels[idx[i]].append((lo, dist[i]))
pixels[0][idx[i]].append(np.sum(self._inc*lo))
pixels[1][idx[i]].append(dist[i])
else:
t1 = mid-lo
idx2 = idx[dist < np.linalg.norm(t1) + self.Rmax]
......@@ -57,18 +54,15 @@ class Pixellister(object):
def _comp_traverse(start, end, shp):
nlos = start.shape[1]
ndim = start.shape[0]
inc = np.full(len(shp), 1)
for i in range(-2, -len(shp)-1, -1):
inc[i] = inc[i+1]*shp[i+1]
out = [None]*nlos
pl = Pixellister(start, end)
pixels = [[] for _ in range(nlos)]
pl = Pixellister(start, end, inc)
pixels = [[array("l") for _ in range(nlos)], [array("d") for _ in range(nlos)]]
pl.build_pixellist(np.arange(nlos), np.array(shp)*0, np.array(shp), pixels)
for i in range(nlos):
pix1 = np.array([np.sum(inc*p[0]) for p in pixels[i]])
wgt = np.array([np.exp(-p[1]*p[1]) for p in pixels[i]])
out[i] = (pix1, wgt)
return out
return [(np.array(p), np.exp(-(np.array(w)**2)))
for p, w in zip(pixels[0], pixels[1])]
class LOSResponse(LinearOperator):
......
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