PyRate_FieldSimulation.py

import math
import numpy as np


def generate_line_spot_map_order(field_details):
    spot_map_order = np.empty((0, 2), int)

    spots_number_vertical = field_details[1]
    spots_number_horizontal = field_details[2]

    centre_vertical = math.floor((spots_number_vertical-1) / 2)
    move = np.array([1,0])
    position = np.array([0,centre_vertical])

    for x in range(spots_number_horizontal):
        spot_map_order = np.concatenate((spot_map_order, [position]), axis=0)
        position += move

    return spot_map_order


def generate_spiral_spot_map_order(field_details, direction):
    spot_map_order = np.empty((0,2),int)
    spots_number_vertical = field_details[1]
    spots_number_horizontal = field_details[2]
    total_spot_number = spots_number_horizontal*spots_number_vertical

    # (di, dj) is a vector - direction in which we move right now
    move = np.array([-1,0])

    # current position (i, j), start top left
    position = np.array([spots_number_horizontal-1,0])
    spot_map_order = np.concatenate((spot_map_order, [position]), axis=0)

    segment_length_h = spots_number_horizontal-1
    segment_length_v = spots_number_vertical-1

    current_segment = 0

    current_spot_number = 1

    # initial line
    for x in range(segment_length_h):
        position += move
        spot_map_order = np.concatenate((spot_map_order, [position]), axis=0)
        current_spot_number += 1

    #change direction
    previous_di = move[0]
    previous_dj = move[1]
    move = [0,1]

    while current_spot_number != total_spot_number:
        position += move
        spot_map_order = np.concatenate((spot_map_order, [position]), axis=0)
        current_segment += 1

        if current_segment == segment_length_v and move[1] != 0:

            previous_dj = move[1]

            if previous_di == -1:
                move = [1, 0]
            if previous_di == 1:
                move = [-1, 0]

            segment_length_v -= 1
            current_segment = 0

        if current_segment == segment_length_h and move[0] != 0:

            previous_di = move[0]

            if previous_dj == -1:
                move = [0, 1]
            if previous_dj == 1:
                move = [0, -1]

            segment_length_h -= 1
            current_segment = 0

        current_spot_number += 1

    if direction == "out":
        spot_map_order = np.flip(spot_map_order)

    return spot_map_order


def generate_raster_dose_map(field_details):

    spot_map_order = np.empty((0,2),int)
    spots_number_vertical = field_details[1]
    spots_number_horizontal = field_details[2]

    left_to_right = True

    for y in range(spots_number_vertical):
        for x in range(spots_number_horizontal):
            if not left_to_right:
                x = (spots_number_horizontal - 1) - x

            position = np.array([x, y])
            spot_map_order = np.concatenate((spot_map_order, [position]), axis=0)

        if left_to_right:
            left_to_right = False
        elif not left_to_right:
            left_to_right = True

    return spot_map_order


def generate_circular_raster(field_details):
    raster_spot_map_order = generate_raster_dose_map(field_details)
    spot_map_order = []

    max_x = 0
    max_y = 0
    for x, y in raster_spot_map_order:
        if x > max_x:
            max_x = x
        if y > max_y:
            max_y = y
    max_x += 1
    max_y += 1

    centre_x = max_x / 2
    centre_y = max_y / 2
    radius = max(max_x, max_y) / 2

    for x, y in raster_spot_map_order:
        vector = [x - centre_x, y - centre_y]
        distance = pow(pow(vector[0], 2) + pow(vector[1], 2), 0.5)
        if distance <= radius:
            spot_map_order.append([x, y])

    return spot_map_order


def generate_spot_map(spot_map_order, spot_details):
    spot_position_store = []
    spot_spacing_pixels = spot_details[2]
    spot_pixels_max = spot_details[1]
    for x, y in spot_map_order:
        spot_position_y = int(y * spot_spacing_pixels) + spot_pixels_max
        spot_position_x = int(x * spot_spacing_pixels) + spot_pixels_max
        spot_position = (spot_position_x, spot_position_y)
        spot_position_store.append(spot_position)

    return spot_position_store


def create_dose_map(target_details, spot_details, overscan, resolution):
    # Work out xy size of target
    # Work out pixel size of target using resolution
    # Create dose_map array using pixel size of target
    # Pad dose_map with spot pixel maximum. This is the size of the spot map needed to ensure a
    #   homogeneous dose across the target.
    # Pad dose_map again with spot pixel maximum so we capture the dose falloff

    target_horizontal_size = target_details[0]
    target_vertical_size = target_details[1]
    spot_spacing_pixels = spot_details[2]
    spot_pixels_max = spot_details[1]

    dose_map_horizontal_target = math.ceil(target_horizontal_size / resolution)
    dose_map_vertical_target = math.ceil(target_vertical_size / resolution)

    dose_map_horizontal_dose_area = dose_map_horizontal_target + 2*overscan[0]  # pad for homogeneous dose
    dose_map_vertical_dose_area = dose_map_vertical_target +  2*overscan[1]  # pad for homogeneous dose


    dose_map_horizontal_display_area = dose_map_horizontal_dose_area + 2*spot_pixels_max  # pad for plotting falloff
    dose_map_vertical_display_area = dose_map_vertical_dose_area + 2*spot_pixels_max  # pad for plotting falloff
    column = dose_map_horizontal_display_area
    row = dose_map_vertical_display_area

    dose_map = np.zeros((row, column))

    # Work out how many pixels between spots by using spot_spacing and resolution
    # Work out how many x & y spot are required using pixels between spots and dose map dose area
    spots_number_horizontal = math.ceil(dose_map_horizontal_dose_area / spot_spacing_pixels) + 1
    spots_number_vertical = math.ceil(dose_map_vertical_dose_area / spot_spacing_pixels) + 1

    details = [dose_map_horizontal_target,
               dose_map_vertical_target,
               dose_map_horizontal_dose_area,
               dose_map_vertical_dose_area,
               spots_number_horizontal,
               spots_number_vertical,
               column,
               row]

    return dose_map, spots_number_vertical, spots_number_horizontal, details


def simulate_spot_with_kernel(dose_map, kernel, spot_position, spot_pixels_max):

    #  Work out the i & j range of dose_map to iterate over by using spot position
    #   and the pixel extent of the spot
    min_i = spot_position[0] - spot_pixels_max  # x
    max_i = spot_position[0] + spot_pixels_max  # x
    min_j = spot_position[1] - spot_pixels_max  # y
    max_j = spot_position[1] + spot_pixels_max  # y


    if min_i < 0:
        min_i = 0
    if min_j < 0:
        min_j = 0
    if max_i > dose_map.shape[1]:
        max_i = dose_map.shape[1]
    if max_j > dose_map.shape[0]:
        max_j = dose_map.shape[0]

    for j in range(min_j, max_j):
        for i in range(min_i, max_i):
            # Generate gaussian values
            dose_map[j][i] += kernel[j-min_j][i-min_i]

    return dose_map


def simulate_field_with_defined_spots(field_details, kernel, spot_map, point_of_interest):
    # Break composite variables into parts
    dose_map = field_details[0]
    spot_pixels_max = int(len(kernel)/2)

    for spot_position_x, spot_position_y in spot_map:

        spot_position = (spot_position_x, spot_position_y)

        # ------------- start looping over dose_map -------------
        dose_map = simulate_spot_with_kernel(dose_map, kernel, spot_position, spot_pixels_max)
        # ------------- done looping over dose_map -------------

    # select midpoint of field
    dose = dose_map[point_of_interest[0][0], point_of_interest[0][1]]

    return dose


def simulate_field_and_time_with_defined_spots(field_details, resolution, spot_delivery_time, scan_speed_horizontal,
                                               scan_speed_vertical, kernel, spot_map, point_of_interest):

    time = 0.0
    data = []
    history = []
    spot_number = 0
    previous_spot_x = 0
    previous_spot_y = 0

    dose_map = field_details[0]
    spot_pixels_max = int(len(kernel)/2)

    for x, y in spot_map:
        diff_x = abs(previous_spot_x - x) * resolution
        diff_y = abs(previous_spot_y - y) * resolution

        time_x = diff_x/scan_speed_horizontal
        time_y = diff_y/scan_speed_vertical
        time += max(time_x, time_y)

        spot_position = (x, y)

        spot_number += 1

        previous_spot_x = x
        previous_spot_y = y


        # ------------- start looping over dose_map -------------
        dose_map = simulate_spot_with_kernel(dose_map, kernel, spot_position, spot_pixels_max)
        # ------------- done looping over dose_map -------------

        time += spot_delivery_time

        line = [time]

        line.append(dose_map[point_of_interest[0][0], point_of_interest[0][1]])
        line.append(dose_map[point_of_interest[1][0], point_of_interest[1][1]])
        line.append(dose_map[point_of_interest[2][0], point_of_interest[2][1]])
        line.append(dose_map[point_of_interest[3][0], point_of_interest[3][1]])


        data.append(line)
        dose_map_copy = np.copy(dose_map)
        history.append([time, dose_map_copy])

    print(f'Total irradiation time: {time} s\n')

    return data, history