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 dded50b..d21f504 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',
@@ -18,5 +18,6 @@
         'pywavelets',
         'setuptools',  # for pkg_resources at runtime
         'tifffile',
+        'zarr',
     ]
 )
diff --git a/src/tomocupy/config.py b/src/tomocupy/config.py
index 6a4a055..33f75cf 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', 'FourierFilter', '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,
@@ -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'] = {
@@ -395,21 +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']},
-    'clear-folder': {
-        'default': 'False',
-        'type': str,
-        'help': "Clear output folder before reconstruction",
-        'choices': ['True', 'False']},
     'fbp-filter': {
         'default': 'parzen',
         'type': str,
@@ -446,6 +447,41 @@ 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']},    
+    'zarr-chunk': {
+        '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,
+        'help': "Clear output folder before reconstruction",
+        'choices': ['True', 'False']}
+}
+
 
 SECTIONS['beam-hardening'] = {
     'beam-hardening-method': {
@@ -564,9 +600,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')
@@ -773,7 +809,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/reader.py b/src/tomocupy/dataio/reader.py
index ac75e66..e047bd2 100644
--- a/src/tomocupy/dataio/reader.py
+++ b/src/tomocupy/dataio/reader.py
@@ -84,7 +84,7 @@ def init_sizes(self):
 
         # read data sizes and projection angles with a reader
         sizes = self.read_sizes()
-        theta = self.read_theta()
+        theta = self.read_theta(sizes['nproji'])
         nproji = sizes['nproji']
         nzi = sizes['nzi']
         ni = sizes['ni']
@@ -306,12 +306,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..a6e70a8 100644
--- a/src/tomocupy/dataio/writer.py
+++ b/src/tomocupy/dataio/writer.py
@@ -41,12 +41,27 @@
 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
 import sys
 import tifffile
 
+#Zarr writer
+import shutil
+import zarr
+import json
+from numcodecs import Blosc
+import threading
+import subprocess
+import time
+
+import cupy as cp
+from pathlib import PosixPath
+from types import SimpleNamespace
+
+
 __author__ = "Viktor Nikitin"
 __copyright__ = "Copyright (c) 2022, UChicago Argonne, LLC."
 __docformat__ = 'restructuredtext en'
@@ -202,9 +217,17 @@ 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
+            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}")
+              
         params.fnameout = fnameout
         log.info(f'Output: {fnameout}')
+        
 
     def write_meta(self, rec_virtual):
 
@@ -231,6 +254,7 @@ def write_meta(self, rec_virtual):
             log.error('write_meta() error: Skip copying meta')
             pass
 
+
     def write_data_chunk(self, rec, st, end, k):
         """Writing the kth data chunk to hard disk"""
 
@@ -250,9 +274,284 @@ 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':  # 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))
+
+                #levels = min(max_levels(params.nz, end-st),6)
+                levels = min(max_levels(int(rec.shape[0]), 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,
+                    chunks=chunks,
+                    dtype=params.dtype,
+                    num_levels=levels,
+                    scale_factors=scale_factors,
+                    compression=args.zarr_compression
+                )
+                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
+                data_chunk=rec[:end - st],  # Data chunk to save
+                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):
         """Write tiff reconstruction with a given name"""
 
         tifffile.imwrite(
             f'{params.fnameout}_slice{id_slice:04d}_center{cid:05.2f}.tiff', rec)
+                        
+            
+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 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 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, datasets, output_path, metadata_args, mode='w'):
+    """
+    Fill metadata for the Zarr multiscale datasets and include additional parameters.
+
+    Parameters:
+    - root_group (zarr.Group): The root Zarr group.
+    - datasets (list): List of datasets with their 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'.
+    """
+    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": "scipy.ndimage.zoom",
+            "args": [True],
+            "kwargs": {
+                "anti_aliasing": True,
+                "preserve_range": True
+            }
+        }
+    }]
+
+    # Update Zarr group attributes
+    if mode == 'w':
+        root_group.attrs.update({"multiscales": multiscales})
+
+        # 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 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 or open a multiscale Zarr container based on the existence of the store.
+
+    Parameters:
+    - output_path (str): Path to the Zarr file.
+    - 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.
+    - compression (str): Compression algorithm.
+
+    Returns:
+    - 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)
+
+    current_shape = base_shape
+    datasets = []
+
+    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,
+            shape=current_shape,
+            chunks=chunks,
+            dtype=dtype,
+            compressor=compressor
+        )
+
+        # Add metadata
+        datasets.append({
+            "path": level_name,
+            "coordinateTransformations": [
+                {"type": "scale", "scale": [scalef] * 3},
+                {"type": "translation", "translation": [2**(level-1) - 0.5, 2**(level-1) - 0.5, 2**(level-1) - 0.5]}
+            ]
+        })
+
+        # Downscale the shape for the next level
+        current_shape = tuple(max(1, s // 2) for s in current_shape)
+
+    return root_group, datasets
+
+
+
+
+def load_zarr(store, output_path, num_levels):
+    """
+    Load an existing Zarr container and return its group and dataset metadata.
+
+    Parameters:
+    - 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:
+    - zarr.Group: The loaded Zarr group.
+    - list: Metadata for the datasets.
+    """
+    # 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}")
+
+    # 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}")
+
+        datasets.append({
+            "path": level_name,
+            "coordinateTransformations": [
+                {"type": "scale", "scale": None},
+                {"type": "translation", "translation": [2**(level-1) - 0.5, 2**(level-1) - 0.5, 2**(level-1) - 0.5]}
+            ]
+        })
+
+    return root_group, datasets
diff --git a/src/tomocupy/processing/proc_functions.py b/src/tomocupy/processing/proc_functions.py
index 1807a5d..5a339a0 100644
--- a/src/tomocupy/processing/proc_functions.py
+++ b/src/tomocupy/processing/proc_functions.py
@@ -171,6 +171,14 @@ 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.retrieve_phase_method == 'farago':        
+                data[:] = retrieve_phase.farago_filter(
+                data,  args.pixel_size*1e-4, args.propagation_distance/10, args.energy,
+                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 b789910..d52911b 100644
--- a/src/tomocupy/processing/retrieve_phase.py
+++ b/src/tomocupy/processing/retrieve_phase.py
@@ -107,6 +107,48 @@ 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
+    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 +209,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):
     """
@@ -258,3 +302,113 @@ 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
+    #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:
+            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
+            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])
+            data[:, m, :] = mat_dec
+
+    if apply_log:
+        data = cp.exp(-data)
+
+    return data
+    
+    
diff --git a/src/tomocupy/utils.py b/src/tomocupy/utils.py
index e7c2afd..0cee182 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 scipy.ndimage import zoom
+import time
+from functools import wraps
+import subprocess
+
+
+
 from tomocupy import logging
 log = logging.getLogger(__name__)
 
@@ -265,3 +274,41 @@ def param_from_dxchange(hdf_file, data_path, attr=None, scalar=True, char_array=
                 return None
         except KeyError:
             return None
+
+
+def downsampleZarr(volume, scale_factor):
+    """
+    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]).
+    - 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
+
+    # 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):
+    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 1d282c7..db43aa5 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,    
+    # '{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,9 +85,33 @@ def test_full_recon(self):
                         f'data_rec/{file_name}_rec/recon_{k:05}.tiff'))
                 except:
                     pass
-            self.assertAlmostEqual(ssum, cmd[1], places=0)
+            #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)
+
+            # 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)