diff --git a/Dockerfile b/Dockerfile index dbeaf3e7468d24158bb63ffaa36908d6feda16e8..157129f7c223252d23f973daa26544175d660beb 100644 --- a/Dockerfile +++ b/Dockerfile @@ -10,7 +10,7 @@ RUN apt-get update && apt-get install -y \ # Testing dependencies python3-pytest-cov jupyter \ # Optional NIFTy dependencies - libfftw3-dev python3-mpi4py python3-matplotlib \ + libfftw3-dev python3-mpi4py python3-matplotlib python3-pynfft \ # more optional NIFTy dependencies && pip3 install pyfftw \ && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \ diff --git a/nifty5/__init__.py b/nifty5/__init__.py index 52d8311475f1f837f9951715b3f164585a9d62bb..a6bd52eafbc76a33654c92478270e52f1dd2ebff 100644 --- a/nifty5/__init__.py +++ b/nifty5/__init__.py @@ -85,6 +85,7 @@ from .library.wiener_filter_curvature import WienerFilterCurvature from .library.correlated_fields import CorrelatedField, MfCorrelatedField from .library.adjust_variances import (make_adjust_variances_hamiltonian, do_adjust_variances) +from .library.nfft import NFFT from . import extra diff --git a/nifty5/library/nfft.py b/nifty5/library/nfft.py new file mode 100644 index 0000000000000000000000000000000000000000..c96866d6ee948c0fbda23a7720c53b792059e5d9 --- /dev/null +++ b/nifty5/library/nfft.py @@ -0,0 +1,73 @@ +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <http://www.gnu.org/licenses/>. +# +# Copyright(C) 2018-2019 Max-Planck-Society +# +# Resolve is being developed at the Max-Planck-Institut fuer Astrophysik. + +import numpy as np + +import nifty5 as ift + + +class NFFT(ift.LinearOperator): + """Performs a non-equidistant Fourier transform, i.e. a Fourier transform + followed by a degridding operation. + + Parameters + ---------- + domain : RGSpace + Domain of the operator. It has to be two-dimensional and have shape + `(2N, 2N)`. The coordinates of the lower left pixel of the dirty image + are `(-N,-N)`, and of the upper right pixel `(N-1,N-1)`. + uv : numpy.ndarray + 2D numpy array of type float64 and shape (M,2), where M is the number + of measurements. uv[i,0] and uv[i,1] contain the u and v coordinates + of measurement #i, respectively. All coordinates must lie in the range + `[-0.5; 0,5[`. + """ + def __init__(self, domain, uv): + from pynfft.nfft import NFFT + npix = domain.shape[0] + assert npix == domain.shape[1] + assert len(domain.shape) == 2 + assert type(npix) == int, "npix must be integer" + assert npix > 0 and ( + npix % 2) == 0, "npix must be an even, positive integer" + assert isinstance(uv, np.ndarray), "uv must be a Numpy array" + assert uv.dtype == np.float64, "uv must be an array of float64" + assert uv.ndim == 2, "uv must be a 2D array" + assert uv.shape[0] > 0, "at least one point needed" + assert uv.shape[1] == 2, "the second dimension of uv must be 2" + assert np.all(uv >= -0.5) and np.all(uv <= 0.5),\ + "all coordinates must lie between -0.5 and 0.5" + + self._domain = ift.DomainTuple.make(domain) + self._target = ift.DomainTuple.make( + ift.UnstructuredDomain(uv.shape[0])) + self._capability = self.TIMES | self.ADJOINT_TIMES + + self.npt = uv.shape[0] + self.plan = NFFT(self.domain.shape, self.npt, m=6) + self.plan.x = uv + self.plan.precompute() + + def apply(self, x, mode): + self._check_input(x, mode) + if mode == self.TIMES: + self.plan.f_hat = x.to_global_data() + res = self.plan.trafo().copy() + else: + self.plan.f = x.to_global_data() + res = self.plan.adjoint().copy() + return ift.Field.from_global_data(self._tgt(mode), res) diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py index 718498d021a50721a96e237784cfc5322dbca7a0..6727b47568a637cb743d4f05bea5fef8cd4e99b1 100644 --- a/test/test_operators/test_adjoint.py +++ b/test/test_operators/test_adjoint.py @@ -279,3 +279,10 @@ def testValueInserter(sp, seed): ind.append(np.random.randint(0, ss-1)) op = ift.ValueInserter(sp, ind) ift.extra.consistency_check(op) + + +def testNFFT(): + dom = ift.RGSpace(2*(16,)) + uv = np.array([[.2, .4], [-.22, .452]]) + op = ift.NFFT(dom, uv) + ift.extra.consistency_check(op)