<jupyter_start><jupyter_text>Generate Object<jupyter_code>object_true = simulation.ucb()
object_size = yp.shape(object_true)<jupyter_output><empty_output><jupyter_text>Define Settings<jupyter_code># Camera Parameters
camera_pixel_size = 6.5

# Illumination parameters
illumination_device_name = 'quasi-dome'
illumination_wavelength = .53  #um
illumination_na = 0.8

# Objective Parameters
objective_na = .25
objective_mag = 10    #objective mag

# System Parameters
system_mag = 2

# Derived constants
upsampling_factor = math.ceil(illumination_na / objective_na)
pixel_size_measurement = camera_pixel_size / objective_mag
pixel_size_object = pixel_size_measurement / upsampling_factor
dky, dkx = [1 / (sz * pixel_size_object) for sz in object_size]
measurement_size = tuple([int(sz / upsampling_factor) for sz in object_size])<jupyter_output><empty_output><jupyter_text>Generate LED Positions and Pupil Function<jupyter_code># Generate LED Positions
led_position_list_full = ledarray.getPositionsNa(illumination_device_name)
led_board_indicies_full = ledarray.getBoardIndicies(illumination_device_name)
led_position_list = [src for src in led_position_list_full if np.sqrt(src[0] ** 2 + src[1] ** 2) <= illumination_na]
led_board_indicies = [board_index for board_index, src in zip(led_board_indicies_full, led_position_list_full) if np.sqrt(src[0] ** 2 + src[1] ** 2) <= illumination_na]

# Generate Pupil
p = yp.asarray(pupilgen.pupil(measurement_size, pixel_size_measurement, illumination_wavelength, objective_na))<jupyter_output><empty_output><jupyter_text>Apply Distortion to Each Board (Except Center Board)<jupyter_code># Get list of boards used
board_indicies_used = np.unique(led_board_indicies).tolist()

# Remove center board (everything will be registered relative to this board)
board_indicies_used.pop(board_indicies_used.index(0))

# Store distorted led positions
led_position_list_distorted = np.asarray(yp.dcopy(led_position_list))

homography_vectors_true = []

# Set maximum perturbation values
shear_magnitude = 0.1
scale_magnitude = 0.1
rotation_magnitude = 4
translation_magnitude = 0.02
# Apply a random homography to each board index
for board_index in board_indicies_used:
    shear = (yp.rand(2, backend='numpy', dtype='float32') - 0.5) * shear_magnitude
    scale = (yp.rand(2, backend='numpy', dtype='float32') - 0.5) * scale_magnitude
    rotation = (yp.rand(1, backend='numpy', dtype='float32') - 0.5) * rotation_magnitude
    translation = (yp.rand(2, backend='numpy', dtype='float32') - 0.5) * shear_magnitude
    
    # Generate Transformation Matricies
    Sh = transformation.shear(shear)
    Sc = transformation.scale(scale)
    Rt = transformation.rotation(rotation)
    Tr = transformation.translation(translation)
    
    # Compute full homography matrix
    H = Rt.dot(Tr).dot(Sh)
    
    # Store homography vectors
    homography_vectors_true.append(H)
    
    # Apply affine transformation
    led_positions_board_na = led_position_list_distorted[np.asarray(led_board_indicies) == board_index, :]
    led_positions_board_na = np.append(led_positions_board_na, np.ones((len(led_positions_board_na), 1)), axis=1)
    led_position_list_distorted[np.asarray(led_board_indicies) == board_index, :] = H.dot(led_positions_board_na.T).T[:,:2]
    
# Show board Indicies
plt.figure(figsize=(10,4))
plt.subplot(121)
plt.title('Correct LED Positions')
plt.scatter(np.asarray(led_position_list)[:,0], np.asarray(led_position_list)[:,1], c=led_board_indicies)
plt.xlabel('$NA_x$')
plt.ylabel('$NA_y$')
plt.tight_layout()
plt.colorbar()

plt.subplot(122)
plt.title('Distorted LED Positions')
plt.scatter(np.asarray(led_position_list_distorted)[:,0], np.asarray(led_position_list_distorted)[:,1], c=led_board_indicies)
plt.xlabel('$NA_x$')
plt.ylabel('$NA_y$')
plt.tight_layout()
plt.colorbar()<jupyter_output><empty_output><jupyter_text>Generate LED Sequences<jupyter_code># Select LED pattern. Should be a list of list of integers; each internal list has the LED id positions which are illuminated.
single = [[i] for i, x in enumerate(led_position_list)]
multiplexed = [[4 * i + j for j in range(4)] for i in range(8)]
patterns = single<jupyter_output><empty_output><jupyter_text>Generate Forward Models<jupyter_code># Generate operators
I = ops.Intensity(measurement_size)
R = ops.PhaseRamp(object_size)
F = ops.FourierTransform(object_size, label='F_{large}')
Fc = ops.FourierTransform(measurement_size,  label='F_{small}')
P = ops.Diagonalize(p, label='\\tilde P')
CR = ops.Crop(object_size, measurement_size, center=True)

A_list, y_list = [], []
for pattern_index, led_pattern in enumerate(yp.display.progressBar(patterns, name='Forward Models Generated')):
    A_pattern = None
    
    # Loop over all LEDs in this pattern
    for led_id in led_pattern:
        led = led_position_list_distorted[led_id]
        
        # Determine the shift of the LED position in pixels
        pixel_shift = yp.asarray((np.asarray(led).astype(np.complex64) / illumination_wavelength) / np.asarray([dkx, dky]))
        
        # Create phase ramp multiplication operator
        S = ops.Diagonalize(R * pixel_shift, label='S_{' + str(pattern_index) + '}')
        
        # Generate complete forward model
        _A = I * (Fc.H * P * CR * F * S)
        
        # Add this forward operator to the list
        A_pattern = _A if not A_pattern else A_pattern + _A
    
    # Append forward model to list
    A_list.append(A_pattern)
    
    # Generate measurement
    y_list.append(A_pattern * object_true)
    
# Combine these operators into a set if vertically stacked measurements
A = ops.Vstack(A_list)

# Create measurements
y = ops.VecStack(y_list)<jupyter_output><empty_output><jupyter_text>Create Objective Function<jupyter_code># Create object
objective = objectivefunctions.L2(A, y)

# Check gradient
objective.gradient_check(eps=1.0)<jupyter_output><empty_output><jupyter_text>Perform Gradient Descent<jupyter_code>initialization = yp.ones_like(object_true) * yp.mean(object_true)
# Define Initialization
initialization = yp.ones(object_size)

# Solve by accelerated gradient descent
x_opt = ops.solvers.Fista(objective).solve(iteration_count=5, display_type='text', step_size=1e-5,
                                                     use_nesterov_acceleration=True, nesterov_restart_enabled=True, initialization=initialization)

# Show recovered object
plt.figure()
plt.subplot(121)
plt.imshow(np.abs(x_opt), cmap = 'gray')
plt.title('Reconstructed Amp')

plt.subplot(122)
plt.imshow(np.angle(x_opt), cmap = 'gray')
plt.title('Reconstructed Phase')

plt.tight_layout()
plt.show()<jupyter_output>Minimizing function:<jupyter_text>Create Objective Function for LED Positions  Define Homography Matricies$ \begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} = \begin{bmatrix}h1 & h2 & h3 \\ h4 & h5 & h6 \\ h7 & h8 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} $$ \begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} = \begin{bmatrix}x_1 & y_1 & 1 & 0 & 0 & 0 & 0 & 0  \\0 & 0 & 0 & x_1 & y_1 & 1 & 0 & 0  \\0 & 0 & 0 & 0 & 0 & 0 & x_1 & -y_1  \\\end{bmatrix} \begin{bmatrix}h1 \\ h2 \\ h3 \\ h4 \\ h5 \\ h6 \\ h7 \\ h8  \end{bmatrix} $  Affine Transform Matricies$ \begin{bmatrix} x' \\ y' \end{bmatrix} = \begin{bmatrix}h1 & h2 & h3 \\ h4 & h5 & h6 \end{bmatrix} \begin{bmatrix} x \\ y\end{bmatrix} $$ \begin{bmatrix} x' \\ y'\end{bmatrix} = \begin{bmatrix}x_1 & y_1 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & x_1 & y_1 & 1 \\\end{bmatrix} \begin{bmatrix}h1 \\ h2 \\ h3 \\ h4 \\ h5 \\ h6 \end{bmatrix} $  Initialize Homography Vectors<jupyter_code>board_calibration_count = 4

# Create initial affine homography matrix
h_single = yp.asarray([1, 0, 0, 0, 1, 0])

# Tile for all boards and ensure dimensions are correct
h = yp.tile(h_single, [board_calibration_count])
h = yp.expandDims(h, new_dim_count=2)

# Create crop operators for each homography matrix
board_selector_op_list = []
for board_index in range(board_calibration_count):
    board_selector_op_list.append(ops.Crop(yp.shape(h), yp.shape(yp.expandDims(h_single, new_dim_count=2)), crop_start=[board_index * yp.size(h_single)], label='H_{%d}' % board_index))
# Generate operators
I = ops.Intensity(measurement_size)
R = ops.PhaseRamp(object_size)
F = ops.FourierTransform(object_size, label='F_{large}')
Fc = ops.FourierTransform(measurement_size,  label='F_{small}')
P = ops.Diagonalize(p, label='\\tilde P')
CR = ops.Crop(object_size, measurement_size, center=True)
X = ops.Diagonalize(object_true, label='x')

# Generate complete forward model
_A = I * (Fc.H * P * CR * F * X * R)

# Generate measurements from led positions
A_list_illum, y_list_illum = [], []
for pattern_index, led_pattern in enumerate(yp.display.progressBar(patterns, name='Forward Models Generated')):
    A_pattern, y_pattern = None, None
    
    # Loop over all LEDs in this pattern
    for led_id in led_pattern:
        
        # Get led position and board index
        led_position_na = led_position_list[led_id]
        board_index = led_board_indicies[led_id]
        
        if board_index > 0:

            # Determine the shift of the LED position in pixels
            pixel_shift = yp.asarray((np.asarray(led_position_na).astype(np.complex64) / illumination_wavelength) / np.asarray([dkx, dky]))

            # Create homography base matrix for this position
            m = transformation.affineHomographyBlocks(pixel_shift)
            M = ops.MatrixMultiply(m, label='M_{%d}' % led_id)

            # Generate full forward operator
            A = _A * M * board_selector_op_list[board_index - 1]

            # Generate measurement
            _y = A * h

            # Add this forward operator to the list
            A_pattern = A if not A_pattern else A_pattern + A
            y_pattern = _y if not y_pattern else y_pattern + _y
    
    if A_pattern is not None:
        # Append forward model to list
        A_list_illum.append(A_pattern)

        # Generate measurement
        y_list_illum.append(y_pattern)
    
# Combine these operators into a set if vertically stacked measurements
A_illum = ops.Vstack(A_list_illum)

# Create measurements
y_illum = ops.VecStack(y_list_illum)<jupyter_output><empty_output><jupyter_text>Create Objective Function<jupyter_code>objective_illum = objectivefunctions.L2(A_illum, y_illum)

objective_illum.latex()<jupyter_output><empty_output><jupyter_text>Solve for Board Positions<jupyter_code>objective_illum.gradient_check()
# Define Initialization
initialization = h + 1e-9 * yp.rand_like(h)

# Solve by accelerated gradient descent
x_opt = ops.solvers.Fista(objective_illum).solve(iteration_count=20, display_type='text', step_size=5e-12,
                                                     use_nesterov_acceleration=False, nesterov_restart_enabled=False, 
                                                     initialization=initialization)

# Show recovered object
plt.figure()
plt.subplot(121)
plt.imshow(np.abs(x_opt), cmap = 'gray')
plt.title('Reconstructed Amp')

plt.subplot(122)
plt.imshow(np.angle(x_opt), cmap = 'gray')
plt.title('Reconstructed Phase')

plt.tight_layout()
plt.show()
board_calibration_count = 4
# Generate list of true homography matricies for each PCB
h_single = yp.asarray([1, 0, 0, 0, 1, 0])
h = yp.reshape(yp.tile(h_single, [board_calibration_count]), yp.size(h_single), 1)

for board_index in range(board_calibration_count):
    CR = ops.Crop(yp.shape(h), [yp.shape(h_single), 1])
A = A_list_illum[-1]

A.gradient_check()
yp.config.setDefaultDatatype('float32')
backend = 'arrayfire'
x = yp.asbackend(x, backend)
h = yp.asbackend(h, backend)

y1 = yp.matmul(h, x)

_x = np.asarray(x).tolist()
x_new = yp.asarray([[_x[0], _x[1], 1, 0, 0, 0],
                    [0, 0, 0, _x[0], _x[1], 1]], backend=backend)

h_new = yp.vec(h)[:6]

y2 = yp.matmul(x_new, h_new)

print(np.asarray(y1))
print(np.asarray(y2))
x,y = led_position_list[0]

homography_block = [[x, y, 1 , 0, 0, 0, ], [0, 0, 0, -x, -y, -1]
# Compute list of positions used for each board
led_position_list<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11871577899452464, 'fraction_numerical': 0.017836402853824457, 'mean_word_length': 4.338795394154118, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '422131', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4430, 'n_tokens_neox': 4184, 'n_words': 928, 'path': 'notebooks/testing/fpm-simulator-led-correction-Copy1.ipynb', 'repo_name': 'zfphil/pyfpm', 'size': '12054'}