From 8272336759f1d71fe1d9934eb8d4263bce604909 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Mon, 25 Mar 2024 17:29:42 -0500 Subject: [PATCH 01/21] add fourier Filter --- src/tomocupy/config.py | 18 +++- src/tomocupy/processing/proc_functions.py | 3 + src/tomocupy/processing/retrieve_phase.py | 107 ++++++++++++++++++++++ 3 files changed, 127 insertions(+), 1 deletion(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index 6a4a055..f36de1b 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -224,7 +224,7 @@ def default_parameter(func, param): 'default': 'none', 'type': str, 'help': "Phase retrieval correction method", - 'choices': ['none', 'paganin', 'Gpaganin']}, + 'choices': ['none', 'paganin', 'Gpaganin','FF']}, 'energy': { 'default': 0, 'type': float, @@ -249,6 +249,22 @@ def default_parameter(func, param): 'type': utils.positive_int, 'default': 1, 'help': "Padding with extra slices in z for phase-retrieval filtering"}, + 'FFratio': { + 'default': 100, + 'type': float, + 'help': "Shape of the Fourier Filter window"}, + 'FFdim': { + 'default': 2, + 'type': float, + 'help': "1 for sinograms, 2 for projections"}, + 'FFlog': { + 'default': 0, + 'type': int, + 'help': "0/1"}, + 'FFpad': { + 'default': 150, + 'type': float, + 'help': "Padding size for the Fourier Filter"}, } SECTIONS['rotate-proj'] = { diff --git a/src/tomocupy/processing/proc_functions.py b/src/tomocupy/processing/proc_functions.py index 1807a5d..b8b5f8a 100644 --- a/src/tomocupy/processing/proc_functions.py +++ b/src/tomocupy/processing/proc_functions.py @@ -171,6 +171,9 @@ def proc_proj(self, data, st=None, end=None, res=None): data, args.pixel_size*1e-4, args.propagation_distance/10, args.energy, args.retrieve_phase_alpha, args.retrieve_phase_method, args.retrieve_phase_delta_beta, args.retrieve_phase_W*1e-4) # units adjusted based on the tomopy implementation + if args.retrieve_phase_method == 'FF': + data[:] = retrieve_phase.fresnel_filter( + data, args.FFratio, args.FFdim, window=None, pad=args.FFpad, apply_log=args.FFlog) if args.rotate_proj_angle != 0: data[:] = self.rotate_proj( data, args.rotate_proj_angle, args.rotate_proj_order) diff --git a/src/tomocupy/processing/retrieve_phase.py b/src/tomocupy/processing/retrieve_phase.py index b789910..db4235c 100644 --- a/src/tomocupy/processing/retrieve_phase.py +++ b/src/tomocupy/processing/retrieve_phase.py @@ -258,3 +258,110 @@ def _reciprocal_coord(pixel_size, num_grid): rc = cp.arange(-n, num_grid, 2, dtype=cp.float32) rc *= 0.5 / (n * pixel_size) return rc + + +def make_fresnel_window(height, width, ratio, dim): + """ + Create a low pass window based on the Fresnel propagator. + It is used to denoise a projection image (dim=2) or a + sinogram image (dim=1). + + Parameters + ---------- + height : int + Image height + width : int + Image width + ratio : float + To define the shape of the window. + dim : {1, 2} + Use "1" if working on a sinogram image and "2" if working on + a projection image. + + Returns + ------- + array_like + 2D array. + """ + ycenter = (height - 1) * 0.5 + xcenter = (width - 1) * 0.5 + if dim == 2: + u = (cp.arange(width) - xcenter) / width + v = (cp.arange(height) - ycenter) / height + u, v = cp.meshgrid(u, v) + window = 1.0 + ratio * (u ** 2 + v ** 2) + else: + u = (cp.arange(width) - xcenter) / width + win1d = 1.0 + ratio * u ** 2 + window = cp.tile(win1d, (height, 1)) + return window + + +def fresnel_filter(data, ratio, dim, window=None, pad=150, apply_log=True): + """ + Apply a low-pass filter based on the Fresnel propagator to an image + (Ref. [1]). It can be used for improving the contrast of an image. + It's simpler than the well-known Paganin's filter (Ref. [2]). + + Parameters + ---------- + mat : array_like + 2D array. Projection image or sinogram image. + ratio : float + Define the shape of the window. Larger is more smoothing. + dim : {1, 2} + Use "1" if working on a sinogram image and "2" if working on + a projection image. + window : array_like, optional + Window for deconvolution. + pad : int + Padding width. + apply_log : bool, optional + Apply the logarithm function to the sinogram before filtering. + + Returns + ------- + array_like + 2D array. Filtered image. + + References + ---------- + [1] : https://doi.org/10.1364/OE.418448 + + [2] : https://tinyurl.com/2f8nv875 + """ + if apply_log: + data = -cp.log(data) + + if dim == 2: + (nrow, ncol, num_jobs) = data.shape + else: + (nrow, num_jobs, ncol) = data.shape + for m in range(num_jobs): + mat = data[:, m, :] if dim != 2 else data[:, :, m] + + if dim == 2: # On projections + if window is None: + window = make_fresnel_window(nrow, ncol, ratio, dim) + mat_pad = cp.pad(mat, pad, mode="edge") + win_pad = cp.pad(window, pad, mode="edge") + + mat_dec = cp.fft.ifft2(cp.fft.fft2(mat_pad) / cp.fft.ifftshift(win_pad)) + mat_dec = cp.real(mat_dec[pad:pad + nrow, pad:pad + ncol]) + data[:, :, m] = mat_dec + else: # On sinograms + if window is None: + window = make_fresnel_window(nrow, ncol, ratio, dim) + mat_pad = cp.pad(mat, ((0, 0), (pad, pad)), mode='edge') + win_pad = cp.pad(window, ((0, 0), (pad, pad)), mode="edge") + mat_fft = cp.fft.fftshift(cp.fft.fft(mat_pad), axes=1) / win_pad + mat_dec = cp.fft.ifft(cp.fft.ifftshift(mat_fft, axes=1)) + mat_dec = cp.real(mat_dec[:, pad:pad + ncol]) + data[:, m, :] = mat_dec + + if apply_log: + data = cp.exp(-data) + + return data + + From 46157cdfcf0e719fa7be8aa3f8faea92f07336e0 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Tue, 16 Apr 2024 12:01:41 -0500 Subject: [PATCH 02/21] bug fix --- src/tomocupy/processing/retrieve_phase.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/tomocupy/processing/retrieve_phase.py b/src/tomocupy/processing/retrieve_phase.py index db4235c..0eff09c 100644 --- a/src/tomocupy/processing/retrieve_phase.py +++ b/src/tomocupy/processing/retrieve_phase.py @@ -337,23 +337,26 @@ def fresnel_filter(data, ratio, dim, window=None, pad=150, apply_log=True): (nrow, ncol, num_jobs) = data.shape else: (nrow, num_jobs, ncol) = data.shape + #create FF window + if window is None: + if dim == 2: # On projections + window = make_fresnel_window(nrow, ncol, ratio, dim) + else: + window = make_fresnel_window(nrow, ncol, ratio, dim) + for m in range(num_jobs): mat = data[:, m, :] if dim != 2 else data[:, :, m] - - if dim == 2: # On projections - if window is None: - window = make_fresnel_window(nrow, ncol, ratio, dim) - mat_pad = cp.pad(mat, pad, mode="edge") + if dim == 2: + mat_pad = cp.pad(data[:, :, m], pad, mode="edge") win_pad = cp.pad(window, pad, mode="edge") - + #win_pad = cp.repeat(win_pad[:,:,cp.newaxis], num_jobs,axis=2) mat_dec = cp.fft.ifft2(cp.fft.fft2(mat_pad) / cp.fft.ifftshift(win_pad)) mat_dec = cp.real(mat_dec[pad:pad + nrow, pad:pad + ncol]) data[:, :, m] = mat_dec else: # On sinograms - if window is None: - window = make_fresnel_window(nrow, ncol, ratio, dim) - mat_pad = cp.pad(mat, ((0, 0), (pad, pad)), mode='edge') + mat_pad = cp.pad(data[:, m, :], ((0, 0), (pad, pad)), mode='edge') win_pad = cp.pad(window, ((0, 0), (pad, pad)), mode="edge") + #win_pad = cp.repeat(win_pad[:,cp.newaxis,:], num_jobs,axis=1) mat_fft = cp.fft.fftshift(cp.fft.fft(mat_pad), axes=1) / win_pad mat_dec = cp.fft.ifft(cp.fft.ifftshift(mat_fft, axes=1)) mat_dec = cp.real(mat_dec[:, pad:pad + ncol]) From af13b220ba45d9aeae2a77f2c3beb629a2267fab Mon Sep 17 00:00:00 2001 From: mittoalb Date: Tue, 10 Sep 2024 16:06:29 -0500 Subject: [PATCH 03/21] add phase method --- src/tomocupy/config.py | 4 +- src/tomocupy/processing/proc_functions.py | 4 ++ src/tomocupy/processing/retrieve_phase.py | 45 +++++++++++++++++++++++ 3 files changed, 51 insertions(+), 2 deletions(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index f36de1b..0d8a80f 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -224,7 +224,7 @@ def default_parameter(func, param): 'default': 'none', 'type': str, 'help': "Phase retrieval correction method", - 'choices': ['none', 'paganin', 'Gpaganin','FF']}, + 'choices': ['none', 'paganin', 'Gpaganin', 'FF', 'farago']}, 'energy': { 'default': 0, 'type': float, @@ -240,7 +240,7 @@ def default_parameter(func, param): 'retrieve-phase-delta-beta': { 'default': 1500.0, 'type': float, - 'help': "delta/beta material for Generalized Paganin"}, + 'help': "delta/beta material for Generalized Paganin & Farago"}, 'retrieve-phase-W': { 'default': 2e-4, 'type': float, diff --git a/src/tomocupy/processing/proc_functions.py b/src/tomocupy/processing/proc_functions.py index b8b5f8a..2e42f35 100644 --- a/src/tomocupy/processing/proc_functions.py +++ b/src/tomocupy/processing/proc_functions.py @@ -174,6 +174,10 @@ def proc_proj(self, data, st=None, end=None, res=None): if args.retrieve_phase_method == 'FF': data[:] = retrieve_phase.fresnel_filter( data, args.FFratio, args.FFdim, window=None, pad=args.FFpad, apply_log=args.FFlog) + if args.retrieve_phase_method == 'farago': + data[:] = retrieve_phase.farago_filter( + data, args.pixel_size*1e-4, args.propagation_distance/10, args.energy, + args.retrieve_phase_method, args.retrieve_phase_delta_beta) if args.rotate_proj_angle != 0: data[:] = self.rotate_proj( data, args.rotate_proj_angle, args.rotate_proj_order) diff --git a/src/tomocupy/processing/retrieve_phase.py b/src/tomocupy/processing/retrieve_phase.py index 0eff09c..4a72af2 100644 --- a/src/tomocupy/processing/retrieve_phase.py +++ b/src/tomocupy/processing/retrieve_phase.py @@ -107,6 +107,49 @@ def paganin_filter( return data +def farago_filter( + data, pixel_size=1e-4, dist=50, energy=20, db=1000, method='farago', pad=True): + """ + Perform single-step phase retrieval from phase-contrast measurements + :cite:`Farago:24`. + + Parameters + ---------- + tomo : ndarray + 3D tomographic data. + pixel_size : float, optional + Detector pixel size in cm. + dist : float, optional + Propagation distance of the wavefront in cm. + energy : float, optional + Energy of incident wave in keV. + method : string + phase retrieval method. Standard Paganin or Generalized Paganin. + db : float, optional + delta/beta for generalized Farago phase retrieval + pad : bool, optional + If True, extend the size of the projections by padding with zeros. + Returns + ------- + ndarray + Approximated 3D tomographic phase data. + """ + + # New dimensions and pad value after padding. + py, pz, val = _calc_pad(data, pixel_size, dist, energy, pad) + + # Compute the reciprocal grid. + dx, dy, dz = data.shape + if method == 'farago': + w2 = _reciprocal_grid(pixel_size, dy + 2 * py, dz + 2 * pz) + phase_filter = cp.fft.fftshift( + _farago_filter_factor(energy, dist, db, w2)) + + prj = cp.full((dy + 2 * py, dz + 2 * pz), val, dtype=data.dtype) + _retrieve_phase(data, phase_filter, py, pz, prj, pad) + + return data + def _retrieve_phase(data, phase_filter, px, py, prj, pad): dx, dy, dz = data.shape @@ -167,6 +210,8 @@ def _calc_pad(data, pixel_size, dist, energy, pad): def _paganin_filter_factor(energy, dist, alpha, w2): return 1 / (_wavelength(energy) * dist * w2 / (4 * PI) + alpha) +def _farago_filter_factor(energy, dist, db, w2): + return 1 / (cp.cos(PI*_wavelength(energy)*dist*w2) + db*cp.sin(PI*_wavelength(energy)*dist*w2)) def _paganin_filter_factorG(energy, dist, kf, pixel_size, db, W): """ From 41c53076ac913d90e11c0b2fb6dc120ba2845a97 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Fri, 25 Oct 2024 12:10:01 -0500 Subject: [PATCH 04/21] zarr Reader, various phase-retrievals --- list.dat | 379 ++++++++++++++++++++++ list_clean.dat | 81 +++++ out | 271 ++++++++++++++++ setup.py | 2 +- src/tomocupy/config.py | 2 +- src/tomocupy/dataio/reader.py | 12 +- src/tomocupy/dataio/writer.py | 93 +++++- src/tomocupy/processing/proc_functions.py | 3 +- src/tomocupy/processing/retrieve_phase.py | 7 +- src/tomocupy/utils.py | 176 ++++++++++ 10 files changed, 1014 insertions(+), 12 deletions(-) create mode 100644 list.dat create mode 100644 list_clean.dat create mode 100644 out diff --git a/list.dat b/list.dat new file mode 100644 index 0000000..3178980 --- /dev/null +++ b/list.dat @@ -0,0 +1,379 @@ +aiohttp==3.9.3 +aiosignal==1.3.1 +alabaster==0.7.16 +annotated-types==0.7.0 +anyio==4.3.0 +app-model==0.2.7 +appdirs==1.4.4 +argcomplete==1.12.0 +argon2-cffi==23.1.0 +argon2-cffi-bindings==21.2.0 +arrow==1.3.0 +asciitree==0.3.3 +asttokens==2.4.1 +async-lru==2.0.4 +async-timeout==4.0.3 +atomicwrites==1.4.1 +attrs==23.2.0 +Babel==2.14.0 +bcc==0.28.0 +bcrypt==4.1.3 +beautifulsoup4==4.12.3 +bleach==6.1.0 +blinker==1.8.2 +blivet==3.6.0 +boto3==1.34.63 +botocore==1.34.63 +Brlapi==0.8.2 +Brotli==1.1.0 +brotlipy==0.7.0 +build==1.2.1 +cachetools==5.3.3 +cachey==0.2.1 +certifi==2024.2.2 +cffi==1.16.0 +chardet==4.0.0 +charset-normalizer==3.3.2 +ciso8601==2.3.1 +click==8.1.7 +clikit==0.6.2 +cloud-files==4.23.0 +cloud-volume==8.31.0 +cloudpickle==3.0.0 +cmdstanpy==0.9.5 +cockpit @ file:///builddir/build/BUILD/cockpit-311.2/tmp/wheel/cockpit-311.2-py3-none-any.whl +colorama==0.4.6 +comm==0.2.2 +compressed-segmentation==2.2.2 +compresso==3.3.0 +configshell-fb==1.1.30 +contourpy==1.2.0 +convertdate==2.4.0 +crackle-codec==0.8.0 +crashtest==0.3.1 +crc32c==2.4 +cryptography==42.0.8 +cupshelpers==1.0 +cycler==0.12.1 +Cython==3.0.8 +dasbus==1.4 +dask==2021.12.0 +dbus-python==1.2.18 +debugpy==1.8.1 +decorator==4.4.2 +deflate==0.5.0 +defusedxml==0.7.1 +dill==0.3.8 +distributed==2021.12.0 +distro==1.5.0 +docstring_parser==0.16 +docutils==0.21.2 +DracoPy==1.4.0 +ephem==4.1.5 +ethtool==0.15 +exceptiongroup==1.2.0 +executing==2.0.1 +fabio==2024.4.0 +fasteners==0.19 +fastjsonschema==2.19.1 +fastremap==1.14.1 +fastrlock==0.8.2 +file-magic==0.4.0 +Flask==3.0.3 +flexcache==0.3 +flexparser==0.3.1 +fonttools==4.48.1 +fpzip==1.2.3 +fqdn==1.5.1 +freetype-py==2.4.0 +frozenlist==1.4.1 +fsspec==2022.5.0 +gevent==24.2.1 +globus-cli==3.29.0 +globus-sdk==3.41.0 +Glymur==0.13.2.post1 +google-api-core==2.17.1 +google-apitools==0.5.32 +google-auth==2.28.2 +google-cloud-core==2.4.1 +google-cloud-storage==2.15.0 +google-crc32c==1.5.0 +google-resumable-media==2.7.0 +googleapis-common-protos==1.63.0 +gpg==1.15.1 +greenlet==3.0.3 +h11==0.14.0 +h5py==3.10.0 +hdf5plugin==4.4.0 +HeapDict==1.0.1 +holidays==0.42 +hsluv==5.0.4 +httpcore==1.0.4 +httplib2==0.22.0 +httpstan==4.12.0 +httpx==0.27.0 +idna==2.10 +imageio==2.34.0 +imagej==0.3.2 +imagesize==1.4.1 +imglyb==2.1.0 +importlib-resources==6.1.1 +importlib_metadata==7.0.2 +in-n-out==0.2.1 +inflection==0.5.1 +iniconfig==2.0.0 +iniparse==0.4 +iotop==0.6 +ipython==8.18.1 +ipywidgets==8.1.2 +isc==2.0 +isoduration==20.11.0 +itsdangerous==2.2.0 +jedi==0.19.1 +jgo==1.0.5 +Jinja2==3.1.3 +jmespath==1.0.1 +JPype1==1.5.0 +json5==0.9.22 +jsonpointer==2.0 +jsonschema==4.21.1 +jsonschema-specifications==2023.12.1 +jupyter==1.0.0 +jupyter-console==6.6.3 +jupyter-events==0.9.1 +jupyter-lsp==2.2.4 +jupyter_client==8.6.1 +jupyter_core==5.7.2 +jupyter_server==2.13.0 +jupyter_server_terminals==0.5.3 +jupyterlab==4.1.5 +jupyterlab_pygments==0.3.0 +jupyterlab_server==2.25.4 +jupyterlab_widgets==3.0.10 +katello-host-tools==4.2.3 +kiwisolver==1.4.5 +kmod==0.1 +lab==8.1 +labeling==0.1.14 +langtable==0.0.54 +lazy_loader==0.4 +libcomps==0.1.18 +libtiff==0.4.2 +libvirt-python==10.0.0 +llvmlite==0.43.0 +locket==1.0.0 +LunarCalendar==0.0.9 +lxml==4.6.5 +magicgui==0.8.3 +Markdown==3.6 +markdown-it-py==3.0.0 +MarkupSafe==2.1.5 +marshmallow==3.20.2 +matplotlib==3.8.2 +matplotlib-inline==0.1.5 +mdurl==0.1.2 +mistune==3.0.2 +msgpack==1.0.8 +multidict==6.0.5 +multiprocess==0.70.16 +napari-console==0.0.9 +napari-plugin-engine==0.2.0 +napari-plugin-manager==0.1.0a2 +napari-svg==0.2.0 +nbclient==0.10.0 +nbconvert==7.16.2 +nbformat==5.10.3 +nest-asyncio==1.6.0 +netaddr==0.8.0 +networkx==3.2.1 +nftables==0.1 +ninja==1.11.1.1 +notebook==7.1.2 +notebook_shim==0.2.4 +npe2==0.7.6 +numba==0.60.0 +numcodecs==0.9.0 +numexpr==2.10.1 +numpy==1.23.0 +numpydoc==1.7.0 +nvmetcli==0.7 +oauth2client==4.1.3 +olefile==0.46 +opencv-python==4.10.0.84 +orjson==3.9.15 +outcome==1.3.0.post0 +overrides==7.7.0 +packaging==24.0 +pandas==2.0.3 +pandocfilters==1.5.1 +paramiko==2.12.0 +parso==0.8.3 +partd==1.4.2 +pastel==0.2.1 +pathos==0.3.2 +pcp==5.0 +pep517==0.13.1 +perf==0.1 +pexpect==4.8.0 +pid==2.2.3 +Pillow==10.0.1 +Pint==0.24.1 +platformdirs==4.2.0 +pluggy==1.4.0 +ply==3.11 +pooch==1.8.2 +posix-ipc==1.1.1 +pox==0.3.4 +ppft==1.7.6.8 +productmd==1.31 +prometheus_client==0.20.0 +prompt-toolkit==3.0.43 +protobuf==4.25.3 +psutil==5.8.0 +psycopg2==2.8.6 +psygnal==0.11.1 +ptyprocess==0.6.0 +pure-eval==0.2.2 +pvaserver @ file:///home/beams0/AMITTONE/Software/pvaServer +pwquality==1.4.4 +pyasn1==0.5.1 +pyasn1-modules==0.3.0 +pybind11==2.11.1 +pycairo==1.20.1 +pyconify==0.1.6 +pycparser==2.21 +pycups==2.0.1 +pycurl==7.43.0.6 +pydantic==2.8.2 +pydantic-compat==0.1.2 +pydantic_core==2.20.1 +pyenchant==3.2.0 +pyFFTW==0.13.1 +Pygments==2.17.2 +PyGObject==3.40.1 +pyimagej==1.5.0 +pyinotify==0.9.6 +PyJWT==2.8.0 +pykickstart==3.32.11 +pylev==1.4.0 +PyMeeus==0.5.12 +PyNaCl==1.5.0 +pyodbc===4.0.0-unsupported +PyOpenGL==3.1.7 +pyparsing==2.4.7 +pyparted==3.12.0 +pyproject_hooks==1.1.0 +pyprop @ file:///home/beams0/AMITTONE/Software/pyprop +PyQt5==5.15.10 +PyQt5-Qt5==5.15.14 +PyQt5_sip==12.15.0 +pysimdjson==5.0.2 +PySocks==1.7.1 +pyspng-seunglab==1.1.1 +pystan==3.9.0 +pytest==8.1.1 +python-dateutil==2.8.2 +python-json-logger==2.0.7 +python-jsonschema-objects==0.5.4 +python-linux-procfs==0.7.3 +python-meh==0.50 +pytz==2021.1 +pyudev==0.22.0 +PyWavelets==1.4.1 +pyxdg==0.27 +PyYAML==5.4.1 +pyzmq==25.1.2 +qtconsole==5.5.1 +QtPy==2.4.1 +referencing==0.33.0 +requests==2.25.1 +requests-file==1.5.1 +requests-ftp==0.3.1 +rfc3339-validator==0.1.4 +rfc3986-validator==0.1.1 +rich==13.7.1 +rpds-py==0.18.0 +rpm==4.16.1.3 +rsa==4.9 +rtslib-fb==2.1.76 +s3transfer==0.10.1 +scikit-build==0.18.1 +scikit-image==0.21.0 +scipy==1.10.1 +scp==0.15.0 +scyjava==1.10.0 +selenium==4.18.1 +selinux==3.6 +Send2Trash==1.8.2 +sepolicy==3.6 +setools==4.4.4 +setroubleshoot @ file:///builddir/build/BUILD/setroubleshoot-3.3.32/src +setuptools-git==1.2 +shellingham==1.5.4 +simplejpeg==1.7.2 +simplejson==3.19.2 +simpleline==1.8.3 +six==1.15.0 +sniffio==1.3.1 +snowballstemmer==2.2.0 +sockjs-tornado==1.0.7 +sortedcontainers==2.4.0 +sos==4.7.1 +soupsieve==2.5 +Sphinx==7.3.7 +sphinxcontrib-applehelp==1.0.8 +sphinxcontrib-devhelp==1.0.6 +sphinxcontrib-htmlhelp==2.0.5 +sphinxcontrib-jsmath==1.0.1 +sphinxcontrib-qthelp==1.0.7 +sphinxcontrib-serializinghtml==1.1.10 +stack-data==0.6.3 +subscription-manager==1.29.40 +superqt==0.6.7 +systemd-python==234 +tabulate==0.9.0 +targetcli-fb==2.1.57 +tblib==3.0.0 +tenacity==8.2.3 +terminado==0.18.1 +tifffile==2024.7.2 +tinycss2==1.2.1 +tk==0.1.0 +tomli==2.0.1 +tomli_w==1.0.0 +toolz==0.12.1 +tornado==6.4 +tqdm==4.66.2 +traitlets==5.14.2 +triangle==20230923 +trio==0.24.0 +trio-websocket==0.11.1 +tsx==0.1.16 +txt2tags==3.9 +typer==0.12.3 +types-python-dateutil==2.9.0.20240315 +typing_extensions==4.10.0 +tzdata==2024.1 +uri-template==1.3.0 +urllib3==1.26.5 +urwid==2.1.2 +vispy==0.14.3 +wcwidth==0.2.13 +webargs==8.4.0 +webcolors==1.13 +webencodings==0.5.1 +websocket-client==1.7.0 +Werkzeug==3.0.3 +widgetsnbextension==4.0.10 +wrapt==1.16.0 +wsproto==1.2.0 +xarray==2024.6.0 +yarl==1.9.4 +zarr==2.11.3 +zfpc==0.1.2 +zfpy==1.0.0 +zict==3.0.0 +zipp==3.17.0 +zope.event==5.0 +zope.interface==6.2 +zstandard==0.22.0 diff --git a/list_clean.dat b/list_clean.dat new file mode 100644 index 0000000..9319c35 --- /dev/null +++ b/list_clean.dat @@ -0,0 +1,81 @@ + +bcc==0.28.0 +blivet==3.6.0 +Brlapi==0.8.2 +chardet==4.0.0 +charset-normalizer==3.3.2 +cockpit @ file:///builddir/build/BUILD/cockpit-311.2/tmp/wheel/cockpit-311.2-py3-none-any.whl +configshell-fb==1.1.30 +cupshelpers==1.0 +dasbus==1.4 +dbus-python==1.2.18 +decorator==4.4.2 +distro==1.5.0 +ethtool==0.15 +file-magic==0.4.0 +gpg==1.15.1 +idna==2.10 +iniparse==0.4 +iotop==0.6 +isc==2.0 +jsonpointer==2.0 +katello-host-tools==4.2.3 +kmod==0.1 +langtable==0.0.54 +libcomps==0.1.18 +libvirt-python==10.0.0 +lxml==4.6.5 +matplotlib-inline==0.1.5 +netaddr==0.8.0 +nftables==0.1 +numpy==1.20.1 +nvmetcli==0.7 +olefile==0.46 +packaging==20.9 +pcp==5.0 +perf==0.1 +pexpect==4.8.0 +pid==2.2.3 +Pillow==10.0.1 +ply==3.11 +productmd==1.31 +psutil==5.8.0 +psycopg2==2.8.6 +ptyprocess==0.6.0 +pwquality==1.4.4 +pycairo==1.20.1 +pycups==2.0.1 +pycurl==7.43.0.6 +pyenchant==3.2.0 +PyGObject==3.40.1 +pyinotify==0.9.6 +pykickstart==3.32.11 +pyodbc===4.0.0-unsupported +pyparsing==2.4.7 +pyparted==3.12.0 +PySocks==1.7.1 +python-dateutil==2.8.1 +python-linux-procfs==0.7.3 +python-meh==0.50 +pytz==2021.1 +pyudev==0.22.0 +pyxdg==0.27 +PyYAML==5.4.1 +requests==2.25.1 +requests-file==1.5.1 +requests-ftp==0.3.1 +rpm==4.16.1.3 +rtslib-fb==2.1.76 +selinux==3.6 +sepolicy==3.6 +setools==4.4.4 +setroubleshoot @ file:///builddir/build/BUILD/setroubleshoot-3.3.32/src +simpleline==1.8.3 +six==1.15.0 +sos==4.7.1 +subscription-manager==1.29.40 +systemd-python==234 +targetcli-fb==2.1.57 +traitlets==5.1.1 +urllib3==1.26.5 +urwid==2.1.2 diff --git a/out b/out new file mode 100644 index 0000000..1522a2d --- /dev/null +++ b/out @@ -0,0 +1,271 @@ +usage: tomocupy recon [-h] [--binning {0,1,2,3}] + [--blocked-views BLOCKED_VIEWS] [--dark-file-name PATH] + [--file-name PATH] [--file-type {standard,double_fov}] + [--flat-file-name PATH] [--out-path-name PATH] + [--remove-stripe-method {none,fw,ti,vo-all}] + [--bright-ratio BRIGHT_RATIO] + [--center-search-step CENTER_SEARCH_STEP] + [--center-search-width CENTER_SEARCH_WIDTH] + [--clear-folder {True,False}] [--dezinger DEZINGER] + [--dezinger-threshold DEZINGER_THRESHOLD] + [--dtype {float32,float16}] [--end-column END_COLUMN] + [--end-proj END_PROJ] [--end-row END_ROW] + [--fbp-filter {none,ramp,shepp,hann,hamming,parzen,cosine,cosine2}] + [--find-center-end-row FIND_CENTER_END_ROW] + [--find-center-start-row FIND_CENTER_START_ROW] + [--flat-linear FLAT_LINEAR] + [--max-read-threads MAX_READ_THREADS] + [--max-write-threads MAX_WRITE_THREADS] + [--minus-log MINUS_LOG] + [--nproj-per-chunk NPROJ_PER_CHUNK] [--nsino NSINO] + [--nsino-per-chunk NSINO_PER_CHUNK] + [--pixel-size PIXEL_SIZE] + [--rotation-axis ROTATION_AXIS] + [--rotation-axis-auto {manual,auto}] + [--rotation-axis-method {sift,vo}] + [--rotation-axis-pairs ROTATION_AXIS_PAIRS] + [--rotation-axis-sift-threshold ROTATION_AXIS_SIFT_THRESHOLD] + [--save-format {tiff,h5,h5sino,h5nolinks}] + [--start-column START_COLUMN] [--start-proj START_PROJ] + [--start-row START_ROW] + [--fw-filter {haar,db5,sym5,sym16}] + [--fw-level FW_LEVEL] [--fw-pad] [--fw-sigma FW_SIGMA] + [--ti-beta TI_BETA] [--ti-mask TI_MASK] + [--vo-all-dim VO_ALL_DIM] + [--vo-all-la-size VO_ALL_LA_SIZE] + [--vo-all-sm-size VO_ALL_SM_SIZE] + [--vo-all-snr VO_ALL_SNR] [--lamino-angle LAMINO_ANGLE] + [--lamino-end-row LAMINO_END_ROW] + [--lamino-search-step LAMINO_SEARCH_STEP] + [--lamino-search-width LAMINO_SEARCH_WIDTH] + [--lamino-start-row LAMINO_START_ROW] + [--reconstruction-algorithm {fourierrec,lprec,linerec}] + [--reconstruction-type {full,try}] + [--b-storage-ring B_STORAGE_RING] + [--beam-hardening-method {none,standard}] + [--calculate-source {none,standard}] + [--e-storage-ring E_STORAGE_RING] + [--filter-1-auto FILTER_1_AUTO] + [--filter-1-density FILTER_1_DENSITY] + [--filter-1-material FILTER_1_MATERIAL] + [--filter-1-thickness FILTER_1_THICKNESS] + [--filter-2-auto FILTER_2_AUTO] + [--filter-2-density FILTER_2_DENSITY] + [--filter-2-material FILTER_2_MATERIAL] + [--filter-2-thickness FILTER_2_THICKNESS] + [--filter-3-auto FILTER_3_AUTO] + [--filter-3-density FILTER_3_DENSITY] + [--filter-3-material FILTER_3_MATERIAL] + [--filter-3-thickness FILTER_3_THICKNESS] + [--maximum-E MAXIMUM_E] + [--maximum-psi-urad MAXIMUM_PSI_URAD] + [--minimum-E MINIMUM_E] [--read-pixel-size] + [--read-scintillator] [--sample-density SAMPLE_DENSITY] + [--sample-material SAMPLE_MATERIAL] + [--scintillator-density SCINTILLATOR_DENSITY] + [--scintillator-material SCINTILLATOR_MATERIAL] + [--scintillator-thickness SCINTILLATOR_THICKNESS] + [--source-distance SOURCE_DISTANCE] [--step-E STEP_E] + [--config FILE] [--config-update] [--logs-home FILE] + [--verbose] + +optional arguments: + -h, --help show this help message and exit + --binning {0,1,2,3} Reconstruction binning factor as power(2, choice) + (default: 0) + --blocked-views BLOCKED_VIEWS + Angle range for blocked views [st,end]. Can be a list + of ranges(e.g. [[0,1.2],[3,3.14]]) (default: none) + --dark-file-name PATH + Name of the hdf file containing dark data (default: + None) + --file-name PATH Name of the last used hdf file or directory containing + multiple hdf files (default: .) + --file-type {standard,double_fov} + Input file type (default: standard) + --flat-file-name PATH + Name of the hdf file containing flat data (default: + None) + --out-path-name PATH Path for output files (default: None) + --remove-stripe-method {none,fw,ti,vo-all} + Remove stripe method: none, fourier-wavelet, titarenko + (default: none) + --bright-ratio BRIGHT_RATIO + exposure time for flat fields divided by the exposure + time of projections (default: 1) + --center-search-step CENTER_SEARCH_STEP + +/- center search step (pixel). (default: 0.5) + --center-search-width CENTER_SEARCH_WIDTH + +/- center search width (pixel). (default: 50.0) + --clear-folder {True,False} + Clear output folder before reconstruction (default: + False) + --dezinger DEZINGER Width of region for removing outliers (default: 0) + --dezinger-threshold DEZINGER_THRESHOLD + Threshold of grayscale above local median to be + considered a zinger pixel (default: 5000) + --dtype {float32,float16} + Data type used for reconstruction. Note float16 works + with power of 2 sizes. (default: float32) + --end-column END_COLUMN + End position in x (default: -1) + --end-proj END_PROJ End projection (default: -1) + --end-row END_ROW End slice (default: -1) + --fbp-filter {none,ramp,shepp,hann,hamming,parzen,cosine,cosine2} + Filter for FBP reconstruction (default: parzen) + --find-center-end-row FIND_CENTER_END_ROW + End row to find the rotation center (default: -1) + --find-center-start-row FIND_CENTER_START_ROW + Start row to find the rotation center (default: 0) + --flat-linear FLAT_LINEAR + Interpolate flat fields for each projections, assumes + the number of flat fields at the beginning of the scan + is as the same as a the end. (default: False) + --max-read-threads MAX_READ_THREADS + Max number of threads for reading by chunks (default: + 4) + --max-write-threads MAX_WRITE_THREADS + Max number of threads for writing by chunks (default: + 8) + --minus-log MINUS_LOG + Take -log or not (default: True) + --nproj-per-chunk NPROJ_PER_CHUNK + Number of sinograms per chunk. Use lower numbers with + computers with lower GPU memory. (default: 8) + --nsino NSINO Location of the sinogram used for slice reconstruction + and find axis (0 top, 1 bottom). Can be given as a + list, e.g. [0,0.9]. (default: 0.5) + --nsino-per-chunk NSINO_PER_CHUNK + Number of sinograms per chunk. Use larger numbers with + computers with larger memory. (default: 8) + --pixel-size PIXEL_SIZE + Pixel size [microns] (default: 0) + --rotation-axis ROTATION_AXIS + Location of rotation axis (default: -1.0) + --rotation-axis-auto {manual,auto} + How to get rotation axis auto calculate ('auto'), or + manually ('manual') (default: manual) + --rotation-axis-method {sift,vo} + Method for automatic rotation search. (default: sift) + --rotation-axis-pairs ROTATION_AXIS_PAIRS + Projection pairs to find rotation axis. Each second + projection in a pair will be flipped and used to find + shifts from the first element in a pair. The shifts + are used to calculate the center. Example [0,1499] for + a 180 deg scan, or [0,1499,749,2249] for 360, etc. + (default: [0,0]) + --rotation-axis-sift-threshold ROTATION_AXIS_SIFT_THRESHOLD + SIFT threshold for rotation search. (default: 0.5) + --save-format {tiff,h5,h5sino,h5nolinks} + Output format (default: tiff) + --start-column START_COLUMN + Start position in x (default: 0) + --start-proj START_PROJ + Start projection (default: 0) + --start-row START_ROW + Start slice (default: 0) + --fw-filter {haar,db5,sym5,sym16} + Fourier-Wavelet remove stripe filter (default: sym16) + --fw-level FW_LEVEL Fourier-Wavelet remove stripe level parameter + (default: 7) + --fw-pad When set, Fourier-Wavelet remove stripe extend the + size of the sinogram by padding with zeros (default: + True) + --fw-sigma FW_SIGMA Fourier-Wavelet remove stripe damping parameter + (default: 1) + --ti-beta TI_BETA Parameter for ring removal (0,1) (default: 0.022) + --ti-mask TI_MASK Mask size for ring removal (0,1) (default: 1) + --vo-all-dim VO_ALL_DIM + Dimension of the window. (default: 1) + --vo-all-la-size VO_ALL_LA_SIZE + Window size of the median filter to remove large + stripes. (default: 61) + --vo-all-sm-size VO_ALL_SM_SIZE + Window size of the median filter to remove small-to- + medium stripes. (default: 21) + --vo-all-snr VO_ALL_SNR + Ratio used to locate large stripes. Greater is less + sensitive. (default: 3) + --lamino-angle LAMINO_ANGLE + Pitch of the stage for laminography (default: 0) + --lamino-end-row LAMINO_END_ROW + End slice for lamino reconstruction (default: -1) + --lamino-search-step LAMINO_SEARCH_STEP + +/- center search step (pixel). (default: 0.25) + --lamino-search-width LAMINO_SEARCH_WIDTH + +/- center search width (pixel). (default: 5.0) + --lamino-start-row LAMINO_START_ROW + Start slice for lamino reconstruction (default: 0) + --reconstruction-algorithm {fourierrec,lprec,linerec} + Reconstruction algorithm (default: fourierrec) + --reconstruction-type {full,try} + Reconstruct full data set. (default: try) + --b-storage-ring B_STORAGE_RING + Magnetic field for BM source in T (default: 0.599) + --beam-hardening-method {none,standard} + Beam hardening method. (default: none) + --calculate-source {none,standard} + Use tabulated (none, default) or calculated source + (default: none) + --e-storage-ring E_STORAGE_RING + e-beam energy for BM source in GeV (default: 7.0) + --filter-1-auto FILTER_1_AUTO + If True, read filter 1 from HDF meta data (default: + False) + --filter-1-density FILTER_1_DENSITY + Filter 1 density in g/cm^3 (default: 1.0) + --filter-1-material FILTER_1_MATERIAL + Filter 1 material for beam hardening (default: none) + --filter-1-thickness FILTER_1_THICKNESS + Filter 1 thickness in microns (default: 0.0) + --filter-2-auto FILTER_2_AUTO + If True, read filter 2 from HDF meta data (default: + False) + --filter-2-density FILTER_2_DENSITY + Filter 2 density in g/cm^3 (default: 1.0) + --filter-2-material FILTER_2_MATERIAL + Filter 2 material for beam hardening (default: none) + --filter-2-thickness FILTER_2_THICKNESS + Filter 2 thickness in microns (default: 0.0) + --filter-3-auto FILTER_3_AUTO + If True, read filter 3 from HDF meta data (default: + False) + --filter-3-density FILTER_3_DENSITY + Filter 3 density in g/cm^3 (default: 1.0) + --filter-3-material FILTER_3_MATERIAL + Filter 3 material in microns (default: none) + --filter-3-thickness FILTER_3_THICKNESS + Filter 3 thickness for beam hardening (default: 0.0) + --maximum-E MAXIMUM_E + Maximum energy to model in eV (default: 200000) + --maximum-psi-urad MAXIMUM_PSI_URAD + Maximum vertical angle from centerline to model in + microradians (default: 40) + --minimum-E MINIMUM_E + Minimum energy to model in eV (default: 1000) + --read-pixel-size When set, read effective pixel size from the HDF file + (default: False) + --read-scintillator When set, read scintillator properties from the HDF + file (default: False) + --sample-density SAMPLE_DENSITY + Density of sample material in g/cm^3 (default: 1.0) + --sample-material SAMPLE_MATERIAL + Sample material for beam hardening (default: Fe) + --scintillator-density SCINTILLATOR_DENSITY + Density of scintillator in g/cm^3 (default: 6.0) + --scintillator-material SCINTILLATOR_MATERIAL + Scintillator material for beam hardening (default: + LuAG_Ce) + --scintillator-thickness SCINTILLATOR_THICKNESS + Scintillator thickness in microns (default: 100.0) + --source-distance SOURCE_DISTANCE + Distance from source to scintillator in m (default: + 36.0) + --step-E STEP_E Energy step in eV (default: 500) + --config FILE File name of configuration file (default: + /home/beams/AMITTONE/tomocupyon.conf) + --config-update When set, the content of the config file is updated + using the current params values (default: False) + --logs-home FILE Log file directory (default: + /home/beams/AMITTONE/logs) + --verbose Verbose output (default: False) diff --git a/setup.py b/setup.py index dded50b..72c2c97 100644 --- a/setup.py +++ b/setup.py @@ -10,7 +10,7 @@ packages=find_packages('src'), zip_safe=False, install_requires=[ - 'cupy', + 'cupy', 'opencv-python', 'h5py', 'numexpr', diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index 0d8a80f..4d6d1fc 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -420,7 +420,7 @@ def default_parameter(func, param): 'default': 'tiff', 'type': str, 'help': "Output format", - 'choices': ['tiff', 'h5', 'h5sino', 'h5nolinks']}, + 'choices': ['tiff', 'h5', 'h5sino', 'h5nolinks', 'zarr']}, 'clear-folder': { 'default': 'False', 'type': str, diff --git a/src/tomocupy/dataio/reader.py b/src/tomocupy/dataio/reader.py index ac75e66..cd3abca 100644 --- a/src/tomocupy/dataio/reader.py +++ b/src/tomocupy/dataio/reader.py @@ -84,7 +84,8 @@ def init_sizes(self): # read data sizes and projection angles with a reader sizes = self.read_sizes() - theta = self.read_theta() + print(sizes['nproji']) + theta = self.read_theta(sizes['nproji']) nproji = sizes['nproji'] nzi = sizes['nzi'] ni = sizes['ni'] @@ -306,12 +307,15 @@ def read_sizes(self): return sizes - def read_theta(self): + def read_theta(self, projections): """Read projection angles (in radians)""" with h5py.File(args.file_name) as file_in: - theta = file_in['/exchange/theta'][:].astype('float32')/180*np.pi - + if '/exchange/theta' in file_in: + theta = file_in['/exchange/theta'][:].astype('float32') / 180 * np.pi + else: + # If 'theta' doesn't exist, calculate it over the range [0, proj] + theta = np.linspace(0, np.pi, projections, dtype='float32') return theta def read_data_chunk_to_queue(self, data_queue, ids_proj, st_z, end_z, st_n, end_n, id_z, in_dtype): diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 4b8c76c..5858d20 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -47,6 +47,14 @@ import sys import tifffile +#Zarr writer +import shutil +import zarr +import json +from numcodecs import Blosc +from tomocupy.utils import downsampleZarr + + __author__ = "Viktor Nikitin" __copyright__ = "Copyright (c) 2022, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' @@ -202,9 +210,14 @@ def init_output_files(self): self.write_meta(rec_virtual) rec_virtual.close() - + if args.save_format == 'zarr': # Zarr format support + fnameout += '.zarr' + self.zarr_output_path = fnameout + log.info(f'Zarr dataset will be created at {fnameout}') + params.fnameout = fnameout log.info(f'Output: {fnameout}') + def write_meta(self, rec_virtual): @@ -250,9 +263,87 @@ def write_data_chunk(self, rec, st, end, k): with h5py.File(filename, "w") as fid: fid.create_dataset("/exchange/data", data=rec, chunks=(params.nproj, 1, params.n)) + elif args.save_format == 'zarr': # Use save_zarr to write Zarr data with pyramid levels + pixel_size = 1.0 # Replace with actual pixel size if available + chunks = (1, params.n, params.n) # Replace with appropriate chunk size + compression = 'zstd' # Zarr compression, adjust as necessary + save_zarr(volume=rec, output_path=self.zarr_output_path, chunks=chunks, compression=compression, pixel_size=pixel_size) + #log.info(f'Wrote chunk {k} to Zarr file with multiscale pyramid') def write_data_try(self, rec, cid, id_slice): """Write tiff reconstruction with a given name""" tifffile.imwrite( f'{params.fnameout}_slice{id_slice:04d}_center{cid:05.2f}.tiff', rec) + + +def save_zarr(volume, output_path, chunks, compression, pixel_size, mode='a', original_dtype=np.float32): + """ + Save a 3D volume to a Zarr store, creating a multiscale pyramid representation. + + Parameters: + - volume (numpy array): The 3D volume data to be saved. + - output_path (str): The path to the output Zarr store. + - chunks (tuple of ints): The chunk size for the Zarr array. + - compression (str): The compression algorithm to use (e.g., 'blosclz', 'lz4', etc.). + - pixel_size (float): The size of the pixels in micrometers. + - mode (str, optional): The mode to open the Zarr store ('w' for write, 'a' for append). Default is 'a'. + - original_dtype (numpy dtype, optional): The original data type of the images. Default is np.uint8. + + Returns: + - None + """ + store = zarr.DirectoryStore(output_path) + compressor = Blosc(cname=compression, clevel=5, shuffle=2) + + if mode == 'w': + if os.path.exists(output_path): + shutil.rmtree(output_path) + root_group = zarr.group(store=store) + else: + root_group = zarr.open(store=store, mode='a') + + # Assuming volume is a chunk of the data + pyramid_levels = downsampleZarr(volume) # Assuming downsample is defined elsewhere for multiscale + + datasets = [] + for level, data in enumerate(pyramid_levels): + data = data.astype(original_dtype) + + dataset_name = f"{level}" + if dataset_name in root_group: + z = root_group[dataset_name] + z.append(data, axis=0) + else: + z = root_group.create_dataset(name=dataset_name, shape=data.shape, chunks=chunks, dtype=data.dtype, compressor=compressor) + z[:] = data + + scale_factor = 2 ** level + datasets.append({ + "path": dataset_name, + "coordinateTransformations": [{"type": "scale", "scale": [pixel_size * scale_factor, pixel_size * scale_factor, pixel_size * scale_factor]}] + }) + + if mode == 'w': + multiscales = [{ + "version": "0.4", + "name": "example", + "axes": [ + {"name": "z", "type": "space", "unit": "micrometer"}, + {"name": "y", "type": "space", "unit": "micrometer"}, + {"name": "x", "type": "space", "unit": "micrometer"} + ], + "datasets": datasets, + "type": "gaussian", + "metadata": { + "method": "skimage.transform.resize", + "version": "0.16.1", + "args": "[true]", + "kwargs": {"anti_aliasing": True, "preserve_range": True} + } + }] + + root_group.attrs.update({"multiscales": multiscales}) + with open(os.path.join(output_path, 'multiscales.json'), 'w') as f: + json.dump({"multiscales": multiscales}, f, indent=2) + #info(f"Metadata saved to {output_path}") diff --git a/src/tomocupy/processing/proc_functions.py b/src/tomocupy/processing/proc_functions.py index 2e42f35..5a339a0 100644 --- a/src/tomocupy/processing/proc_functions.py +++ b/src/tomocupy/processing/proc_functions.py @@ -177,7 +177,8 @@ def proc_proj(self, data, st=None, end=None, res=None): if args.retrieve_phase_method == 'farago': data[:] = retrieve_phase.farago_filter( data, args.pixel_size*1e-4, args.propagation_distance/10, args.energy, - args.retrieve_phase_method, args.retrieve_phase_delta_beta) + args.retrieve_phase_delta_beta, args.retrieve_phase_method) + if args.rotate_proj_angle != 0: data[:] = self.rotate_proj( data, args.rotate_proj_angle, args.rotate_proj_order) diff --git a/src/tomocupy/processing/retrieve_phase.py b/src/tomocupy/processing/retrieve_phase.py index 4a72af2..d52911b 100644 --- a/src/tomocupy/processing/retrieve_phase.py +++ b/src/tomocupy/processing/retrieve_phase.py @@ -140,10 +140,9 @@ def farago_filter( # Compute the reciprocal grid. dx, dy, dz = data.shape - if method == 'farago': - w2 = _reciprocal_grid(pixel_size, dy + 2 * py, dz + 2 * pz) - phase_filter = cp.fft.fftshift( - _farago_filter_factor(energy, dist, db, w2)) + w2 = _reciprocal_grid(pixel_size, dy + 2 * py, dz + 2 * pz) + phase_filter = cp.fft.fftshift( + _farago_filter_factor(energy, dist, db, w2)) prj = cp.full((dy + 2 * py, dz + 2 * pz), val, dtype=data.dtype) _retrieve_phase(data, phase_filter, py, pz, prj, pad) diff --git a/src/tomocupy/utils.py b/src/tomocupy/utils.py index e7c2afd..eadf96d 100644 --- a/src/tomocupy/utils.py +++ b/src/tomocupy/utils.py @@ -47,6 +47,15 @@ import time import numexpr as ne import sys +import os +import tifffile as tiff +from skimage.transform import resize +import time +from functools import wraps + + + + from tomocupy import logging log = logging.getLogger(__name__) @@ -265,3 +274,170 @@ def param_from_dxchange(hdf_file, data_path, attr=None, scalar=True, char_array= return None except KeyError: return None + + +def progress_bar(func): + @wraps(func) + def wrapper(*args, **kwargs): + total = kwargs.pop('total', 100) # Total number of iterations (default is 100) + bar_length = 40 # Length of the progress bar + + for i in range(total): + result = func(*args, **kwargs) # Call the original function + percent_complete = (i + 1) / total + bar = '#' * int(percent_complete * bar_length) + '-' * (bar_length - int(percent_complete * bar_length)) + sys.stdout.write(f'\rProgress: [{bar}] {percent_complete:.2%}') + sys.stdout.flush() + time.sleep(0.01) # Simulate work being done + + print("\nCompleted!") + return result + + return wrapper + + +@progress_bar +def calculate_global_min_max(input_dir, bin_factor=4): + """ + Calculate the global minimum and maximum pixel values of all TIFF images in a directory after downsampling (binning). + + Parameters: + - input_dir (str): Path to the directory containing TIFF images. + - bin_factor (int, optional): Factor by which to downsample the images before calculating min and max values. Default is 4. + + Returns: + - global_min (float): The minimum pixel value across all images. + - global_max (float): The maximum pixel value across all images. + + Raises: + - ValueError: If no TIFF files are found in the directory. + """ + file_list = sorted([os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.endswith(('.tiff','.tif'))]) + if not file_list: + raise ValueError(f"No TIFF files found in the directory: {input_dir}") + + global_min = np.inf + global_max = -np.inf + + for file in file_list: + image = tiff.imread(file) + binned_image = resize(image, + (image.shape[0] // bin_factor, image.shape[1] // bin_factor), + order=1, preserve_range=True, anti_aliasing=False).astype(image.dtype) + global_min = min(global_min, binned_image.min()) + global_max = max(global_max, binned_image.max()) + + #info(f"Global min and max found: {global_min}, dtype: {global_max}") + return global_min, global_max + + + +def minmaxHisto(input_dir, thr=1e-5, num_bins=1000): + # Read the image files + file_list = sorted([os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.endswith(('.tiff','.tif'))]) + if not file_list: + raise ValueError(f"No TIFF files found in the directory: {input_dir}") + + # Choose the middle image from the list + middle_image_path = file_list[len(file_list) // 2] + + # Read the image + image = tiff.imread(middle_image_path) + if image is None: + raise ValueError(f"Image not found at {middle_image_path}") + + # Calculate the histogram + hist, bin_edges = np.histogram(image, bins=num_bins) + + # Find the start and end indices based on a threshold + threshold = np.max(hist) * thr + stend = np.where(hist > threshold) + if len(stend[0]) == 0: + raise ValueError("No significant histogram bins found.") + + st = stend[0][0] + end = stend[0][-1] + + # Determine min and max values + mmin = bin_edges[st] + mmax = bin_edges[end + 1] + + # Ensure min and max are not too close + if np.isclose(mmin, mmax): + raise ValueError("The minimum and maximum values are too close. Adjust the threshold or bin count.") + + return mmin, mmax + + + +def load_tiff_chunked(input_dir, dtype, chunk_size, start_index=0, global_min=None, global_max=None): + """ + Load TIFF images from a directory in chunks and convert them to a specified data type, optionally normalizing the values. + + Parameters: + - input_dir (str): Path to the directory containing TIFF images. + - dtype (numpy dtype): Target data type for the output array. + - chunk_size (int): Number of images to load in each chunk. + - start_index (int, optional): Starting index for loading images. Default is 0. + - global_min (float, optional): Minimum pixel value for normalization. If None, no normalization is performed. Default is None. + - global_max (float, optional): Maximum pixel value for normalization. If None, no normalization is performed. Default is None. + + Returns: + - zarr_chunk (numpy array): A chunk of loaded images converted to the specified data type. + - end_index (int): The end index for the current chunk of images. + + Raises: + - ValueError: If no TIFF files are found in the directory. + """ + file_list = sorted([os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.endswith(('.tiff','.tif'))]) + if not file_list: + raise ValueError(f"No TIFF files found in the directory: {input_dir}") + + end_index = min(start_index + chunk_size, len(file_list)) + chunk_files = file_list[start_index:end_index] + if not chunk_files: + return None, end_index + + sample_image = tiff.imread(chunk_files[0]) + chunk_shape = (len(chunk_files),) + sample_image.shape + + zarr_chunk = np.zeros(chunk_shape, dtype=dtype) + + for i, file in enumerate(chunk_files): + image = tiff.imread(file) + if global_min is not None and global_max is not None: + image = (image - global_min) / (global_max - global_min) # Normalize to [0, 1] + if dtype == np.uint16 or dtype == np.int16: + image = image * (2**15 - 1) # Scale to int16 range + elif dtype == np.uint8 or dtype == np.int8: + image = image * (2**8 - 1) # Scale to int8 range + + zarr_chunk[i] = image.astype(dtype) + + #info(f"Loaded TIFF chunk with shape: {zarr_chunk.shape}, dtype: {zarr_chunk.dtype}") + return zarr_chunk, end_index + +def downsampleZarr(data, scale_factor=2, max_levels=6): + """ + Create a multi-level downsampled version of the input data. + + Parameters: + - data (numpy array): The input image data to be downsampled. + - scale_factor (int, optional): Factor by which to downsample the data at each level. Default is 2. + - max_levels (int, optional): Maximum number of downsampled levels to generate. Default is 6. + + Returns: + - levels (list of numpy arrays): A list containing the original data and each downsampled level. + + Logs: + - Information about the shape and data type of each downsampled level. + """ + current_level = data + levels = [current_level] + for _ in range(max_levels): + new_shape = tuple(max(1, dim // 2) for dim in current_level.shape) + if min(new_shape) <= 1: + break + current_level = resize(current_level, new_shape, order=0, preserve_range=True, anti_aliasing=True) + levels.append(current_level) + return levels From 9a8308068523561edcc343ec21fac17d6c14361d Mon Sep 17 00:00:00 2001 From: mittoalb Date: Fri, 25 Oct 2024 15:25:42 -0500 Subject: [PATCH 05/21] config fix --- src/tomocupy/config.py | 5 +++++ src/tomocupy/dataio/writer.py | 7 ++----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index 4d6d1fc..b77f08d 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -421,6 +421,11 @@ def default_parameter(func, param): 'type': str, 'help': "Output format", 'choices': ['tiff', 'h5', 'h5sino', 'h5nolinks', 'zarr']}, + 'zarr-compression': { + 'default': 'blosclz', + 'type': str, + 'help': "ZARR compression format", + 'choices': ['blosclz', 'lz4', 'zstd']}, 'clear-folder': { 'default': 'False', 'type': str, diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 5858d20..b2d2c2e 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -263,12 +263,9 @@ def write_data_chunk(self, rec, st, end, k): with h5py.File(filename, "w") as fid: fid.create_dataset("/exchange/data", data=rec, chunks=(params.nproj, 1, params.n)) - elif args.save_format == 'zarr': # Use save_zarr to write Zarr data with pyramid levels - pixel_size = 1.0 # Replace with actual pixel size if available + elif args.save_format == 'zarr': # chunks = (1, params.n, params.n) # Replace with appropriate chunk size - compression = 'zstd' # Zarr compression, adjust as necessary - save_zarr(volume=rec, output_path=self.zarr_output_path, chunks=chunks, compression=compression, pixel_size=pixel_size) - #log.info(f'Wrote chunk {k} to Zarr file with multiscale pyramid') + save_zarr(volume=rec, output_path=self.zarr_output_path, chunks=chunks, compression=args.zarr_compression, pixel_size=args.pixel_size) def write_data_try(self, rec, cid, id_slice): """Write tiff reconstruction with a given name""" From 355899a6e131ddaf66f25c4c80ecb6bfc14ea2d4 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Tue, 29 Oct 2024 12:23:36 -0500 Subject: [PATCH 06/21] code cleanup --- src/tomocupy/config.py | 49 +++++++------- src/tomocupy/utils.py | 143 +---------------------------------------- 2 files changed, 27 insertions(+), 165 deletions(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index b77f08d..a09c992 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -224,7 +224,7 @@ def default_parameter(func, param): 'default': 'none', 'type': str, 'help': "Phase retrieval correction method", - 'choices': ['none', 'paganin', 'Gpaganin', 'FF', 'farago']}, + 'choices': ['none', 'paganin', 'Gpaganin', 'FourierFilter', 'farago']}, 'energy': { 'default': 0, 'type': float, @@ -411,26 +411,6 @@ def default_parameter(func, param): 'type': int, 'default': -1, 'help': "End row to find the rotation center"}, - 'dtype': { - 'default': 'float32', - 'type': str, - 'choices': ['float32', 'float16'], - 'help': "Data type used for reconstruction. Note float16 works with power of 2 sizes.", }, - 'save-format': { - 'default': 'tiff', - 'type': str, - 'help': "Output format", - 'choices': ['tiff', 'h5', 'h5sino', 'h5nolinks', 'zarr']}, - 'zarr-compression': { - 'default': 'blosclz', - 'type': str, - 'help': "ZARR compression format", - 'choices': ['blosclz', 'lz4', 'zstd']}, - 'clear-folder': { - 'default': 'False', - 'type': str, - 'help': "Clear output folder before reconstruction", - 'choices': ['True', 'False']}, 'fbp-filter': { 'default': 'parzen', 'type': str, @@ -467,6 +447,29 @@ def default_parameter(func, param): 'type': float, 'help': "Pixel size [microns]"}, } +SECTIONS['output'] ={ + 'dtype': { + 'default': 'float32', + 'type': str, + 'choices': ['float32', 'float16'], + 'help': "Data type used for reconstruction. Note float16 works with power of 2 sizes.", }, + 'save-format': { + 'default': 'tiff', + 'type': str, + 'help': "Output format", + 'choices': ['tiff', 'h5', 'h5sino', 'h5nolinks', 'zarr']}, + 'zarr-compression': { + 'default': 'blosclz', + 'type': str, + 'help': "ZARR compression format", + 'choices': ['blosclz', 'lz4', 'zstd']}, + 'clear-folder': { + 'default': 'False', + 'type': str, + 'help': "Clear output folder before reconstruction", + 'choices': ['True', 'False']} +} + SECTIONS['beam-hardening'] = { 'beam-hardening-method': { @@ -585,9 +588,9 @@ def default_parameter(func, param): RECON_PARAMS = ('file-reading', 'remove-stripe', - 'reconstruction', 'fw', 'ti', 'vo-all', 'lamino', 'reconstruction-types', 'beam-hardening') + 'reconstruction', 'fw', 'ti', 'vo-all', 'lamino', 'reconstruction-types', 'beam-hardening', 'output') RECON_STEPS_PARAMS = ('file-reading', 'remove-stripe', 'reconstruction', - 'retrieve-phase', 'fw', 'ti', 'vo-all', 'lamino', 'reconstruction-steps-types', 'rotate-proj', 'beam-hardening') + 'retrieve-phase', 'fw', 'ti', 'vo-all', 'lamino', 'reconstruction-steps-types', 'rotate-proj', 'beam-hardening', 'output') NICE_NAMES = ('General', 'File reading', 'Remove stripe', 'Remove stripe FW', 'Remove stripe Titarenko', 'Remove stripe Vo' 'Retrieve phase', 'Reconstruction') diff --git a/src/tomocupy/utils.py b/src/tomocupy/utils.py index eadf96d..9bae74c 100644 --- a/src/tomocupy/utils.py +++ b/src/tomocupy/utils.py @@ -275,148 +275,7 @@ def param_from_dxchange(hdf_file, data_path, attr=None, scalar=True, char_array= except KeyError: return None - -def progress_bar(func): - @wraps(func) - def wrapper(*args, **kwargs): - total = kwargs.pop('total', 100) # Total number of iterations (default is 100) - bar_length = 40 # Length of the progress bar - - for i in range(total): - result = func(*args, **kwargs) # Call the original function - percent_complete = (i + 1) / total - bar = '#' * int(percent_complete * bar_length) + '-' * (bar_length - int(percent_complete * bar_length)) - sys.stdout.write(f'\rProgress: [{bar}] {percent_complete:.2%}') - sys.stdout.flush() - time.sleep(0.01) # Simulate work being done - - print("\nCompleted!") - return result - - return wrapper - - -@progress_bar -def calculate_global_min_max(input_dir, bin_factor=4): - """ - Calculate the global minimum and maximum pixel values of all TIFF images in a directory after downsampling (binning). - - Parameters: - - input_dir (str): Path to the directory containing TIFF images. - - bin_factor (int, optional): Factor by which to downsample the images before calculating min and max values. Default is 4. - - Returns: - - global_min (float): The minimum pixel value across all images. - - global_max (float): The maximum pixel value across all images. - - Raises: - - ValueError: If no TIFF files are found in the directory. - """ - file_list = sorted([os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.endswith(('.tiff','.tif'))]) - if not file_list: - raise ValueError(f"No TIFF files found in the directory: {input_dir}") - - global_min = np.inf - global_max = -np.inf - - for file in file_list: - image = tiff.imread(file) - binned_image = resize(image, - (image.shape[0] // bin_factor, image.shape[1] // bin_factor), - order=1, preserve_range=True, anti_aliasing=False).astype(image.dtype) - global_min = min(global_min, binned_image.min()) - global_max = max(global_max, binned_image.max()) - - #info(f"Global min and max found: {global_min}, dtype: {global_max}") - return global_min, global_max - - - -def minmaxHisto(input_dir, thr=1e-5, num_bins=1000): - # Read the image files - file_list = sorted([os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.endswith(('.tiff','.tif'))]) - if not file_list: - raise ValueError(f"No TIFF files found in the directory: {input_dir}") - - # Choose the middle image from the list - middle_image_path = file_list[len(file_list) // 2] - - # Read the image - image = tiff.imread(middle_image_path) - if image is None: - raise ValueError(f"Image not found at {middle_image_path}") - - # Calculate the histogram - hist, bin_edges = np.histogram(image, bins=num_bins) - - # Find the start and end indices based on a threshold - threshold = np.max(hist) * thr - stend = np.where(hist > threshold) - if len(stend[0]) == 0: - raise ValueError("No significant histogram bins found.") - - st = stend[0][0] - end = stend[0][-1] - - # Determine min and max values - mmin = bin_edges[st] - mmax = bin_edges[end + 1] - - # Ensure min and max are not too close - if np.isclose(mmin, mmax): - raise ValueError("The minimum and maximum values are too close. Adjust the threshold or bin count.") - - return mmin, mmax - - - -def load_tiff_chunked(input_dir, dtype, chunk_size, start_index=0, global_min=None, global_max=None): - """ - Load TIFF images from a directory in chunks and convert them to a specified data type, optionally normalizing the values. - - Parameters: - - input_dir (str): Path to the directory containing TIFF images. - - dtype (numpy dtype): Target data type for the output array. - - chunk_size (int): Number of images to load in each chunk. - - start_index (int, optional): Starting index for loading images. Default is 0. - - global_min (float, optional): Minimum pixel value for normalization. If None, no normalization is performed. Default is None. - - global_max (float, optional): Maximum pixel value for normalization. If None, no normalization is performed. Default is None. - - Returns: - - zarr_chunk (numpy array): A chunk of loaded images converted to the specified data type. - - end_index (int): The end index for the current chunk of images. - - Raises: - - ValueError: If no TIFF files are found in the directory. - """ - file_list = sorted([os.path.join(input_dir, f) for f in os.listdir(input_dir) if f.endswith(('.tiff','.tif'))]) - if not file_list: - raise ValueError(f"No TIFF files found in the directory: {input_dir}") - - end_index = min(start_index + chunk_size, len(file_list)) - chunk_files = file_list[start_index:end_index] - if not chunk_files: - return None, end_index - - sample_image = tiff.imread(chunk_files[0]) - chunk_shape = (len(chunk_files),) + sample_image.shape - - zarr_chunk = np.zeros(chunk_shape, dtype=dtype) - - for i, file in enumerate(chunk_files): - image = tiff.imread(file) - if global_min is not None and global_max is not None: - image = (image - global_min) / (global_max - global_min) # Normalize to [0, 1] - if dtype == np.uint16 or dtype == np.int16: - image = image * (2**15 - 1) # Scale to int16 range - elif dtype == np.uint8 or dtype == np.int8: - image = image * (2**8 - 1) # Scale to int8 range - - zarr_chunk[i] = image.astype(dtype) - - #info(f"Loaded TIFF chunk with shape: {zarr_chunk.shape}, dtype: {zarr_chunk.dtype}") - return zarr_chunk, end_index - + def downsampleZarr(data, scale_factor=2, max_levels=6): """ Create a multi-level downsampled version of the input data. From 291dfbbd1d62db73dd853908ca6d6aa26d31098a Mon Sep 17 00:00:00 2001 From: mittoalb Date: Wed, 30 Oct 2024 13:43:21 -0500 Subject: [PATCH 07/21] code cleanup --- list.dat | 379 ------------------------------------------------- list_clean.dat | 81 ----------- out | 271 ----------------------------------- 3 files changed, 731 deletions(-) delete mode 100644 list.dat delete mode 100644 list_clean.dat delete mode 100644 out diff --git a/list.dat b/list.dat deleted file mode 100644 index 3178980..0000000 --- a/list.dat +++ /dev/null @@ -1,379 +0,0 @@ -aiohttp==3.9.3 -aiosignal==1.3.1 -alabaster==0.7.16 -annotated-types==0.7.0 -anyio==4.3.0 -app-model==0.2.7 -appdirs==1.4.4 -argcomplete==1.12.0 -argon2-cffi==23.1.0 -argon2-cffi-bindings==21.2.0 -arrow==1.3.0 -asciitree==0.3.3 -asttokens==2.4.1 -async-lru==2.0.4 -async-timeout==4.0.3 -atomicwrites==1.4.1 -attrs==23.2.0 -Babel==2.14.0 -bcc==0.28.0 -bcrypt==4.1.3 -beautifulsoup4==4.12.3 -bleach==6.1.0 -blinker==1.8.2 -blivet==3.6.0 -boto3==1.34.63 -botocore==1.34.63 -Brlapi==0.8.2 -Brotli==1.1.0 -brotlipy==0.7.0 -build==1.2.1 -cachetools==5.3.3 -cachey==0.2.1 -certifi==2024.2.2 -cffi==1.16.0 -chardet==4.0.0 -charset-normalizer==3.3.2 -ciso8601==2.3.1 -click==8.1.7 -clikit==0.6.2 -cloud-files==4.23.0 -cloud-volume==8.31.0 -cloudpickle==3.0.0 -cmdstanpy==0.9.5 -cockpit @ file:///builddir/build/BUILD/cockpit-311.2/tmp/wheel/cockpit-311.2-py3-none-any.whl -colorama==0.4.6 -comm==0.2.2 -compressed-segmentation==2.2.2 -compresso==3.3.0 -configshell-fb==1.1.30 -contourpy==1.2.0 -convertdate==2.4.0 -crackle-codec==0.8.0 -crashtest==0.3.1 -crc32c==2.4 -cryptography==42.0.8 -cupshelpers==1.0 -cycler==0.12.1 -Cython==3.0.8 -dasbus==1.4 -dask==2021.12.0 -dbus-python==1.2.18 -debugpy==1.8.1 -decorator==4.4.2 -deflate==0.5.0 -defusedxml==0.7.1 -dill==0.3.8 -distributed==2021.12.0 -distro==1.5.0 -docstring_parser==0.16 -docutils==0.21.2 -DracoPy==1.4.0 -ephem==4.1.5 -ethtool==0.15 -exceptiongroup==1.2.0 -executing==2.0.1 -fabio==2024.4.0 -fasteners==0.19 -fastjsonschema==2.19.1 -fastremap==1.14.1 -fastrlock==0.8.2 -file-magic==0.4.0 -Flask==3.0.3 -flexcache==0.3 -flexparser==0.3.1 -fonttools==4.48.1 -fpzip==1.2.3 -fqdn==1.5.1 -freetype-py==2.4.0 -frozenlist==1.4.1 -fsspec==2022.5.0 -gevent==24.2.1 -globus-cli==3.29.0 -globus-sdk==3.41.0 -Glymur==0.13.2.post1 -google-api-core==2.17.1 -google-apitools==0.5.32 -google-auth==2.28.2 -google-cloud-core==2.4.1 -google-cloud-storage==2.15.0 -google-crc32c==1.5.0 -google-resumable-media==2.7.0 -googleapis-common-protos==1.63.0 -gpg==1.15.1 -greenlet==3.0.3 -h11==0.14.0 -h5py==3.10.0 -hdf5plugin==4.4.0 -HeapDict==1.0.1 -holidays==0.42 -hsluv==5.0.4 -httpcore==1.0.4 -httplib2==0.22.0 -httpstan==4.12.0 -httpx==0.27.0 -idna==2.10 -imageio==2.34.0 -imagej==0.3.2 -imagesize==1.4.1 -imglyb==2.1.0 -importlib-resources==6.1.1 -importlib_metadata==7.0.2 -in-n-out==0.2.1 -inflection==0.5.1 -iniconfig==2.0.0 -iniparse==0.4 -iotop==0.6 -ipython==8.18.1 -ipywidgets==8.1.2 -isc==2.0 -isoduration==20.11.0 -itsdangerous==2.2.0 -jedi==0.19.1 -jgo==1.0.5 -Jinja2==3.1.3 -jmespath==1.0.1 -JPype1==1.5.0 -json5==0.9.22 -jsonpointer==2.0 -jsonschema==4.21.1 -jsonschema-specifications==2023.12.1 -jupyter==1.0.0 -jupyter-console==6.6.3 -jupyter-events==0.9.1 -jupyter-lsp==2.2.4 -jupyter_client==8.6.1 -jupyter_core==5.7.2 -jupyter_server==2.13.0 -jupyter_server_terminals==0.5.3 -jupyterlab==4.1.5 -jupyterlab_pygments==0.3.0 -jupyterlab_server==2.25.4 -jupyterlab_widgets==3.0.10 -katello-host-tools==4.2.3 -kiwisolver==1.4.5 -kmod==0.1 -lab==8.1 -labeling==0.1.14 -langtable==0.0.54 -lazy_loader==0.4 -libcomps==0.1.18 -libtiff==0.4.2 -libvirt-python==10.0.0 -llvmlite==0.43.0 -locket==1.0.0 -LunarCalendar==0.0.9 -lxml==4.6.5 -magicgui==0.8.3 -Markdown==3.6 -markdown-it-py==3.0.0 -MarkupSafe==2.1.5 -marshmallow==3.20.2 -matplotlib==3.8.2 -matplotlib-inline==0.1.5 -mdurl==0.1.2 -mistune==3.0.2 -msgpack==1.0.8 -multidict==6.0.5 -multiprocess==0.70.16 -napari-console==0.0.9 -napari-plugin-engine==0.2.0 -napari-plugin-manager==0.1.0a2 -napari-svg==0.2.0 -nbclient==0.10.0 -nbconvert==7.16.2 -nbformat==5.10.3 -nest-asyncio==1.6.0 -netaddr==0.8.0 -networkx==3.2.1 -nftables==0.1 -ninja==1.11.1.1 -notebook==7.1.2 -notebook_shim==0.2.4 -npe2==0.7.6 -numba==0.60.0 -numcodecs==0.9.0 -numexpr==2.10.1 -numpy==1.23.0 -numpydoc==1.7.0 -nvmetcli==0.7 -oauth2client==4.1.3 -olefile==0.46 -opencv-python==4.10.0.84 -orjson==3.9.15 -outcome==1.3.0.post0 -overrides==7.7.0 -packaging==24.0 -pandas==2.0.3 -pandocfilters==1.5.1 -paramiko==2.12.0 -parso==0.8.3 -partd==1.4.2 -pastel==0.2.1 -pathos==0.3.2 -pcp==5.0 -pep517==0.13.1 -perf==0.1 -pexpect==4.8.0 -pid==2.2.3 -Pillow==10.0.1 -Pint==0.24.1 -platformdirs==4.2.0 -pluggy==1.4.0 -ply==3.11 -pooch==1.8.2 -posix-ipc==1.1.1 -pox==0.3.4 -ppft==1.7.6.8 -productmd==1.31 -prometheus_client==0.20.0 -prompt-toolkit==3.0.43 -protobuf==4.25.3 -psutil==5.8.0 -psycopg2==2.8.6 -psygnal==0.11.1 -ptyprocess==0.6.0 -pure-eval==0.2.2 -pvaserver @ file:///home/beams0/AMITTONE/Software/pvaServer -pwquality==1.4.4 -pyasn1==0.5.1 -pyasn1-modules==0.3.0 -pybind11==2.11.1 -pycairo==1.20.1 -pyconify==0.1.6 -pycparser==2.21 -pycups==2.0.1 -pycurl==7.43.0.6 -pydantic==2.8.2 -pydantic-compat==0.1.2 -pydantic_core==2.20.1 -pyenchant==3.2.0 -pyFFTW==0.13.1 -Pygments==2.17.2 -PyGObject==3.40.1 -pyimagej==1.5.0 -pyinotify==0.9.6 -PyJWT==2.8.0 -pykickstart==3.32.11 -pylev==1.4.0 -PyMeeus==0.5.12 -PyNaCl==1.5.0 -pyodbc===4.0.0-unsupported -PyOpenGL==3.1.7 -pyparsing==2.4.7 -pyparted==3.12.0 -pyproject_hooks==1.1.0 -pyprop @ file:///home/beams0/AMITTONE/Software/pyprop -PyQt5==5.15.10 -PyQt5-Qt5==5.15.14 -PyQt5_sip==12.15.0 -pysimdjson==5.0.2 -PySocks==1.7.1 -pyspng-seunglab==1.1.1 -pystan==3.9.0 -pytest==8.1.1 -python-dateutil==2.8.2 -python-json-logger==2.0.7 -python-jsonschema-objects==0.5.4 -python-linux-procfs==0.7.3 -python-meh==0.50 -pytz==2021.1 -pyudev==0.22.0 -PyWavelets==1.4.1 -pyxdg==0.27 -PyYAML==5.4.1 -pyzmq==25.1.2 -qtconsole==5.5.1 -QtPy==2.4.1 -referencing==0.33.0 -requests==2.25.1 -requests-file==1.5.1 -requests-ftp==0.3.1 -rfc3339-validator==0.1.4 -rfc3986-validator==0.1.1 -rich==13.7.1 -rpds-py==0.18.0 -rpm==4.16.1.3 -rsa==4.9 -rtslib-fb==2.1.76 -s3transfer==0.10.1 -scikit-build==0.18.1 -scikit-image==0.21.0 -scipy==1.10.1 -scp==0.15.0 -scyjava==1.10.0 -selenium==4.18.1 -selinux==3.6 -Send2Trash==1.8.2 -sepolicy==3.6 -setools==4.4.4 -setroubleshoot @ file:///builddir/build/BUILD/setroubleshoot-3.3.32/src -setuptools-git==1.2 -shellingham==1.5.4 -simplejpeg==1.7.2 -simplejson==3.19.2 -simpleline==1.8.3 -six==1.15.0 -sniffio==1.3.1 -snowballstemmer==2.2.0 -sockjs-tornado==1.0.7 -sortedcontainers==2.4.0 -sos==4.7.1 -soupsieve==2.5 -Sphinx==7.3.7 -sphinxcontrib-applehelp==1.0.8 -sphinxcontrib-devhelp==1.0.6 -sphinxcontrib-htmlhelp==2.0.5 -sphinxcontrib-jsmath==1.0.1 -sphinxcontrib-qthelp==1.0.7 -sphinxcontrib-serializinghtml==1.1.10 -stack-data==0.6.3 -subscription-manager==1.29.40 -superqt==0.6.7 -systemd-python==234 -tabulate==0.9.0 -targetcli-fb==2.1.57 -tblib==3.0.0 -tenacity==8.2.3 -terminado==0.18.1 -tifffile==2024.7.2 -tinycss2==1.2.1 -tk==0.1.0 -tomli==2.0.1 -tomli_w==1.0.0 -toolz==0.12.1 -tornado==6.4 -tqdm==4.66.2 -traitlets==5.14.2 -triangle==20230923 -trio==0.24.0 -trio-websocket==0.11.1 -tsx==0.1.16 -txt2tags==3.9 -typer==0.12.3 -types-python-dateutil==2.9.0.20240315 -typing_extensions==4.10.0 -tzdata==2024.1 -uri-template==1.3.0 -urllib3==1.26.5 -urwid==2.1.2 -vispy==0.14.3 -wcwidth==0.2.13 -webargs==8.4.0 -webcolors==1.13 -webencodings==0.5.1 -websocket-client==1.7.0 -Werkzeug==3.0.3 -widgetsnbextension==4.0.10 -wrapt==1.16.0 -wsproto==1.2.0 -xarray==2024.6.0 -yarl==1.9.4 -zarr==2.11.3 -zfpc==0.1.2 -zfpy==1.0.0 -zict==3.0.0 -zipp==3.17.0 -zope.event==5.0 -zope.interface==6.2 -zstandard==0.22.0 diff --git a/list_clean.dat b/list_clean.dat deleted file mode 100644 index 9319c35..0000000 --- a/list_clean.dat +++ /dev/null @@ -1,81 +0,0 @@ - -bcc==0.28.0 -blivet==3.6.0 -Brlapi==0.8.2 -chardet==4.0.0 -charset-normalizer==3.3.2 -cockpit @ file:///builddir/build/BUILD/cockpit-311.2/tmp/wheel/cockpit-311.2-py3-none-any.whl -configshell-fb==1.1.30 -cupshelpers==1.0 -dasbus==1.4 -dbus-python==1.2.18 -decorator==4.4.2 -distro==1.5.0 -ethtool==0.15 -file-magic==0.4.0 -gpg==1.15.1 -idna==2.10 -iniparse==0.4 -iotop==0.6 -isc==2.0 -jsonpointer==2.0 -katello-host-tools==4.2.3 -kmod==0.1 -langtable==0.0.54 -libcomps==0.1.18 -libvirt-python==10.0.0 -lxml==4.6.5 -matplotlib-inline==0.1.5 -netaddr==0.8.0 -nftables==0.1 -numpy==1.20.1 -nvmetcli==0.7 -olefile==0.46 -packaging==20.9 -pcp==5.0 -perf==0.1 -pexpect==4.8.0 -pid==2.2.3 -Pillow==10.0.1 -ply==3.11 -productmd==1.31 -psutil==5.8.0 -psycopg2==2.8.6 -ptyprocess==0.6.0 -pwquality==1.4.4 -pycairo==1.20.1 -pycups==2.0.1 -pycurl==7.43.0.6 -pyenchant==3.2.0 -PyGObject==3.40.1 -pyinotify==0.9.6 -pykickstart==3.32.11 -pyodbc===4.0.0-unsupported -pyparsing==2.4.7 -pyparted==3.12.0 -PySocks==1.7.1 -python-dateutil==2.8.1 -python-linux-procfs==0.7.3 -python-meh==0.50 -pytz==2021.1 -pyudev==0.22.0 -pyxdg==0.27 -PyYAML==5.4.1 -requests==2.25.1 -requests-file==1.5.1 -requests-ftp==0.3.1 -rpm==4.16.1.3 -rtslib-fb==2.1.76 -selinux==3.6 -sepolicy==3.6 -setools==4.4.4 -setroubleshoot @ file:///builddir/build/BUILD/setroubleshoot-3.3.32/src -simpleline==1.8.3 -six==1.15.0 -sos==4.7.1 -subscription-manager==1.29.40 -systemd-python==234 -targetcli-fb==2.1.57 -traitlets==5.1.1 -urllib3==1.26.5 -urwid==2.1.2 diff --git a/out b/out deleted file mode 100644 index 1522a2d..0000000 --- a/out +++ /dev/null @@ -1,271 +0,0 @@ -usage: tomocupy recon [-h] [--binning {0,1,2,3}] - [--blocked-views BLOCKED_VIEWS] [--dark-file-name PATH] - [--file-name PATH] [--file-type {standard,double_fov}] - [--flat-file-name PATH] [--out-path-name PATH] - [--remove-stripe-method {none,fw,ti,vo-all}] - [--bright-ratio BRIGHT_RATIO] - [--center-search-step CENTER_SEARCH_STEP] - [--center-search-width CENTER_SEARCH_WIDTH] - [--clear-folder {True,False}] [--dezinger DEZINGER] - [--dezinger-threshold DEZINGER_THRESHOLD] - [--dtype {float32,float16}] [--end-column END_COLUMN] - [--end-proj END_PROJ] [--end-row END_ROW] - [--fbp-filter {none,ramp,shepp,hann,hamming,parzen,cosine,cosine2}] - [--find-center-end-row FIND_CENTER_END_ROW] - [--find-center-start-row FIND_CENTER_START_ROW] - [--flat-linear FLAT_LINEAR] - [--max-read-threads MAX_READ_THREADS] - [--max-write-threads MAX_WRITE_THREADS] - [--minus-log MINUS_LOG] - [--nproj-per-chunk NPROJ_PER_CHUNK] [--nsino NSINO] - [--nsino-per-chunk NSINO_PER_CHUNK] - [--pixel-size PIXEL_SIZE] - [--rotation-axis ROTATION_AXIS] - [--rotation-axis-auto {manual,auto}] - [--rotation-axis-method {sift,vo}] - [--rotation-axis-pairs ROTATION_AXIS_PAIRS] - [--rotation-axis-sift-threshold ROTATION_AXIS_SIFT_THRESHOLD] - [--save-format {tiff,h5,h5sino,h5nolinks}] - [--start-column START_COLUMN] [--start-proj START_PROJ] - [--start-row START_ROW] - [--fw-filter {haar,db5,sym5,sym16}] - [--fw-level FW_LEVEL] [--fw-pad] [--fw-sigma FW_SIGMA] - [--ti-beta TI_BETA] [--ti-mask TI_MASK] - [--vo-all-dim VO_ALL_DIM] - [--vo-all-la-size VO_ALL_LA_SIZE] - [--vo-all-sm-size VO_ALL_SM_SIZE] - [--vo-all-snr VO_ALL_SNR] [--lamino-angle LAMINO_ANGLE] - [--lamino-end-row LAMINO_END_ROW] - [--lamino-search-step LAMINO_SEARCH_STEP] - [--lamino-search-width LAMINO_SEARCH_WIDTH] - [--lamino-start-row LAMINO_START_ROW] - [--reconstruction-algorithm {fourierrec,lprec,linerec}] - [--reconstruction-type {full,try}] - [--b-storage-ring B_STORAGE_RING] - [--beam-hardening-method {none,standard}] - [--calculate-source {none,standard}] - [--e-storage-ring E_STORAGE_RING] - [--filter-1-auto FILTER_1_AUTO] - [--filter-1-density FILTER_1_DENSITY] - [--filter-1-material FILTER_1_MATERIAL] - [--filter-1-thickness FILTER_1_THICKNESS] - [--filter-2-auto FILTER_2_AUTO] - [--filter-2-density FILTER_2_DENSITY] - [--filter-2-material FILTER_2_MATERIAL] - [--filter-2-thickness FILTER_2_THICKNESS] - [--filter-3-auto FILTER_3_AUTO] - [--filter-3-density FILTER_3_DENSITY] - [--filter-3-material FILTER_3_MATERIAL] - [--filter-3-thickness FILTER_3_THICKNESS] - [--maximum-E MAXIMUM_E] - [--maximum-psi-urad MAXIMUM_PSI_URAD] - [--minimum-E MINIMUM_E] [--read-pixel-size] - [--read-scintillator] [--sample-density SAMPLE_DENSITY] - [--sample-material SAMPLE_MATERIAL] - [--scintillator-density SCINTILLATOR_DENSITY] - [--scintillator-material SCINTILLATOR_MATERIAL] - [--scintillator-thickness SCINTILLATOR_THICKNESS] - [--source-distance SOURCE_DISTANCE] [--step-E STEP_E] - [--config FILE] [--config-update] [--logs-home FILE] - [--verbose] - -optional arguments: - -h, --help show this help message and exit - --binning {0,1,2,3} Reconstruction binning factor as power(2, choice) - (default: 0) - --blocked-views BLOCKED_VIEWS - Angle range for blocked views [st,end]. Can be a list - of ranges(e.g. [[0,1.2],[3,3.14]]) (default: none) - --dark-file-name PATH - Name of the hdf file containing dark data (default: - None) - --file-name PATH Name of the last used hdf file or directory containing - multiple hdf files (default: .) - --file-type {standard,double_fov} - Input file type (default: standard) - --flat-file-name PATH - Name of the hdf file containing flat data (default: - None) - --out-path-name PATH Path for output files (default: None) - --remove-stripe-method {none,fw,ti,vo-all} - Remove stripe method: none, fourier-wavelet, titarenko - (default: none) - --bright-ratio BRIGHT_RATIO - exposure time for flat fields divided by the exposure - time of projections (default: 1) - --center-search-step CENTER_SEARCH_STEP - +/- center search step (pixel). (default: 0.5) - --center-search-width CENTER_SEARCH_WIDTH - +/- center search width (pixel). (default: 50.0) - --clear-folder {True,False} - Clear output folder before reconstruction (default: - False) - --dezinger DEZINGER Width of region for removing outliers (default: 0) - --dezinger-threshold DEZINGER_THRESHOLD - Threshold of grayscale above local median to be - considered a zinger pixel (default: 5000) - --dtype {float32,float16} - Data type used for reconstruction. Note float16 works - with power of 2 sizes. (default: float32) - --end-column END_COLUMN - End position in x (default: -1) - --end-proj END_PROJ End projection (default: -1) - --end-row END_ROW End slice (default: -1) - --fbp-filter {none,ramp,shepp,hann,hamming,parzen,cosine,cosine2} - Filter for FBP reconstruction (default: parzen) - --find-center-end-row FIND_CENTER_END_ROW - End row to find the rotation center (default: -1) - --find-center-start-row FIND_CENTER_START_ROW - Start row to find the rotation center (default: 0) - --flat-linear FLAT_LINEAR - Interpolate flat fields for each projections, assumes - the number of flat fields at the beginning of the scan - is as the same as a the end. (default: False) - --max-read-threads MAX_READ_THREADS - Max number of threads for reading by chunks (default: - 4) - --max-write-threads MAX_WRITE_THREADS - Max number of threads for writing by chunks (default: - 8) - --minus-log MINUS_LOG - Take -log or not (default: True) - --nproj-per-chunk NPROJ_PER_CHUNK - Number of sinograms per chunk. Use lower numbers with - computers with lower GPU memory. (default: 8) - --nsino NSINO Location of the sinogram used for slice reconstruction - and find axis (0 top, 1 bottom). Can be given as a - list, e.g. [0,0.9]. (default: 0.5) - --nsino-per-chunk NSINO_PER_CHUNK - Number of sinograms per chunk. Use larger numbers with - computers with larger memory. (default: 8) - --pixel-size PIXEL_SIZE - Pixel size [microns] (default: 0) - --rotation-axis ROTATION_AXIS - Location of rotation axis (default: -1.0) - --rotation-axis-auto {manual,auto} - How to get rotation axis auto calculate ('auto'), or - manually ('manual') (default: manual) - --rotation-axis-method {sift,vo} - Method for automatic rotation search. (default: sift) - --rotation-axis-pairs ROTATION_AXIS_PAIRS - Projection pairs to find rotation axis. Each second - projection in a pair will be flipped and used to find - shifts from the first element in a pair. The shifts - are used to calculate the center. Example [0,1499] for - a 180 deg scan, or [0,1499,749,2249] for 360, etc. - (default: [0,0]) - --rotation-axis-sift-threshold ROTATION_AXIS_SIFT_THRESHOLD - SIFT threshold for rotation search. (default: 0.5) - --save-format {tiff,h5,h5sino,h5nolinks} - Output format (default: tiff) - --start-column START_COLUMN - Start position in x (default: 0) - --start-proj START_PROJ - Start projection (default: 0) - --start-row START_ROW - Start slice (default: 0) - --fw-filter {haar,db5,sym5,sym16} - Fourier-Wavelet remove stripe filter (default: sym16) - --fw-level FW_LEVEL Fourier-Wavelet remove stripe level parameter - (default: 7) - --fw-pad When set, Fourier-Wavelet remove stripe extend the - size of the sinogram by padding with zeros (default: - True) - --fw-sigma FW_SIGMA Fourier-Wavelet remove stripe damping parameter - (default: 1) - --ti-beta TI_BETA Parameter for ring removal (0,1) (default: 0.022) - --ti-mask TI_MASK Mask size for ring removal (0,1) (default: 1) - --vo-all-dim VO_ALL_DIM - Dimension of the window. (default: 1) - --vo-all-la-size VO_ALL_LA_SIZE - Window size of the median filter to remove large - stripes. (default: 61) - --vo-all-sm-size VO_ALL_SM_SIZE - Window size of the median filter to remove small-to- - medium stripes. (default: 21) - --vo-all-snr VO_ALL_SNR - Ratio used to locate large stripes. Greater is less - sensitive. (default: 3) - --lamino-angle LAMINO_ANGLE - Pitch of the stage for laminography (default: 0) - --lamino-end-row LAMINO_END_ROW - End slice for lamino reconstruction (default: -1) - --lamino-search-step LAMINO_SEARCH_STEP - +/- center search step (pixel). (default: 0.25) - --lamino-search-width LAMINO_SEARCH_WIDTH - +/- center search width (pixel). (default: 5.0) - --lamino-start-row LAMINO_START_ROW - Start slice for lamino reconstruction (default: 0) - --reconstruction-algorithm {fourierrec,lprec,linerec} - Reconstruction algorithm (default: fourierrec) - --reconstruction-type {full,try} - Reconstruct full data set. (default: try) - --b-storage-ring B_STORAGE_RING - Magnetic field for BM source in T (default: 0.599) - --beam-hardening-method {none,standard} - Beam hardening method. (default: none) - --calculate-source {none,standard} - Use tabulated (none, default) or calculated source - (default: none) - --e-storage-ring E_STORAGE_RING - e-beam energy for BM source in GeV (default: 7.0) - --filter-1-auto FILTER_1_AUTO - If True, read filter 1 from HDF meta data (default: - False) - --filter-1-density FILTER_1_DENSITY - Filter 1 density in g/cm^3 (default: 1.0) - --filter-1-material FILTER_1_MATERIAL - Filter 1 material for beam hardening (default: none) - --filter-1-thickness FILTER_1_THICKNESS - Filter 1 thickness in microns (default: 0.0) - --filter-2-auto FILTER_2_AUTO - If True, read filter 2 from HDF meta data (default: - False) - --filter-2-density FILTER_2_DENSITY - Filter 2 density in g/cm^3 (default: 1.0) - --filter-2-material FILTER_2_MATERIAL - Filter 2 material for beam hardening (default: none) - --filter-2-thickness FILTER_2_THICKNESS - Filter 2 thickness in microns (default: 0.0) - --filter-3-auto FILTER_3_AUTO - If True, read filter 3 from HDF meta data (default: - False) - --filter-3-density FILTER_3_DENSITY - Filter 3 density in g/cm^3 (default: 1.0) - --filter-3-material FILTER_3_MATERIAL - Filter 3 material in microns (default: none) - --filter-3-thickness FILTER_3_THICKNESS - Filter 3 thickness for beam hardening (default: 0.0) - --maximum-E MAXIMUM_E - Maximum energy to model in eV (default: 200000) - --maximum-psi-urad MAXIMUM_PSI_URAD - Maximum vertical angle from centerline to model in - microradians (default: 40) - --minimum-E MINIMUM_E - Minimum energy to model in eV (default: 1000) - --read-pixel-size When set, read effective pixel size from the HDF file - (default: False) - --read-scintillator When set, read scintillator properties from the HDF - file (default: False) - --sample-density SAMPLE_DENSITY - Density of sample material in g/cm^3 (default: 1.0) - --sample-material SAMPLE_MATERIAL - Sample material for beam hardening (default: Fe) - --scintillator-density SCINTILLATOR_DENSITY - Density of scintillator in g/cm^3 (default: 6.0) - --scintillator-material SCINTILLATOR_MATERIAL - Scintillator material for beam hardening (default: - LuAG_Ce) - --scintillator-thickness SCINTILLATOR_THICKNESS - Scintillator thickness in microns (default: 100.0) - --source-distance SOURCE_DISTANCE - Distance from source to scintillator in m (default: - 36.0) - --step-E STEP_E Energy step in eV (default: 500) - --config FILE File name of configuration file (default: - /home/beams/AMITTONE/tomocupyon.conf) - --config-update When set, the content of the config file is updated - using the current params values (default: False) - --logs-home FILE Log file directory (default: - /home/beams/AMITTONE/logs) - --verbose Verbose output (default: False) From e8c3abb7da900b5409cf71600b0b9ea10acbdbc7 Mon Sep 17 00:00:00 2001 From: Viktor Nikitin Date: Mon, 11 Nov 2024 14:34:44 -0600 Subject: [PATCH 08/21] add zarr dependency --- VERSION | 2 +- recipe/meta.yaml | 3 ++- setup.py | 1 + src/tomocupy/dataio/reader.py | 1 - 4 files changed, 4 insertions(+), 3 deletions(-) diff --git a/VERSION b/VERSION index a6a3a43..ece61c6 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.0.4 \ No newline at end of file +1.0.6 \ No newline at end of file diff --git a/recipe/meta.yaml b/recipe/meta.yaml index e26bb4e..ee6d203 100644 --- a/recipe/meta.yaml +++ b/recipe/meta.yaml @@ -41,7 +41,8 @@ requirements: - opencv - opencv-contrib-python - pywavelets - - tifffile + - tifffile + - zarr about: home: https://github.com/tomography/tomocupy diff --git a/setup.py b/setup.py index 72c2c97..d21f504 100644 --- a/setup.py +++ b/setup.py @@ -18,5 +18,6 @@ 'pywavelets', 'setuptools', # for pkg_resources at runtime 'tifffile', + 'zarr', ] ) diff --git a/src/tomocupy/dataio/reader.py b/src/tomocupy/dataio/reader.py index cd3abca..e047bd2 100644 --- a/src/tomocupy/dataio/reader.py +++ b/src/tomocupy/dataio/reader.py @@ -84,7 +84,6 @@ def init_sizes(self): # read data sizes and projection angles with a reader sizes = self.read_sizes() - print(sizes['nproji']) theta = self.read_theta(sizes['nproji']) nproji = sizes['nproji'] nzi = sizes['nzi'] From a2b57194641c8768141e2afccc792c6842ea1d32 Mon Sep 17 00:00:00 2001 From: Viktor Nikitin Date: Mon, 11 Nov 2024 14:35:17 -0600 Subject: [PATCH 09/21] checking write with zarr --- src/tomocupy/dataio/writer.py | 3 +- tests/test_full.py | 73 +++++++++++++++++++++-------------- 2 files changed, 45 insertions(+), 31 deletions(-) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index b2d2c2e..0c3af02 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -265,7 +265,8 @@ def write_data_chunk(self, rec, st, end, k): chunks=(params.nproj, 1, params.n)) elif args.save_format == 'zarr': # chunks = (1, params.n, params.n) # Replace with appropriate chunk size - save_zarr(volume=rec, output_path=self.zarr_output_path, chunks=chunks, compression=args.zarr_compression, pixel_size=args.pixel_size) + print(f'save chunk {rec[:end-st].shape}') + save_zarr(volume=rec[:end-st], output_path=self.zarr_output_path, chunks=chunks, compression=args.zarr_compression, pixel_size=args.pixel_size) def write_data_try(self, rec, cid, id_slice): """Write tiff reconstruction with a given name""" diff --git a/tests/test_full.py b/tests/test_full.py index 1d282c7..54b89df 100644 --- a/tests/test_full.py +++ b/tests/test_full.py @@ -5,6 +5,7 @@ import tifffile import inspect import h5py +import zarr import shutil prefix = 'tomocupy recon --file-name data/test_data.h5 --reconstruction-type full --rotation-axis 782.5 --nsino-per-chunk 4' @@ -13,35 +14,36 @@ prefix4 = '--filter-1-density 1.85 --filter-2-density 8.9 --filter-3-density 8.9' prefix5 = '--filter-1-density 0.0 --filter-2-density 0.0 --filter-3-density 0.0' cmd_dict = { - f'{prefix} ': 28.307, - f'{prefix} --reconstruction-algorithm lprec ': 27.992, - f'{prefix} --reconstruction-algorithm linerec ': 28.341, - f'{prefix} --dtype float16': 24.186, - f'{prefix} --reconstruction-algorithm lprec --dtype float16': 24.050, - f'{prefix} --reconstruction-algorithm linerec --dtype float16': 25.543, - f'{prefix} --binning 1': 12.286, - f'{prefix} --reconstruction-algorithm lprec --binning 1': 12.252, - f'{prefix} --reconstruction-algorithm linerec --binning 1': 12.259, - f'{prefix} --start-row 3 --end-row 15 --start-proj 200 --end-proj 700': 17.589, - f'{prefix} --save-format h5': 28.307, - f'{prefix} --nsino-per-chunk 2 --file-type double_fov': 15.552, - f'{prefix} --nsino-per-chunk 2 --blocked-views [0.2,1]': 30.790, - f'{prefix} --nsino-per-chunk 2 --blocked-views [[0.2,1],[2,3]]': 40.849, - f'{prefix} --remove-stripe-method fw': 28.167, - f'{prefix} --remove-stripe-method fw --dtype float16': 23.945, - f'{prefix} --start-column 200 --end-column 1000': 18.248, - f'{prefix} --start-column 200 --end-column 1000 --binning 1': 7.945, - f'{prefix} --flat-linear True': 28.308, - f'{prefix} --rotation-axis-auto auto --rotation-axis-method sift --reconstruction-type full' : 28.305, - f'{prefix} --rotation-axis-auto auto --rotation-axis-method vo --center-search-step 0.1 --nsino 0.5 --center-search-width 100 --reconstruction-type full' : 28.303, - f'{prefix} --remove-stripe-method vo-all ': 27.993, - f'{prefix} --bright-ratio 10': 32.631, - f'{prefix} --end-column 1535': 28.293, - f'{prefix} --end-column 1535 --binning 3': 1.82, - f'{prefix2} {prefix3} {prefix5} --beam-hardening-method standard --calculate-source standard': 3255.912, - f'{prefix2} {prefix3} {prefix4} --beam-hardening-method standard': 3248.832, - f'{prefix2} {prefix3} {prefix4} --beam-hardening-method standard --calculate-source standard': 3254.634, - f'{prefix2} {prefix3} {prefix4} --beam-hardening-method standard --calculate-source standard --e-storage-ring 3.0 --b-storage-ring 0.3': 822.178, + # f'{prefix} ': 28.307, + # f'{prefix} --reconstruction-algorithm lprec ': 27.992, + # f'{prefix} --reconstruction-algorithm linerec ': 28.341, + # f'{prefix} --dtype float16': 24.186, + # f'{prefix} --reconstruction-algorithm lprec --dtype float16': 24.050, + # f'{prefix} --reconstruction-algorithm linerec --dtype float16': 25.543, + # f'{prefix} --binning 1': 12.286, + # f'{prefix} --reconstruction-algorithm lprec --binning 1': 12.252, + # f'{prefix} --reconstruction-algorithm linerec --binning 1': 12.259, + # f'{prefix} --start-row 3 --end-row 15 --start-proj 200 --end-proj 700': 17.589, + # f'{prefix} --save-format h5': 28.307, + # f'{prefix} --nsino-per-chunk 2 --file-type double_fov': 15.552, + # f'{prefix} --nsino-per-chunk 2 --blocked-views [0.2,1]': 30.790, + # f'{prefix} --nsino-per-chunk 2 --blocked-views [[0.2,1],[2,3]]': 40.849, + # f'{prefix} --remove-stripe-method fw': 28.167, + # f'{prefix} --remove-stripe-method fw --dtype float16': 23.945, + # f'{prefix} --start-column 200 --end-column 1000': 18.248, + # f'{prefix} --start-column 200 --end-column 1000 --binning 1': 7.945, + # f'{prefix} --flat-linear True': 28.308, + # f'{prefix} --rotation-axis-auto auto --rotation-axis-method sift --reconstruction-type full' : 28.305, + # f'{prefix} --rotation-axis-auto auto --rotation-axis-method vo --center-search-step 0.1 --nsino 0.5 --center-search-width 100 --reconstruction-type full' : 28.303, + # f'{prefix} --remove-stripe-method vo-all ': 27.993, + # f'{prefix} --bright-ratio 10': 32.631, + # f'{prefix} --end-column 1535': 28.293, + # f'{prefix} --end-column 1535 --binning 3': 1.82, + # f'{prefix2} {prefix3} {prefix5} --beam-hardening-method standard --calculate-source standard': 3255.912, + # f'{prefix2} {prefix3} {prefix4} --beam-hardening-method standard': 3248.832, + # f'{prefix2} {prefix3} {prefix4} --beam-hardening-method standard --calculate-source standard': 3254.634, + # f'{prefix2} {prefix3} {prefix4} --beam-hardening-method standard --calculate-source standard --e-storage-ring 3.0 --b-storage-ring 0.3': 822.178, + f'{prefix} --save-format zarr': 28.307, } class SequentialTestLoader(unittest.TestLoader): @@ -73,7 +75,6 @@ def test_full_recon(self): data_file = Path('data_rec').joinpath(file_name) with h5py.File('data_rec/test_data_rec.h5', 'r') as fid: data = fid['exchange/data'] - print(data.shape) ssum = np.sum(np.linalg.norm(data[:], axis=(1, 2))) except: pass @@ -84,6 +85,18 @@ def test_full_recon(self): f'data_rec/{file_name}_rec/recon_{k:05}.tiff')) except: pass + #try: + import time + time.sleep(1) + file_name = cmd[0].split("--file-name ")[1].split('.')[0].split('/')[-1] + data_file = Path('data_rec').joinpath(file_name) + with zarr.open('data_rec/test_data_rec.zarr','r') as fid: + print(fid[0][:].shape,fid[1][:].shape) + data = fid[0][:].copy() + print(data.shape) + ssum = np.sum(np.linalg.norm(data[:], axis=(1, 2))) + #except: + #pass self.assertAlmostEqual(ssum, cmd[1], places=0) From af7e1848ff1f99358d496153ee192433fa9e0039 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Thu, 21 Nov 2024 16:06:18 -0600 Subject: [PATCH 10/21] solved bug in zarr data saving --- src/tomocupy/config.py | 2 +- src/tomocupy/dataio/writer.py | 198 ++++++++++++++++++++++------------ src/tomocupy/rec.py | 2 +- src/tomocupy/utils.py | 23 +++- tests/test_full.py | 30 ++++-- 5 files changed, 171 insertions(+), 84 deletions(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index a09c992..fbbae33 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -797,7 +797,7 @@ def update_hdf_process(fname, args=None, sections=None): hdf_file.require_dataset(dataset, shape=(1,), dtype=dt) log.info(name + ': ' + str(value)) try: - hdf_file[dataset][0] = np.string_(str(value)) + hdf_file[dataset][0] = np.bytes_(str(value)) except TypeError: log.error( "Could not convert value {}".format(value)) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 0c3af02..b44fe1b 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -52,8 +52,11 @@ import zarr import json from numcodecs import Blosc -from tomocupy.utils import downsampleZarr +import threading +import subprocess +import time +from skimage.transform import downscale_local_mean __author__ = "Viktor Nikitin" __copyright__ = "Copyright (c) 2022, UChicago Argonne, LLC." @@ -213,8 +216,9 @@ def init_output_files(self): if args.save_format == 'zarr': # Zarr format support fnameout += '.zarr' self.zarr_output_path = fnameout + clean_zarr(self.zarr_output_path) log.info(f'Zarr dataset will be created at {fnameout}') - + params.fnameout = fnameout log.info(f'Output: {fnameout}') @@ -244,7 +248,8 @@ def write_meta(self, rec_virtual): log.error('write_meta() error: Skip copying meta') pass - def write_data_chunk(self, rec, st, end, k): + + def write_data_chunk(self, rec, st, end, k, shift_index): """Writing the kth data chunk to hard disk""" if args.save_format == 'tiff': @@ -263,85 +268,140 @@ def write_data_chunk(self, rec, st, end, k): with h5py.File(filename, "w") as fid: fid.create_dataset("/exchange/data", data=rec, chunks=(params.nproj, 1, params.n)) - elif args.save_format == 'zarr': # - chunks = (1, params.n, params.n) # Replace with appropriate chunk size - print(f'save chunk {rec[:end-st].shape}') - save_zarr(volume=rec[:end-st], output_path=self.zarr_output_path, chunks=chunks, compression=args.zarr_compression, pixel_size=args.pixel_size) + elif args.save_format == 'zarr': # Zarr format support + # Define chunk size + chunks = (1, params.n, params.n) + + if not hasattr(self, 'zarr_array'): + shape = (int(params.nz / 2**args.binning), params.n, params.n) # Full dataset shape + + max_levels = lambda X, Y: (lambda r: (int(r).bit_length() - 1) if r != 0 else (int(X // Y).bit_length() - 1))(int(X) % int(Y)) + levels = min(max_levels(params.nz, end-st),6) + + self.zarr_array = initialize_zarr( + output_path=self.zarr_output_path, + base_shape=shape, + chunks=chunks, + dtype=params.dtype, + num_levels=levels, + compression=args.zarr_compression + ) + + # Write the current chunk to the Zarr container + write_zarr_chunk( + zarr_group=self.zarr_array, # Pre-initialized Zarr container + data_chunk=rec[:end - st], # Data chunk to save + start=st-shift_index, # Starting index for this chunk along the z-axis + end=end-shift_index # Ending index for this chunk along the z-axis + ) + def write_data_try(self, rec, cid, id_slice): """Write tiff reconstruction with a given name""" tifffile.imwrite( f'{params.fnameout}_slice{id_slice:04d}_center{cid:05.2f}.tiff', rec) + - -def save_zarr(volume, output_path, chunks, compression, pixel_size, mode='a', original_dtype=np.float32): +def clean_zarr(output_path): + if os.path.exists(output_path): + try: + subprocess.run(["rm", "-rf", output_path], check=True) + log.info(f"Successfully removed directory: {output_path}") + except subprocess.CalledProcessError as e: + log.error(f"Error removing directory {output_path}: {e}") + raise + else: + log.warning(f"Path does not exist: {output_path}") + + +def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, compression='blosclz'): """ - Save a 3D volume to a Zarr store, creating a multiscale pyramid representation. - + Initialize a multiscale Zarr container with specified levels, dimensions, and compression. + Parameters: - - volume (numpy array): The 3D volume data to be saved. - - output_path (str): The path to the output Zarr store. - - chunks (tuple of ints): The chunk size for the Zarr array. - - compression (str): The compression algorithm to use (e.g., 'blosclz', 'lz4', etc.). - - pixel_size (float): The size of the pixels in micrometers. - - mode (str, optional): The mode to open the Zarr store ('w' for write, 'a' for append). Default is 'a'. - - original_dtype (numpy dtype, optional): The original data type of the images. Default is np.uint8. + - output_path (str): Path to the Zarr file. + - base_shape (tuple): Shape of the full dataset at the highest resolution (e.g., (z, y, x)). + - chunks (tuple): Chunk size for the dataset (e.g., (1, y, x)). + - dtype: Data type of the dataset (e.g., np.float32). + - num_levels (int): Number of multiresolution levels. + - compression (str): Compression algorithm (default 'blosclz'). Returns: - - None - """ + - zarr.Group: The initialized Zarr group containing multiscale datasets. + """ store = zarr.DirectoryStore(output_path) compressor = Blosc(cname=compression, clevel=5, shuffle=2) + root_group = zarr.group(store=store) - if mode == 'w': - if os.path.exists(output_path): - shutil.rmtree(output_path) - root_group = zarr.group(store=store) - else: - root_group = zarr.open(store=store, mode='a') - - # Assuming volume is a chunk of the data - pyramid_levels = downsampleZarr(volume) # Assuming downsample is defined elsewhere for multiscale - - datasets = [] - for level, data in enumerate(pyramid_levels): - data = data.astype(original_dtype) + current_shape = base_shape + for level in range(num_levels): + level_name = f"{level}" + log.info(f"Initializing level {level} with shape {current_shape}") - dataset_name = f"{level}" - if dataset_name in root_group: - z = root_group[dataset_name] - z.append(data, axis=0) - else: - z = root_group.create_dataset(name=dataset_name, shape=data.shape, chunks=chunks, dtype=data.dtype, compressor=compressor) - z[:] = data + # Dynamically scale chunks based on the level + level_chunks = tuple(max(1, c // (2 ** level)) for c in chunks) - scale_factor = 2 ** level - datasets.append({ - "path": dataset_name, - "coordinateTransformations": [{"type": "scale", "scale": [pixel_size * scale_factor, pixel_size * scale_factor, pixel_size * scale_factor]}] - }) - - if mode == 'w': - multiscales = [{ - "version": "0.4", - "name": "example", - "axes": [ - {"name": "z", "type": "space", "unit": "micrometer"}, - {"name": "y", "type": "space", "unit": "micrometer"}, - {"name": "x", "type": "space", "unit": "micrometer"} - ], - "datasets": datasets, - "type": "gaussian", - "metadata": { - "method": "skimage.transform.resize", - "version": "0.16.1", - "args": "[true]", - "kwargs": {"anti_aliasing": True, "preserve_range": True} - } - }] - - root_group.attrs.update({"multiscales": multiscales}) - with open(os.path.join(output_path, 'multiscales.json'), 'w') as f: - json.dump({"multiscales": multiscales}, f, indent=2) - #info(f"Metadata saved to {output_path}") + # Create dataset for the current level + root_group.create_dataset( + name=level_name, + shape=current_shape, + chunks=level_chunks, + dtype=dtype, + compressor=compressor + ) + current_shape = tuple(max(1, s // 2) for s in current_shape) + + return root_group + +def write_zarr_chunk(zarr_group, data_chunk, start, end): + """ + Write a chunk of data into the Zarr container at all resolutions. + + Parameters: + - zarr_group (zarr.Group): The initialized Zarr group containing multiscale datasets. + - data_chunk (np.ndarray): The data chunk to write (highest resolution). + - start (int): Start index in the first dimension (z-axis) for the highest resolution. + - end (int): End index in the first dimension (z-axis) for the highest resolution. + """ + for level in sorted(zarr_group.keys(), key=int): # Process levels in order (0, 1, ...) + zarr_array = zarr_group[level] # Access the dataset for the current level + + # Calculate the downscaling factor for this resolution level + scale_factor = 2 ** int(level) + # Downsample data chunk for this resolution level + if scale_factor > 1: + downsampled_chunk = downsample_volume(data_chunk, scale_factor) + else: + downsampled_chunk = data_chunk + + level_start = start // scale_factor + level_end = end // scale_factor + + expected_z_size = level_end - level_start + actual_z_size = downsampled_chunk.shape[0] + + if actual_z_size != expected_z_size: + downsampled_chunk = downsampled_chunk[:expected_z_size] + + # Write the downsampled chunk into the Zarr dataset + zarr_array[level_start:level_end, :, :] = downsampled_chunk + #log.info(f"Saved chunk to level {level} [{level_start}:{level_end}] with shape {downsampled_chunk.shape}") + + + +def downsample_volume(volume, scale_factor): + """ + Downsample a 3D volume by a given scale factor. + + Parameters: + - volume (numpy array): Input 3D volume. + - scale_factor (int): Factor by which to downsample (e.g., 2 for halving). + + Returns: + - numpy array: Downsampled volume. + """ + if scale_factor == 1: + return volume # No downsampling needed for the highest resolution + factors = (1, scale_factor, scale_factor) # Only downsample spatial dimensions + return downscale_local_mean(volume, factors) diff --git a/src/tomocupy/rec.py b/src/tomocupy/rec.py index aeed5e8..4f7090b 100644 --- a/src/tomocupy/rec.py +++ b/src/tomocupy/rec.py @@ -194,7 +194,7 @@ def recon_all(self): st = ids[k-2]*ncz+args.start_row//2**args.binning end = st+lzchunk[ids[k-2]] self.write_threads[ithread].run( - self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, ids[k-2])) + self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, ids[k-2], args.start_row)) self.stream1.synchronize() self.stream2.synchronize() diff --git a/src/tomocupy/utils.py b/src/tomocupy/utils.py index 9bae74c..396fb17 100644 --- a/src/tomocupy/utils.py +++ b/src/tomocupy/utils.py @@ -52,7 +52,7 @@ from skimage.transform import resize import time from functools import wraps - +import subprocess @@ -276,7 +276,7 @@ def param_from_dxchange(hdf_file, data_path, attr=None, scalar=True, char_array= return None -def downsampleZarr(data, scale_factor=2, max_levels=6): +def downsampleZarr(data, scale_factor=2, max_levels=1): """ Create a multi-level downsampled version of the input data. @@ -294,9 +294,24 @@ def downsampleZarr(data, scale_factor=2, max_levels=6): current_level = data levels = [current_level] for _ in range(max_levels): - new_shape = tuple(max(1, dim // 2) for dim in current_level.shape) + new_shape = tuple(max(1, dim // 1) for dim in current_level.shape) if min(new_shape) <= 1: break - current_level = resize(current_level, new_shape, order=0, preserve_range=True, anti_aliasing=True) + current_level = current_level#resize(current_level, new_shape, order=0, preserve_range=True, anti_aliasing=True) + print(current_level.shape) levels.append(current_level) return levels + + +def clean_zarr(output_path): + if os.path.exists(output_path): + try: + subprocess.run(["mv", output_path, output_path + "_tmp"], check=True) + subprocess.run(["ls -lrt", output_path + "_tmp"], check=True) + subprocess.run(["rm", "-rf", output_path + "_tmp"], check=True) + log.info(f"Successfully removed directory: {output_path}") + except subprocess.CalledProcessError as e: + log.error(f"Error removing directory {output_path}: {e}") + raise + else: + log.warning(f"Path does not exist: {output_path}") diff --git a/tests/test_full.py b/tests/test_full.py index 54b89df..db43aa5 100644 --- a/tests/test_full.py +++ b/tests/test_full.py @@ -14,7 +14,7 @@ prefix4 = '--filter-1-density 1.85 --filter-2-density 8.9 --filter-3-density 8.9' prefix5 = '--filter-1-density 0.0 --filter-2-density 0.0 --filter-3-density 0.0' cmd_dict = { - # f'{prefix} ': 28.307, + # '{prefix} ': 28.307, # f'{prefix} --reconstruction-algorithm lprec ': 27.992, # f'{prefix} --reconstruction-algorithm linerec ': 28.341, # f'{prefix} --dtype float16': 24.186, @@ -90,16 +90,28 @@ def test_full_recon(self): time.sleep(1) file_name = cmd[0].split("--file-name ")[1].split('.')[0].split('/')[-1] data_file = Path('data_rec').joinpath(file_name) - with zarr.open('data_rec/test_data_rec.zarr','r') as fid: - print(fid[0][:].shape,fid[1][:].shape) - data = fid[0][:].copy() - print(data.shape) - ssum = np.sum(np.linalg.norm(data[:], axis=(1, 2))) - #except: - #pass - self.assertAlmostEqual(ssum, cmd[1], places=0) + # Open the Zarr dataset + fid = zarr.open('data_rec/test_data_rec.zarr', mode='r') + + # Log dataset shapes + print(fid[0][:].shape, fid[1][:].shape) + + # Access and copy data + data = fid[0][:].astype(np.float64).copy() + print(f"Data shape: {data.shape}") + print(f"Data sample: {data[:5]}") + # Normalize data + data = np.abs(data) + + # Calculate the sum of norms + ssum = np.sum(np.linalg.norm(data, axis=1)) # Adjust axis if needed + print(f"Computed ssum: {ssum}") + print(f"Expected value: {cmd[1]}") + + # Perform the test comparison + self.assertAlmostEqual(ssum, cmd[1], places=0) if __name__ == '__main__': unittest.main(testLoader=SequentialTestLoader(), failfast=True) From 824a800ccd2a91e099d670075ca8ab47d09b36fa Mon Sep 17 00:00:00 2001 From: mittoalb Date: Fri, 22 Nov 2024 14:09:01 -0600 Subject: [PATCH 11/21] recon steps zarr data saving --- src/tomocupy/reconstruction/backproj_lamfourier_parallel.py | 2 +- src/tomocupy/reconstruction/backproj_parallel.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py b/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py index a1649cd..b463906 100644 --- a/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py +++ b/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py @@ -290,7 +290,7 @@ def write_parallel(self, u, nthreads=16): st = k*nchunk+params.lamino_start_row end = min((k+1)*nchunk, u.shape[0])+params.lamino_start_row th = Thread(target=self.cl_writer.write_data_chunk, args=( - u[k*nchunk:min((k+1)*nchunk, u.shape[0])], st, end, k)) + u[k*nchunk:min((k+1)*nchunk, u.shape[0])], st, end, k, args.start_row)) mthreads.append(th) th.start() for th in mthreads: diff --git a/src/tomocupy/reconstruction/backproj_parallel.py b/src/tomocupy/reconstruction/backproj_parallel.py index f94249e..4355b01 100644 --- a/src/tomocupy/reconstruction/backproj_parallel.py +++ b/src/tomocupy/reconstruction/backproj_parallel.py @@ -159,7 +159,7 @@ def recon_sino_proj_parallel(self, data): st = (kr-2)*ncz+args.lamino_start_row//2**args.binning end = st+lrchunk[kr-2] self.write_threads[ithread].run( - self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, kr-2)) + self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, kr-2,args.start_row)) self.stream1.synchronize() self.stream2.synchronize() @@ -364,7 +364,7 @@ def recon_sino_parallel(self, data): st = (k-2)*ncz+args.start_row//2**args.binning end = st+lzchunk[k-2] self.write_threads[ithread].run( - self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, k-2)) + self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, k-2, args.start_row)) self.stream1.synchronize() self.stream2.synchronize() From 9c374cb3cdc8c2b3d8ae34a6585c2f1e950efd67 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Mon, 25 Nov 2024 17:52:40 -0600 Subject: [PATCH 12/21] Meta saving --- src/tomocupy/dataio/writer.py | 86 ++++++++++++++++++++++++++++++++--- 1 file changed, 80 insertions(+), 6 deletions(-) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index b44fe1b..40a7334 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -286,7 +286,7 @@ def write_data_chunk(self, rec, st, end, k, shift_index): num_levels=levels, compression=args.zarr_compression ) - + fill_zarr_meta(self.zarr_array, self.zarr_output_path, args) # Write the current chunk to the Zarr container write_zarr_chunk( zarr_group=self.zarr_array, # Pre-initialized Zarr container @@ -314,8 +314,79 @@ def clean_zarr(output_path): else: log.warning(f"Path does not exist: {output_path}") +from pathlib import PosixPath +from types import SimpleNamespace + + +def args2json(data): + """ + Recursively convert all unsupported types (e.g., PosixPath, Namespace) to JSON-serializable types. + + Parameters: + - data: The input data (can be dict, list, PosixPath, Namespace, etc.). + + Returns: + - A JSON-serializable version of the data. + """ + if isinstance(data, PosixPath): + return str(data) # Convert PosixPath to string + elif isinstance(data, SimpleNamespace): + return {k: args2json(v) for k, v in vars(data).items()} # Convert Namespace to dict + elif isinstance(data, dict): + return {k: args2json(v) for k, v in data.items()} # Recurse into dict + elif isinstance(data, list): + return [args2json(item) for item in data] # Recurse into list + elif isinstance(data, tuple): + return tuple(args2json(item) for item in data) # Recurse into tuple + else: + return data + + +def fill_zarr_meta(root_group, output_path, args, mode='w'): + """ + Fill metadata for the Zarr multiscale datasets and include additional parameters. -def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, compression='blosclz'): + Parameters: + - output_path (str): Path to save the metadata file. + - datasets (list): List of datasets with their metadata. + - args (dict): Additional parameters to include in the metadata. + - mode (str): Mode for metadata handling. Default is 'w'. + """ + + if not isinstance(args, dict): + metadata = vars(args) # Convert Namespace to dictionary + + multiscales = [{ + "version": "0.4", + "name": args.file_name, + "axes": [ + {"name": "z", "type": "space", "unit": "micrometer"}, + {"name": "y", "type": "space", "unit": "micrometer"}, + {"name": "x", "type": "space", "unit": "micrometer"} + ], + "type": "gaussian", + "metadata": { + "method": "skimage.transform.downscale_local_mean", + "version": "0.16.1", + "args": "[true]", + "kwargs": {"anti_aliasing": True, "preserve_range": True}, + "Reco Data": str(args) # Add additional parameters from the params dictionary + } + }] + + + multiscales = args2json(multiscales) + + if mode == 'w': + root_group.attrs.update({"multiscales": multiscales}) + + # Write human-readable JSON to file + metadata_file = os.path.join(output_path, str(args.file_name)[:-2] +'json') + with open(metadata_file, 'w') as f: + json.dump({"multiscales": multiscales}, f, indent=4) + + +def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, compression='blosclz'):#blosclz'): """ Initialize a multiscale Zarr container with specified levels, dimensions, and compression. @@ -329,11 +400,12 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, compress Returns: - zarr.Group: The initialized Zarr group containing multiscale datasets. - """ + """ + store = zarr.DirectoryStore(output_path) compressor = Blosc(cname=compression, clevel=5, shuffle=2) root_group = zarr.group(store=store) - + current_shape = base_shape for level in range(num_levels): level_name = f"{level}" @@ -349,11 +421,13 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, compress chunks=level_chunks, dtype=dtype, compressor=compressor - ) + ) + current_shape = tuple(max(1, s // 2) for s in current_shape) - + return root_group + def write_zarr_chunk(zarr_group, data_chunk, start, end): """ Write a chunk of data into the Zarr container at all resolutions. From c32cfadb95c046f9bef506eb8bc96d2ed32ff2eb Mon Sep 17 00:00:00 2001 From: mittoalb Date: Mon, 25 Nov 2024 18:30:04 -0600 Subject: [PATCH 13/21] Updated meta saving --- src/tomocupy/dataio/writer.py | 39 ++++++++++++++++++++++++++++++++--- 1 file changed, 36 insertions(+), 3 deletions(-) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 40a7334..49fa00c 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -56,7 +56,9 @@ import subprocess import time -from skimage.transform import downscale_local_mean +#from skimage.transform import downscale_local_mean +import cupy as cp + __author__ = "Viktor Nikitin" __copyright__ = "Copyright (c) 2022, UChicago Argonne, LLC." @@ -378,7 +380,7 @@ def fill_zarr_meta(root_group, output_path, args, mode='w'): multiscales = args2json(multiscales) if mode == 'w': - root_group.attrs.update({"multiscales": multiscales}) + root_group.attrs.update({"multiresolution": multiscales}) # Write human-readable JSON to file metadata_file = os.path.join(output_path, str(args.file_name)[:-2] +'json') @@ -463,7 +465,7 @@ def write_zarr_chunk(zarr_group, data_chunk, start, end): #log.info(f"Saved chunk to level {level} [{level_start}:{level_end}] with shape {downsampled_chunk.shape}") - +""" def downsample_volume(volume, scale_factor): """ Downsample a 3D volume by a given scale factor. @@ -479,3 +481,34 @@ def downsample_volume(volume, scale_factor): return volume # No downsampling needed for the highest resolution factors = (1, scale_factor, scale_factor) # Only downsample spatial dimensions return downscale_local_mean(volume, factors) +""" + + +def downsample_volume(volume, scale_factor): + """ + Downsample a 3D volume by a given scale factor using NumPy. + + Parameters: + - volume (numpy array): Input 3D volume (e.g., [z, y, x]). + - scale_factor (int): Factor by which to downsample (e.g., 2 for halving). + + Returns: + - numpy array: Downsampled volume. + """ + if scale_factor == 1: + return volume # No downsampling needed for the highest resolution + + # Ensure the input dimensions are divisible by the scale factor + if volume.shape[1] % scale_factor != 0 or volume.shape[2] % scale_factor != 0: + raise ValueError("Volume dimensions must be divisible by the scale factor.") + + # Reshape the spatial dimensions for downsampling + z, y, x = volume.shape + downsampled = volume.reshape( + z, + y // scale_factor, scale_factor, + x // scale_factor, scale_factor + ) + + # Compute the mean along the downsampling axes + return downsampled.mean(axis=(2, 4)) From 36b1d4a959eac414ec216a4b59971aec03f3b042 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Tue, 26 Nov 2024 10:38:52 -0600 Subject: [PATCH 14/21] Updated meta saving --- src/tomocupy/dataio/writer.py | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 49fa00c..2bf0f34 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -328,7 +328,7 @@ def args2json(data): - data: The input data (can be dict, list, PosixPath, Namespace, etc.). Returns: - - A JSON-serializable version of the data. + - A JSON data. """ if isinstance(data, PosixPath): return str(data) # Convert PosixPath to string @@ -465,25 +465,6 @@ def write_zarr_chunk(zarr_group, data_chunk, start, end): #log.info(f"Saved chunk to level {level} [{level_start}:{level_end}] with shape {downsampled_chunk.shape}") -""" -def downsample_volume(volume, scale_factor): - """ - Downsample a 3D volume by a given scale factor. - - Parameters: - - volume (numpy array): Input 3D volume. - - scale_factor (int): Factor by which to downsample (e.g., 2 for halving). - - Returns: - - numpy array: Downsampled volume. - """ - if scale_factor == 1: - return volume # No downsampling needed for the highest resolution - factors = (1, scale_factor, scale_factor) # Only downsample spatial dimensions - return downscale_local_mean(volume, factors) -""" - - def downsample_volume(volume, scale_factor): """ Downsample a 3D volume by a given scale factor using NumPy. From b7d9cc3a7aa27f5bbf023154fba43445028995e7 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Tue, 26 Nov 2024 19:09:19 -0600 Subject: [PATCH 15/21] Fixed zarr data neuroglancer compatibility --- src/tomocupy/dataio/writer.py | 84 ++++++++++++++++++----------------- 1 file changed, 44 insertions(+), 40 deletions(-) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 2bf0f34..2715049 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -56,8 +56,9 @@ import subprocess import time -#from skimage.transform import downscale_local_mean import cupy as cp +from pathlib import PosixPath +from types import SimpleNamespace __author__ = "Viktor Nikitin" @@ -272,23 +273,26 @@ def write_data_chunk(self, rec, st, end, k, shift_index): chunks=(params.nproj, 1, params.n)) elif args.save_format == 'zarr': # Zarr format support # Define chunk size - chunks = (1, params.n, params.n) + #chunks = (1, params.n, params.n) + chunks = (params.nz, 256, 256) if not hasattr(self, 'zarr_array'): shape = (int(params.nz / 2**args.binning), params.n, params.n) # Full dataset shape max_levels = lambda X, Y: (lambda r: (int(r).bit_length() - 1) if r != 0 else (int(X // Y).bit_length() - 1))(int(X) % int(Y)) levels = min(max_levels(params.nz, end-st),6) + scale_factors = [float(args.pixel_size) * (i + 1) for i in range(levels)] - self.zarr_array = initialize_zarr( + self.zarr_array, datasets = initialize_zarr( output_path=self.zarr_output_path, base_shape=shape, chunks=chunks, dtype=params.dtype, num_levels=levels, + scale_factors=scale_factors, compression=args.zarr_compression ) - fill_zarr_meta(self.zarr_array, self.zarr_output_path, args) + fill_zarr_meta(self.zarr_array, datasets, self.zarr_output_path, args) # Write the current chunk to the Zarr container write_zarr_chunk( zarr_group=self.zarr_array, # Pre-initialized Zarr container @@ -316,9 +320,6 @@ def clean_zarr(output_path): else: log.warning(f"Path does not exist: {output_path}") -from pathlib import PosixPath -from types import SimpleNamespace - def args2json(data): """ @@ -344,79 +345,76 @@ def args2json(data): return data -def fill_zarr_meta(root_group, output_path, args, mode='w'): +def fill_zarr_meta(root_group, datasets, output_path, metadata_args, mode='w'): """ Fill metadata for the Zarr multiscale datasets and include additional parameters. Parameters: - - output_path (str): Path to save the metadata file. + - root_group (zarr.Group): The root Zarr group. - datasets (list): List of datasets with their metadata. - - args (dict): Additional parameters to include in the metadata. + - output_path (str): Path to save the metadata file. + - metadata_args (dict): Metadata arguments for custom configurations. - mode (str): Mode for metadata handling. Default is 'w'. """ - - if not isinstance(args, dict): - metadata = vars(args) # Convert Namespace to dictionary - multiscales = [{ "version": "0.4", - "name": args.file_name, + "name": "example", "axes": [ {"name": "z", "type": "space", "unit": "micrometer"}, {"name": "y", "type": "space", "unit": "micrometer"}, {"name": "x", "type": "space", "unit": "micrometer"} ], + "datasets": datasets, "type": "gaussian", "metadata": { - "method": "skimage.transform.downscale_local_mean", + "method": "skimage.transform.resize", "version": "0.16.1", - "args": "[true]", - "kwargs": {"anti_aliasing": True, "preserve_range": True}, - "Reco Data": str(args) # Add additional parameters from the params dictionary + "args": [True], + "kwargs": { + "anti_aliasing": True, + "preserve_range": True + } } }] - - multiscales = args2json(multiscales) - + # Update Zarr group attributes if mode == 'w': - root_group.attrs.update({"multiresolution": multiscales}) + root_group.attrs.update({"multiscales": multiscales}) - # Write human-readable JSON to file - metadata_file = os.path.join(output_path, str(args.file_name)[:-2] +'json') + # Save metadata as JSON + metadata_file = os.path.join(output_path, 'multiscales.json') with open(metadata_file, 'w') as f: json.dump({"multiscales": multiscales}, f, indent=4) -def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, compression='blosclz'):#blosclz'): +def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_factors, compression='blosclz'): """ Initialize a multiscale Zarr container with specified levels, dimensions, and compression. Parameters: - output_path (str): Path to the Zarr file. - - base_shape (tuple): Shape of the full dataset at the highest resolution (e.g., (z, y, x)). - - chunks (tuple): Chunk size for the dataset (e.g., (1, y, x)). - - dtype: Data type of the dataset (e.g., np.float32). + - base_shape (tuple): Shape of the full dataset at the highest resolution. + - chunks (tuple): Chunk size for the dataset. + - dtype: Data type of the dataset. - num_levels (int): Number of multiresolution levels. - - compression (str): Compression algorithm (default 'blosclz'). + - scale_factors (list): List of scale factors for each level. + - compression (str): Compression algorithm. Returns: - zarr.Group: The initialized Zarr group containing multiscale datasets. + - list: Dataset metadata for multiscales. """ - store = zarr.DirectoryStore(output_path) compressor = Blosc(cname=compression, clevel=5, shuffle=2) root_group = zarr.group(store=store) + datasets = [] current_shape = base_shape - for level in range(num_levels): + + for level, scale in enumerate(scale_factors): level_name = f"{level}" - log.info(f"Initializing level {level} with shape {current_shape}") - - # Dynamically scale chunks based on the level level_chunks = tuple(max(1, c // (2 ** level)) for c in chunks) - - # Create dataset for the current level + root_group.create_dataset( name=level_name, shape=current_shape, @@ -424,10 +422,17 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, compress dtype=dtype, compressor=compressor ) - + + datasets.append({ + "path": level_name, + "coordinateTransformations": [ + {"type": "scale", "scale": [scale] * 3} + ] + }) + current_shape = tuple(max(1, s // 2) for s in current_shape) - return root_group + return root_group, datasets def write_zarr_chunk(zarr_group, data_chunk, start, end): @@ -491,5 +496,4 @@ def downsample_volume(volume, scale_factor): x // scale_factor, scale_factor ) - # Compute the mean along the downsampling axes return downsampled.mean(axis=(2, 4)) From 78671045afc8e8fef3a0d6782f82bef1ddf3a620 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Wed, 4 Dec 2024 12:17:29 -0600 Subject: [PATCH 16/21] Updated meta saving --- src/tomocupy/config.py | 4 ++++ src/tomocupy/dataio/writer.py | 22 ++++++++++++---------- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index fbbae33..f57efdc 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -463,6 +463,10 @@ def default_parameter(func, param): 'type': str, 'help': "ZARR compression format", 'choices': ['blosclz', 'lz4', 'zstd']}, + 'zarr-chunk': { + 'default': '8,64,64', + 'type': str, + 'help': "ZARR chunk size"}, 'clear-folder': { 'default': 'False', 'type': str, diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 2715049..92b767f 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -221,6 +221,7 @@ def init_output_files(self): self.zarr_output_path = fnameout clean_zarr(self.zarr_output_path) log.info(f'Zarr dataset will be created at {fnameout}') + log.info(f"ZARR chunk structure: {args.zarr_chunk}") params.fnameout = fnameout log.info(f'Output: {fnameout}') @@ -272,17 +273,20 @@ def write_data_chunk(self, rec, st, end, k, shift_index): fid.create_dataset("/exchange/data", data=rec, chunks=(params.nproj, 1, params.n)) elif args.save_format == 'zarr': # Zarr format support - # Define chunk size - #chunks = (1, params.n, params.n) - chunks = (params.nz, 256, 256) + + chunks = [int(c.strip()) for c in args.zarr_chunk.split(',')] if not hasattr(self, 'zarr_array'): + print(rec.shape) + exit(0) shape = (int(params.nz / 2**args.binning), params.n, params.n) # Full dataset shape - + print('initialize') + print(params.nz) max_levels = lambda X, Y: (lambda r: (int(r).bit_length() - 1) if r != 0 else (int(X // Y).bit_length() - 1))(int(X) % int(Y)) levels = min(max_levels(params.nz, end-st),6) + log.info(f"Resolution levels: {levels}") + scale_factors = [float(args.pixel_size) * (i + 1) for i in range(levels)] - self.zarr_array, datasets = initialize_zarr( output_path=self.zarr_output_path, base_shape=shape, @@ -300,7 +304,6 @@ def write_data_chunk(self, rec, st, end, k, shift_index): start=st-shift_index, # Starting index for this chunk along the z-axis end=end-shift_index # Ending index for this chunk along the z-axis ) - def write_data_try(self, rec, cid, id_slice): """Write tiff reconstruction with a given name""" @@ -387,7 +390,7 @@ def fill_zarr_meta(root_group, datasets, output_path, metadata_args, mode='w'): json.dump({"multiscales": multiscales}, f, indent=4) -def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_factors, compression='blosclz'): +def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_factors, compression='None'): """ Initialize a multiscale Zarr container with specified levels, dimensions, and compression. @@ -413,7 +416,7 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_fa for level, scale in enumerate(scale_factors): level_name = f"{level}" - level_chunks = tuple(max(1, c // (2 ** level)) for c in chunks) + level_chunks = tuple(max(1, int(c) // (2 ** level)) for c in chunks) root_group.create_dataset( name=level_name, @@ -434,7 +437,6 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_fa return root_group, datasets - def write_zarr_chunk(zarr_group, data_chunk, start, end): """ Write a chunk of data into the Zarr container at all resolutions. @@ -496,4 +498,4 @@ def downsample_volume(volume, scale_factor): x // scale_factor, scale_factor ) - return downsampled.mean(axis=(2, 4)) + return downsampled.mean(axis=(2, 4)) From abd2dbf431514afc0e7485e4d86b332d8d4cd4b2 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Thu, 5 Dec 2024 18:13:34 -0600 Subject: [PATCH 17/21] bugs over bugs --- src/tomocupy/config.py | 4 ---- src/tomocupy/dataio/reader.py | 1 + src/tomocupy/dataio/writer.py | 36 ++++++++++++++--------------------- 3 files changed, 15 insertions(+), 26 deletions(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index f57efdc..0dc8e70 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -381,10 +381,6 @@ def default_parameter(func, param): 'type': int, 'default': -1, 'help': "End projection"}, - 'nproj-per-chunk': { - 'type': int, - 'default': 8, - 'help': "Number of sinograms per chunk. Use lower numbers with computers with lower GPU memory.", }, 'rotation-axis-auto': { 'default': 'manual', 'type': str, diff --git a/src/tomocupy/dataio/reader.py b/src/tomocupy/dataio/reader.py index e047bd2..86118cf 100644 --- a/src/tomocupy/dataio/reader.py +++ b/src/tomocupy/dataio/reader.py @@ -189,6 +189,7 @@ def init_sizes(self): ltchunk = np.minimum( ncproj, np.int32(nproj-np.arange(ntchunk)*ncproj)) # chunk sizes in proj + tmp = literal_eval(args.nsino) if not isinstance(tmp, list): tmp = [tmp] diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 92b767f..b1a0cff 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -59,6 +59,7 @@ import cupy as cp from pathlib import PosixPath from types import SimpleNamespace +from scipy.ndimage import zoom __author__ = "Viktor Nikitin" @@ -277,13 +278,13 @@ def write_data_chunk(self, rec, st, end, k, shift_index): chunks = [int(c.strip()) for c in args.zarr_chunk.split(',')] if not hasattr(self, 'zarr_array'): - print(rec.shape) - exit(0) + shape = (int(params.nz / 2**args.binning), params.n, params.n) # Full dataset shape print('initialize') - print(params.nz) + max_levels = lambda X, Y: (lambda r: (int(r).bit_length() - 1) if r != 0 else (int(X // Y).bit_length() - 1))(int(X) % int(Y)) - levels = min(max_levels(params.nz, end-st),6) + + levels = min(max_levels(int(params.nz / 2**args.binning), end-st),6) log.info(f"Resolution levels: {levels}") scale_factors = [float(args.pixel_size) * (i + 1) for i in range(levels)] @@ -414,14 +415,14 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_fa datasets = [] current_shape = base_shape - for level, scale in enumerate(scale_factors): + for level, _ in enumerate(scale_factors): level_name = f"{level}" - level_chunks = tuple(max(1, int(c) // (2 ** level)) for c in chunks) + scale = float(args.pixel_size) * (pow(2,level)) root_group.create_dataset( name=level_name, shape=current_shape, - chunks=level_chunks, + chunks=chunks, dtype=dtype, compressor=compressor ) @@ -464,9 +465,6 @@ def write_zarr_chunk(zarr_group, data_chunk, start, end): expected_z_size = level_end - level_start actual_z_size = downsampled_chunk.shape[0] - if actual_z_size != expected_z_size: - downsampled_chunk = downsampled_chunk[:expected_z_size] - # Write the downsampled chunk into the Zarr dataset zarr_array[level_start:level_end, :, :] = downsampled_chunk #log.info(f"Saved chunk to level {level} [{level_start}:{level_end}] with shape {downsampled_chunk.shape}") @@ -474,7 +472,7 @@ def write_zarr_chunk(zarr_group, data_chunk, start, end): def downsample_volume(volume, scale_factor): """ - Downsample a 3D volume by a given scale factor using NumPy. + Downsample a 3D volume by a given scale factor using scipy.ndimage.zoom. Parameters: - volume (numpy array): Input 3D volume (e.g., [z, y, x]). @@ -486,16 +484,10 @@ def downsample_volume(volume, scale_factor): if scale_factor == 1: return volume # No downsampling needed for the highest resolution - # Ensure the input dimensions are divisible by the scale factor - if volume.shape[1] % scale_factor != 0 or volume.shape[2] % scale_factor != 0: - raise ValueError("Volume dimensions must be divisible by the scale factor.") + # Calculate the zoom factors for each axis + zoom_factors = (1 / scale_factor, 1 / scale_factor, 1 / scale_factor) - # Reshape the spatial dimensions for downsampling - z, y, x = volume.shape - downsampled = volume.reshape( - z, - y // scale_factor, scale_factor, - x // scale_factor, scale_factor - ) + # Perform downsampling using interpolation + downsampled = zoom(volume, zoom_factors, order=1) # Use order=1 for bilinear interpolation - return downsampled.mean(axis=(2, 4)) + return downsampled From 02de7d86862128f10862c6c61eb4998c3cbdfa68 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Thu, 5 Dec 2024 18:26:31 -0600 Subject: [PATCH 18/21] clean up --- src/tomocupy/dataio/writer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index b1a0cff..fe33aa6 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -287,7 +287,7 @@ def write_data_chunk(self, rec, st, end, k, shift_index): levels = min(max_levels(int(params.nz / 2**args.binning), end-st),6) log.info(f"Resolution levels: {levels}") - scale_factors = [float(args.pixel_size) * (i + 1) for i in range(levels)] + scale_factors = [float(args.pixel_size) * pow(2,i) for i in range(levels)] self.zarr_array, datasets = initialize_zarr( output_path=self.zarr_output_path, base_shape=shape, From 3a92c83329681be2e00a6b250a80aa41e8badaf0 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Thu, 5 Dec 2024 18:27:01 -0600 Subject: [PATCH 19/21] clean up --- src/tomocupy/dataio/writer.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index fe33aa6..5aa59da 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -415,9 +415,8 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_fa datasets = [] current_shape = base_shape - for level, _ in enumerate(scale_factors): + for level, scale in enumerate(scale_factors): level_name = f"{level}" - scale = float(args.pixel_size) * (pow(2,level)) root_group.create_dataset( name=level_name, From 7f9d08ae1b2bcf69e50a65c33492730fc0297dd4 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Wed, 11 Dec 2024 12:06:57 -0600 Subject: [PATCH 20/21] Restructured/add zarr sum data --- src/tomocupy/config.py | 12 ++ src/tomocupy/dataio/writer.py | 182 ++++++++++++------ src/tomocupy/rec.py | 2 +- .../backproj_lamfourier_parallel.py | 2 +- .../reconstruction/backproj_parallel.py | 4 +- src/tomocupy/utils.py | 40 ++-- 6 files changed, 156 insertions(+), 86 deletions(-) diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py index 0dc8e70..33f75cf 100644 --- a/src/tomocupy/config.py +++ b/src/tomocupy/config.py @@ -381,6 +381,10 @@ def default_parameter(func, param): 'type': int, 'default': -1, 'help': "End projection"}, + 'nproj-per-chunk': { + 'type': int, + 'default': 8, + 'help': "Number of sinograms per chunk. Use lower numbers with computers with lower GPU memory.", }, 'rotation-axis-auto': { 'default': 'manual', 'type': str, @@ -463,6 +467,14 @@ def default_parameter(func, param): 'default': '8,64,64', 'type': str, 'help': "ZARR chunk size"}, + 'large-data': { + 'default': False, + 'type': bool, + 'help': "If Active it computes ldchunk chunks of angular projections"}, + 'ldchunk': { + 'default': 10, + 'type': int, + 'help': "Number of angular chunks for large-data"}, 'clear-folder': { 'default': 'False', 'type': str, diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 5aa59da..1bdbe83 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -41,6 +41,7 @@ from tomocupy import config from tomocupy import logging from tomocupy.global_vars import args, params +from tomocupy.utils import downsampleZarr import numpy as np import h5py import os @@ -59,7 +60,6 @@ import cupy as cp from pathlib import PosixPath from types import SimpleNamespace -from scipy.ndimage import zoom __author__ = "Viktor Nikitin" @@ -220,7 +220,8 @@ def init_output_files(self): if args.save_format == 'zarr': # Zarr format support fnameout += '.zarr' self.zarr_output_path = fnameout - clean_zarr(self.zarr_output_path) + if not args.large_data: + clean_zarr(self.zarr_output_path) log.info(f'Zarr dataset will be created at {fnameout}') log.info(f"ZARR chunk structure: {args.zarr_chunk}") @@ -254,7 +255,7 @@ def write_meta(self, rec_virtual): pass - def write_data_chunk(self, rec, st, end, k, shift_index): + def write_data_chunk(self, rec, st, end, k): """Writing the kth data chunk to hard disk""" if args.save_format == 'tiff': @@ -276,18 +277,17 @@ def write_data_chunk(self, rec, st, end, k, shift_index): elif args.save_format == 'zarr': # Zarr format support chunks = [int(c.strip()) for c in args.zarr_chunk.split(',')] - if not hasattr(self, 'zarr_array'): - shape = (int(params.nz / 2**args.binning), params.n, params.n) # Full dataset shape print('initialize') + max_levels = lambda X, Y: (lambda r: (int(r).bit_length() - 1) if r != 0 else (X.bit_length() - 1))(int(X) % int(Y)) - max_levels = lambda X, Y: (lambda r: (int(r).bit_length() - 1) if r != 0 else (int(X // Y).bit_length() - 1))(int(X) % int(Y)) + #levels = min(max_levels(params.nz, end-st),6) + levels = min(max_levels(int(rec.shape[0]), end-st),6) - levels = min(max_levels(int(params.nz / 2**args.binning), end-st),6) log.info(f"Resolution levels: {levels}") - scale_factors = [float(args.pixel_size) * pow(2,i) for i in range(levels)] + scale_factors = [float(args.pixel_size) * (i + 1) for i in range(levels)] self.zarr_array, datasets = initialize_zarr( output_path=self.zarr_output_path, base_shape=shape, @@ -302,8 +302,8 @@ def write_data_chunk(self, rec, st, end, k, shift_index): write_zarr_chunk( zarr_group=self.zarr_array, # Pre-initialized Zarr container data_chunk=rec[:end - st], # Data chunk to save - start=st-shift_index, # Starting index for this chunk along the z-axis - end=end-shift_index # Ending index for this chunk along the z-axis + start=st-args.start_row, # Starting index for this chunk along the z-axis + end=end-args.start_row # Ending index for this chunk along the z-axis ) def write_data_try(self, rec, cid, id_slice): @@ -371,8 +371,7 @@ def fill_zarr_meta(root_group, datasets, output_path, metadata_args, mode='w'): "datasets": datasets, "type": "gaussian", "metadata": { - "method": "skimage.transform.resize", - "version": "0.16.1", + "method": "scipy.ndimage.zoom", "args": [True], "kwargs": { "anti_aliasing": True, @@ -389,11 +388,58 @@ def fill_zarr_meta(root_group, datasets, output_path, metadata_args, mode='w'): metadata_file = os.path.join(output_path, 'multiscales.json') with open(metadata_file, 'w') as f: json.dump({"multiscales": multiscales}, f, indent=4) + +def write_zarr_chunk(zarr_group, data_chunk, start, end): + """ + Write a chunk of data into the Zarr container at all resolutions, summing with existing data if 'large_data' is active. + Parameters: + - zarr_group (zarr.Group): The initialized Zarr group containing multiscale datasets. + - data_chunk (np.ndarray): The data chunk to write (highest resolution). + - start (int): Start index in the first dimension (z-axis) for the highest resolution. + - end (int): End index in the first dimension (z-axis) for the highest resolution. + - large_data (bool): If True, sum the incoming data chunk with pre-existing data in the Zarr array. + """ + for level in sorted(zarr_group.keys(), key=int): # Process levels in order (0, 1, ...) + zarr_array = zarr_group[level] # Access the dataset for the current level + + # Calculate the downscaling factor for this resolution level + scale_factor = 2 ** int(level) + # Downsample data chunk for this resolution level + if scale_factor > 1: + downsampled_chunk = downsampleZarr(data_chunk, scale_factor) + else: + downsampled_chunk = data_chunk + + # Calculate the start and end indices for the current resolution + level_start = start // scale_factor + level_end = end // scale_factor + + expected_z_size = level_end - level_start + actual_z_size = downsampled_chunk.shape[0] + + if expected_z_size != actual_z_size: + raise ValueError( + f"Mismatch between expected z-size ({expected_z_size}) " + f"and actual z-size ({actual_z_size}) at level {level}." + ) + + # Fetch existing data if 'large_data' is active + if args.large_data: + existing_data = zarr_array[level_start:level_end, :, :] + downsampled_chunk = existing_data + downsampled_chunk + + # Write the (summed or replaced) downsampled chunk into the Zarr dataset + zarr_array[level_start:level_end, :, :] = downsampled_chunk + + # Optionally log the operation + # log.info(f"Saved chunk to level {level} [{level_start}:{level_end}] with shape {downsampled_chunk.shape}") + + def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_factors, compression='None'): """ - Initialize a multiscale Zarr container with specified levels, dimensions, and compression. + Initialize or open a multiscale Zarr container based on the existence of the store. Parameters: - output_path (str): Path to the Zarr file. @@ -403,21 +449,53 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_fa - num_levels (int): Number of multiresolution levels. - scale_factors (list): List of scale factors for each level. - compression (str): Compression algorithm. - + Returns: - - zarr.Group: The initialized Zarr group containing multiscale datasets. + - zarr.Group: The initialized or opened Zarr group containing multiscale datasets. - list: Dataset metadata for multiscales. """ + # Check if the output path (store) already exists + store_exists = os.path.exists(output_path) and os.path.isdir(output_path) + + # Prepare the store reference store = zarr.DirectoryStore(output_path) compressor = Blosc(cname=compression, clevel=5, shuffle=2) + + if store_exists: + return load_zarr(store, output_path, num_levels) + else: + return create_zarr(store, output_path, base_shape, chunks, dtype, num_levels, scale_factors, compressor) + + +def create_zarr(store, output_path, base_shape, chunks, dtype, num_levels, scale_factors, compressor): + """ + Create the entire structure of a new Zarr container. + + Parameters: + - store (zarr.DirectoryStore): Zarr store to initialize. + - output_path (str): Path to the Zarr container for logging purposes. + - base_shape (tuple): Shape of the full dataset at the highest resolution. + - chunks (tuple): Chunk size for the dataset. + - dtype: Data type of the dataset. + - num_levels (int): Number of multiresolution levels. + - scale_factors (list): List of scale factors for each level. + - compressor: Compressor instance for the datasets. + + Returns: + - zarr.Group: The created Zarr group. + - list: Metadata for the datasets. + """ + log.info(f"Creating a new Zarr container at {output_path}") root_group = zarr.group(store=store) - - datasets = [] + current_shape = base_shape + datasets = [] - for level, scale in enumerate(scale_factors): + for level, scale_factor in enumerate(scale_factors): level_name = f"{level}" + scale = float(scale_factor) + # Create the dataset for this resolution level root_group.create_dataset( name=level_name, shape=current_shape, @@ -426,67 +504,51 @@ def initialize_zarr(output_path, base_shape, chunks, dtype, num_levels, scale_fa compressor=compressor ) + # Add metadata datasets.append({ "path": level_name, "coordinateTransformations": [ - {"type": "scale", "scale": [scale] * 3} + {"type": "scale", "scale": [scale] * len(base_shape)} ] }) + # Downscale the shape for the next level current_shape = tuple(max(1, s // 2) for s in current_shape) return root_group, datasets -def write_zarr_chunk(zarr_group, data_chunk, start, end): - """ - Write a chunk of data into the Zarr container at all resolutions. - - Parameters: - - zarr_group (zarr.Group): The initialized Zarr group containing multiscale datasets. - - data_chunk (np.ndarray): The data chunk to write (highest resolution). - - start (int): Start index in the first dimension (z-axis) for the highest resolution. - - end (int): End index in the first dimension (z-axis) for the highest resolution. - """ - for level in sorted(zarr_group.keys(), key=int): # Process levels in order (0, 1, ...) - zarr_array = zarr_group[level] # Access the dataset for the current level - - # Calculate the downscaling factor for this resolution level - scale_factor = 2 ** int(level) - # Downsample data chunk for this resolution level - if scale_factor > 1: - downsampled_chunk = downsample_volume(data_chunk, scale_factor) - else: - downsampled_chunk = data_chunk - - level_start = start // scale_factor - level_end = end // scale_factor - - expected_z_size = level_end - level_start - actual_z_size = downsampled_chunk.shape[0] - # Write the downsampled chunk into the Zarr dataset - zarr_array[level_start:level_end, :, :] = downsampled_chunk - #log.info(f"Saved chunk to level {level} [{level_start}:{level_end}] with shape {downsampled_chunk.shape}") -def downsample_volume(volume, scale_factor): +def load_zarr(store, output_path, num_levels): """ - Downsample a 3D volume by a given scale factor using scipy.ndimage.zoom. + Load an existing Zarr container and return its group and dataset metadata. Parameters: - - volume (numpy array): Input 3D volume (e.g., [z, y, x]). - - scale_factor (int): Factor by which to downsample (e.g., 2 for halving). + - store (zarr.DirectoryStore): The Zarr store to load. + - output_path (str): Path to the Zarr file (for logging purposes). + - num_levels (int): Number of multiresolution levels to validate. Returns: - - numpy array: Downsampled volume. + - zarr.Group: The loaded Zarr group. + - list: Metadata for the datasets. """ - if scale_factor == 1: - return volume # No downsampling needed for the highest resolution + # Open the existing Zarr group in read-write mode + root_group = zarr.open(store=store, mode='r+') + log.info(f"Opened existing Zarr container at {output_path}") - # Calculate the zoom factors for each axis - zoom_factors = (1 / scale_factor, 1 / scale_factor, 1 / scale_factor) + # Gather metadata from the existing structure + datasets = [] + for level in range(num_levels): + level_name = f"{level}" + if level_name not in root_group: + raise ValueError(f"Level {level} not found in the existing Zarr group at {output_path}") - # Perform downsampling using interpolation - downsampled = zoom(volume, zoom_factors, order=1) # Use order=1 for bilinear interpolation + datasets.append({ + "path": level_name, + "coordinateTransformations": [ + {"type": "scale", "scale": None} # Add scale details if needed + ] + }) - return downsampled + return root_group, datasets diff --git a/src/tomocupy/rec.py b/src/tomocupy/rec.py index 4f7090b..aeed5e8 100644 --- a/src/tomocupy/rec.py +++ b/src/tomocupy/rec.py @@ -194,7 +194,7 @@ def recon_all(self): st = ids[k-2]*ncz+args.start_row//2**args.binning end = st+lzchunk[ids[k-2]] self.write_threads[ithread].run( - self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, ids[k-2], args.start_row)) + self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, ids[k-2])) self.stream1.synchronize() self.stream2.synchronize() diff --git a/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py b/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py index b463906..a1649cd 100644 --- a/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py +++ b/src/tomocupy/reconstruction/backproj_lamfourier_parallel.py @@ -290,7 +290,7 @@ def write_parallel(self, u, nthreads=16): st = k*nchunk+params.lamino_start_row end = min((k+1)*nchunk, u.shape[0])+params.lamino_start_row th = Thread(target=self.cl_writer.write_data_chunk, args=( - u[k*nchunk:min((k+1)*nchunk, u.shape[0])], st, end, k, args.start_row)) + u[k*nchunk:min((k+1)*nchunk, u.shape[0])], st, end, k)) mthreads.append(th) th.start() for th in mthreads: diff --git a/src/tomocupy/reconstruction/backproj_parallel.py b/src/tomocupy/reconstruction/backproj_parallel.py index 4355b01..f94249e 100644 --- a/src/tomocupy/reconstruction/backproj_parallel.py +++ b/src/tomocupy/reconstruction/backproj_parallel.py @@ -159,7 +159,7 @@ def recon_sino_proj_parallel(self, data): st = (kr-2)*ncz+args.lamino_start_row//2**args.binning end = st+lrchunk[kr-2] self.write_threads[ithread].run( - self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, kr-2,args.start_row)) + self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, kr-2)) self.stream1.synchronize() self.stream2.synchronize() @@ -364,7 +364,7 @@ def recon_sino_parallel(self, data): st = (k-2)*ncz+args.start_row//2**args.binning end = st+lzchunk[k-2] self.write_threads[ithread].run( - self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, k-2, args.start_row)) + self.cl_writer.write_data_chunk, (rec_pinned[ithread], st, end, k-2)) self.stream1.synchronize() self.stream2.synchronize() diff --git a/src/tomocupy/utils.py b/src/tomocupy/utils.py index 396fb17..64dfb45 100644 --- a/src/tomocupy/utils.py +++ b/src/tomocupy/utils.py @@ -49,7 +49,7 @@ import sys import os import tifffile as tiff -from skimage.transform import resize +from scipy.ndimage import zoom import time from functools import wraps import subprocess @@ -274,33 +274,29 @@ def param_from_dxchange(hdf_file, data_path, attr=None, scalar=True, char_array= return None except KeyError: return None - - -def downsampleZarr(data, scale_factor=2, max_levels=1): + + +def downsampleZarr(volume, scale_factor): """ - Create a multi-level downsampled version of the input data. + Downsample a 3D volume by a given scale factor using scipy.ndimage.zoom. Parameters: - - data (numpy array): The input image data to be downsampled. - - scale_factor (int, optional): Factor by which to downsample the data at each level. Default is 2. - - max_levels (int, optional): Maximum number of downsampled levels to generate. Default is 6. + - volume (numpy array): Input 3D volume (e.g., [z, y, x]). + - scale_factor (int): Factor by which to downsample (e.g., 2 for halving). Returns: - - levels (list of numpy arrays): A list containing the original data and each downsampled level. - - Logs: - - Information about the shape and data type of each downsampled level. + - numpy array: Downsampled volume. """ - current_level = data - levels = [current_level] - for _ in range(max_levels): - new_shape = tuple(max(1, dim // 1) for dim in current_level.shape) - if min(new_shape) <= 1: - break - current_level = current_level#resize(current_level, new_shape, order=0, preserve_range=True, anti_aliasing=True) - print(current_level.shape) - levels.append(current_level) - return levels + if scale_factor == 1: + return volume # No downsampling needed for the highest resolution + + # Calculate the zoom factors for each axis + zoom_factors = (1 / scale_factor, 1 / scale_factor, 1 / scale_factor) + + # Perform downsampling using interpolation + downsampled = zoom(volume, zoom_factors, order=1) # Use order=1 for bilinear interpolation + + return downsampled def clean_zarr(output_path): From 1d2461d9f2c41eaba5a7c6de2de07c2886d7d548 Mon Sep 17 00:00:00 2001 From: mittoalb Date: Sat, 8 Feb 2025 04:58:07 -0600 Subject: [PATCH 21/21] MIP alignment --- src/tomocupy/dataio/reader.py | 1 - src/tomocupy/dataio/writer.py | 9 ++++++--- src/tomocupy/utils.py | 1 + 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/tomocupy/dataio/reader.py b/src/tomocupy/dataio/reader.py index 86118cf..e047bd2 100644 --- a/src/tomocupy/dataio/reader.py +++ b/src/tomocupy/dataio/reader.py @@ -189,7 +189,6 @@ def init_sizes(self): ltchunk = np.minimum( ncproj, np.int32(nproj-np.arange(ntchunk)*ncproj)) # chunk sizes in proj - tmp = literal_eval(args.nsino) if not isinstance(tmp, list): tmp = [tmp] diff --git a/src/tomocupy/dataio/writer.py b/src/tomocupy/dataio/writer.py index 1bdbe83..a6e70a8 100644 --- a/src/tomocupy/dataio/writer.py +++ b/src/tomocupy/dataio/writer.py @@ -494,7 +494,8 @@ def create_zarr(store, output_path, base_shape, chunks, dtype, num_levels, scale for level, scale_factor in enumerate(scale_factors): level_name = f"{level}" scale = float(scale_factor) - + scalef = 2 ** level + # Create the dataset for this resolution level root_group.create_dataset( name=level_name, @@ -508,7 +509,8 @@ def create_zarr(store, output_path, base_shape, chunks, dtype, num_levels, scale datasets.append({ "path": level_name, "coordinateTransformations": [ - {"type": "scale", "scale": [scale] * len(base_shape)} + {"type": "scale", "scale": [scalef] * 3}, + {"type": "translation", "translation": [2**(level-1) - 0.5, 2**(level-1) - 0.5, 2**(level-1) - 0.5]} ] }) @@ -547,7 +549,8 @@ def load_zarr(store, output_path, num_levels): datasets.append({ "path": level_name, "coordinateTransformations": [ - {"type": "scale", "scale": None} # Add scale details if needed + {"type": "scale", "scale": None}, + {"type": "translation", "translation": [2**(level-1) - 0.5, 2**(level-1) - 0.5, 2**(level-1) - 0.5]} ] }) diff --git a/src/tomocupy/utils.py b/src/tomocupy/utils.py index 64dfb45..0cee182 100644 --- a/src/tomocupy/utils.py +++ b/src/tomocupy/utils.py @@ -299,6 +299,7 @@ def downsampleZarr(volume, scale_factor): return downsampled + def clean_zarr(output_path): if os.path.exists(output_path): try: