Skip to content

Input and Output

PtyLab.io

getExampleDataFolder()

Returns a Path with the example data folder.

Source code in PtyLab/io/__init__.py 
4
5
6
7
8
9
def getExampleDataFolder():
    """
    Returns a Path with the example data folder.
    """
    fracPy_folder = Path(__file__).parent.parent.parent
    return fracPy_folder / "example_data"

generateSimulationData

Generate synthetic CPM simulation data and save as simu.hdf5.

This module contains the simulation logic from example_scripts/simulateData.py as a callable function with no side effects (no plotting, no os.chdir).

generate_simu_hdf5(output_path)

Generate a synthetic CPM ptychography dataset and save it as an HDF5 file.

The simulation produces a focused Gaussian probe and a complex spiral object, computes diffraction patterns using Fraunhofer propagation, adds Poisson noise, and writes the result to output_path.

Parameters

output_path : Path Destination file path (e.g. example_data/simu.hdf5). The parent directory must already exist.

Source code in PtyLab/io/generateSimulationData.py 
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def generate_simu_hdf5(output_path: Path) -> None:
    """
    Generate a synthetic CPM ptychography dataset and save it as an HDF5 file.

    The simulation produces a focused Gaussian probe and a complex spiral object,
    computes diffraction patterns using Fraunhofer propagation, adds Poisson noise,
    and writes the result to *output_path*.

    Parameters
    ----------
    output_path : Path
        Destination file path (e.g. ``example_data/simu.hdf5``).
        The parent directory must already exist.
    """
    output_path = Path(output_path)

    # Physical properties
    wavelength = 632.8e-9
    zo = 5e-2
    binningFactor = 1

    # Detector coordinates
    Nd = 2**7
    dxd = 2**11 / Nd * 4.5e-6
    Ld = Nd * dxd

    # Probe coordinates
    dxp = wavelength * zo / Ld
    Np = Nd
    Lp = dxp * Np
    xp = np.arange(-Np // 2, Np // 2) * dxp
    Xp, Yp = np.meshgrid(xp, xp)

    # Object coordinates
    No = 2**10 + 2**9
    dxo = dxp
    Lo = dxo * No
    xo = np.arange(-No // 2, No // 2) * dxo
    Xo, Yo = np.meshgrid(xo, xo)

    # Generate illumination: focused beam via pinhole + lens
    f = 8e-3
    pinhole = circ(Xp, Yp, Lp / 2)
    pinhole = convolve2d(pinhole, gaussian2D(5, 1).astype(np.float32), mode="same")

    probe = aspw(pinhole, 2 * f, wavelength, Lp)[0]

    aperture = circ(Xp, Yp, 3 * Lp / 4)
    aperture = convolve2d(aperture, gaussian2D(5, 3).astype(np.float32), mode="same")
    probe = (
        probe
        * np.exp(-1.0j * 2 * np.pi / wavelength * (Xp**2 + Yp**2) / (2 * f))
        * aperture
    )
    probe = aspw(probe, 2 * f, wavelength, Lp)[0]

    # Generate object: complex spiral pattern
    d = 1e-3
    b = 33
    theta, rho = cart2pol(Xo, Yo)
    t = (1 + np.sign(np.sin(b * theta + 2 * np.pi * (rho / d) ** 2))) / 2
    phaseFun = np.exp(1.0j * (theta + 2 * np.pi * (rho / d) ** 2))
    t = t * circ(Xo, Yo, Lo) * (1 - circ(Xo, Yo, 200 * dxo)) * phaseFun + circ(
        Xo, Yo, 130 * dxo
    )
    obj = convolve2d(t, gaussian2D(5, 3), mode="same")
    object_ = obj * phaseFun

    # Generate scan positions (non-uniform Fermat spiral)
    numPoints = 100
    radius = 150
    p = 1
    R, C = GenerateNonUniformFermat(numPoints, radius=radius, power=p)

    encoder = np.vstack((R * dxo, C * dxo)).T
    positions = np.round(encoder / dxo)
    offset = np.array([50, 20])
    positions = (positions + No // 2 - Np // 2 + offset).astype(int)
    numFrames = len(R)

    # Estimate beam size for entrancePupilDiameter
    beamSize = (
        np.sqrt(np.sum((Xp**2 + Yp**2) * np.abs(probe) ** 2) / np.sum(abs(probe) ** 2))
        * 2.355
    )

    # Compute ptychogram
    ptychogram = np.zeros((numFrames, Nd, Nd))
    for loop in np.arange(numFrames):
        row, col = positions[loop]
        sy = slice(row, row + Np)
        sx = slice(col, col + Np)
        objectPatch = object_[..., sy, sx].copy()
        esw = objectPatch * probe
        ESW = fft2c(esw)
        ptychogram[loop] = abs(ESW) ** 2

    # Simulate Poisson noise
    bitDepth = 14
    maxNumCountsPerDiff = 2**bitDepth
    ptychogram = ptychogram / np.max(ptychogram) * maxNumCountsPerDiff
    noise = np.random.poisson(ptychogram)
    ptychogram += noise
    ptychogram[ptychogram < 0] = 0

    # Save to HDF5
    with h5py.File(output_path, "w") as hf:
        hf.create_dataset("ptychogram", data=ptychogram, dtype="f")
        hf.create_dataset("encoder", data=encoder, dtype="f")
        hf.create_dataset("binningFactor", data=binningFactor, dtype="i")
        hf.create_dataset("dxd", data=(dxd,), dtype="f")
        hf.create_dataset("Nd", data=(Nd,), dtype="i")
        hf.create_dataset("No", data=(No,), dtype="i")
        hf.create_dataset("zo", data=(zo,), dtype="f")
        hf.create_dataset("wavelength", data=(wavelength,), dtype="f")
        hf.create_dataset("entrancePupilDiameter", data=(beamSize,), dtype="f")
        hf.create_dataset("orientation", data=0)

readExample

examplePath(key)

Return the full path to a particular example file. If the key is not found, the key is returned as path. :param key: which dataset to look for. :return:

Source code in PtyLab/io/readExample.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def examplePath(key: str):
    """
    Return the full path to a particular example file. If the key is not found, the key is returned as path.
    :param key: which dataset to look for.
    :return:
    """
    try:
        key = key.split("example:")[1].lower()
        filename = exampleFiles[key]
    except KeyError:
        raise KeyError(
            f"Could not find example {key}. Allowed filenames are: {list(exampleFiles.keys())}"
        )

    full_filename = getExampleDataFolder() / filename

    return full_filename

loadExample(key, *args, **kwargs)

Load an example from the example_data folder. To check which datasets are available, either run listExamples() or look in the folder example_data.

:param key: Key to look for. Will try to look up the filename in readExample.exampleFiles, otherwise it will attempt to find the file in example_data. :param args: will be passed to loadInputData :param kwargs: idem :return:

Source code in PtyLab/io/readExample.py 
52
53
54
55
56
57
58
59
60
61
62
def loadExample(key: str, *args, **kwargs):
    """Load an example from the example_data folder. To check which datasets are available, either run
    listExamples() or look in the folder example_data.

    :param key: Key to look for. Will try to look up the filename in readExample.exampleFiles, otherwise it will attempt to
        find the file in example_data.
    :param args: will be passed to loadInputData
    :param kwargs: idem
    :return:
    """
    return loadInputData(examplePath(key), *args, **kwargs)

readHdf5

scalify(l)

hdf5 file storing (especially when using matlab) can store integers as Numpy arrays of size [1,1]. Convert to scalar if that's the case

Source code in PtyLab/io/readHdf5.py 
15
16
17
18
19
20
21
22
23
24
def scalify(l):
    """
    hdf5 file storing (especially when using matlab) can store integers as
    Numpy arrays of size [1,1]. Convert to scalar if that's the case
    """
    l = l.squeeze()
    try:
        return l.item()
    except ValueError:
        return l

loadInputData(filename, requiredFields, optionalFields)

Load all values from an hdf5 file into a dictionary, but only with the required fields :param filename: the .hdf5 file that has to be loaded. If it's a .mat file it will attempt to load it :param python_order: Weather to read in the files in a way that is common in python, aka for a list of images the first index is the image and not the pixel. :return:

Source code in PtyLab/io/readHdf5.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def loadInputData(filename: Path, requiredFields, optionalFields):
    """
    Load all values from an hdf5 file into a dictionary, but only with the required fields
    :param filename: the .hdf5 file that has to be loaded. If it's a .mat file it will attempt to load it
    :param python_order:
            Weather to read in the files in a way that is common in python, aka for a list of images the first index
             is the image and not the pixel.
    :return:
    """
    filename = Path(filename)
    if not filename.exists():
        raise FileNotFoundError(f"Could not find file {filename}.")
    logger.debug("Loading input data: %s", filename)

    # sanity checks
    if filename.suffix not in allowed_extensions:
        raise NotImplementedError(
            "%s is not a valid extension. Currently, only these extensions are allowed: %s."
            % (filename.suffix, ["   ".join(allowed_extensions)][0])
        )

    # start h5 loading, but check data fields first (defined above)
    dataset = dict()
    try:
        with tables.open_file(str(filename), mode="r") as hdf5File:
            # load the required fields
            for key in requiredFields:
                value = hdf5File.root[key].read()
                dataset[key] = scalify(value)

            # load optional fields, otherwise set to None and compute later
            for key in optionalFields:
                # check if the optional field exists otherwise set to None
                if key in hdf5File.root:
                    value = hdf5File.root[key].read()
                    dataset[key] = scalify(value)
                else:
                    dataset[key] = None

    except Exception as e:
        logger.error("Error reading hdf5 file!")
        raise e
    if "encoder" in dataset:
        print(f"Found encoder with shape {dataset['encoder'].shape}")
        if dataset["encoder"].shape[0] < dataset["encoder"].shape[1]:
            dataset["encoder"] = dataset["encoder"].T
            print("Warning: Automatically changing the shape of encoder. ")
            dataset["encoder"] -= dataset["encoder"].mean(axis=0, keepdims=True)
            print(dataset["encoder"].shape, dataset["encoder"].mean(axis=0))
            # dataset['encoder'] *= -1
    # dirty hack for now

    # upsample

    from PtyLab.utils.utils import fft2c, ifft2c

    # dataset['dxd'] = dataset['dxd'] / 2
    # padwidth = dataset['ptychogram'].shape[-1]//2
    # padwidth = [[0,0], [padwidth, padwidth], [padwidth,padwidth]]
    # dataset['ptychogram'] = abs(ifft2c(np.pad(fft2c(dataset['ptychogram']), padwidth))).astype(np.float32)
    # dataset['ptychogram'] = np.repeat(np.repeat(dataset['ptychogram'], axis=-2, repeats=2), axis=-1, repeats=2)

    return dataset

getOrientation(filename)

If orientation is given, return it. Otherwise, return None

Source code in PtyLab/io/readHdf5.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def getOrientation(filename):
    """
    If orientation is given, return it. Otherwise, return None
    """
    orientation = None
    with h5py.File(str(filename), "r") as archive:
        if "orientation" in archive.keys():
            raw_value = np.array(archive["orientation"]).ravel()[0]
            orientation = scalify(raw_value)
    return orientation

checkDataFields(filename, requiredFields)

Make sure that all the fields in a given .hdf5 file are supported and do some sanity checks.

This is run before loading the file just to make sure that the file is correctly formatted. :param filename: '.hdf5' file with all the necessary attributes. :return: None if correct :raise: KeyError if one of the attributes is missing.

Source code in PtyLab/io/readHdf5.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def checkDataFields(filename, requiredFields):
    """
    Make sure that all the fields in a given .hdf5 file are supported and do some sanity checks.

    This is run before loading the file just to make sure that the file is correctly formatted.
    :param filename: '.hdf5' file with all the necessary attributes.
    :return: None if correct
    :raise: KeyError if one of the attributes is missing.
    """
    with tables.open_file(str(filename), mode="r") as hdf5_file:
        # get a list of nodes
        nodes = hdf5_file.list_nodes("/")
        # get the names of each node which will be the field names stored
        # within the hdf5 file
        fileFields = [node.name for node in nodes]

    # check if all the required fields are within the file
    for k in requiredFields:
        if k not in fileFields:
            raise KeyError("hdf5 file misses key %s" % k)

    return None