import time
import cv2
import numpy as np
import importlib_resources
import reikna.cluda as cluda
import reikna.cluda.dtypes as dtypes
from reikna.core import Transformation, Parameter, Annotation, Type
from reikna.cluda import functions, dtypes
from reikna.fft import FFT, FFTShift
from reikna.cluda.tempalloc import ZeroOffsetManager
import dpivsoft.DPIV as DPIV # Python PIV implementation
from dpivsoft.Classes import Parameters
from dpivsoft.Classes import grid
from dpivsoft.Classes import GPU
[docs]
def compile_Kernels(thr):
"""
Compile all kernels needed for GPU calculation. Only needs to be
called once per device.
Kernel sources live in the package folder "GPU_Kernels". SubMean
and find_peak are compiled later, in initialization(), because
their local-memory sizes depend on the box geometry.
Parameters
----------
thr : reikna.cluda Thread
Device thread returned by select_Platform.
"""
# Package installation path
path = importlib_resources.files("dpivsoft")
# Split image Kernel
with open(path/"GPU_Kernels/Slice.cl", "r") as file:
program = thr.compile(file.read())
GPU.Slice = program.Slice
# Initialize u_index_1 and v_index_1
with open(path/"GPU_Kernels/ini_index.cl", "r") as file:
program = thr.compile(file.read())
GPU.ini_index = program.ini_index
# SubMean is compiled in initialization() because it needs box sizes
# to allocate local memory at compile time.
with open(path/"GPU_Kernels/Normalize_Img.cl", "r") as file:
program = thr.compile(file.read())
GPU.Normalize = program.Normalize
# Multiplication Kernel
with open(path/"GPU_Kernels/multiply_them.cl", "r") as file:
program = thr.compile(file.read(),
render_kwds=dict(ctype1=dtypes.ctype(np.complex64),
ctype2=dtypes.ctype(np.complex64),
mul=functions.mul(np.complex64, np.complex64)))
GPU.multiply_them = program.multiply_them
# Apply mask Kernel
with open(path/"GPU_Kernels/multiply_them.cl", "r") as file:
program = thr.compile(file.read(),
render_kwds=dict(ctype1=dtypes.ctype(np.float32),
ctype2=dtypes.ctype(bool),
mul=functions.mul(np.float32, bool)))
GPU.masking = program.multiply_them
# find_peak is compiled in initialization() — box/window geometry must be
# known at compile time for local memory sizing (Mako template).
# Interpolation Kernel
with open(path/"GPU_Kernels/interpolation.cl", "r") as file:
program = thr.compile(file.read())
GPU.interpolate = program.Interpolation
# Jacobian Kernel
with open(path/"GPU_Kernels/Jacobian.cl", "r") as file:
program = thr.compile(file.read())
GPU.jacobian = program.Jacobian
# Box blur filter
with open(path/"GPU_Kernels/box_blur.cl", "r") as file:
program = thr.compile(file.read())
GPU.box_blur = program.box_blur
# Deform image kernel
with open(path/"GPU_Kernels/deform_image.cl", "r") as file:
program = thr.compile(file.read())
GPU.deform_image = program.Deform_image
# Median Filter
with open(path/"GPU_Kernels/median_filter.cl", "r") as file:
program = thr.compile(file.read())
GPU.Median_Filter = program.Median_Filter
# Weighting function
with open(path/"GPU_Kernels/Weighting.cl", "r") as file:
program = thr.compile(file.read())
GPU.Weighting = program.Weighting
# Gaussian blur
with open(path/"GPU_Kernels/gaussian_filter.cl", "r") as file:
program = thr.compile(file.read())
GPU.gaussian_filter = program.gaussian_filter
# In-place conjugate kernel
with open(path/"GPU_Kernels/conjugate_inplace.cl", "r") as file:
program = thr.compile(file.read())
GPU.conjugate = program.conjugate
# directCorrelation
with open(path/"GPU_Kernels/directCorrelation.cl", "r") as file:
program = thr.compile(file.read())
GPU.directCorrelation = program.directCorrelation
with open(path/"GPU_Kernels/find_peak_direct.cl", "r") as file:
program = thr.compile(file.read())
GPU.find_peak_direct = program.find_peak_direct
thr.synchronize()
return 0
[docs]
def initialization(width, height, thr):
"""
Allocate all GPU buffers and compile the geometry-dependent kernels
(SubMean, find_peak). Call once after setting Parameters and before
processing(); call again if the image size or the PIV parameters
change.
Parameters
----------
width, height : int
Image size in pixels.
thr : reikna.cluda Thread
Device thread returned by select_Platform.
"""
# Obtain PIV mesh
grid.generate_mesh(width, height)
grid.pixels = width * height
Parameters.Data.height = height
Parameters.Data.width = width
# Total number of boxes for global size
N_boxes_1 = Parameters.no_boxes_1_x * Parameters.no_boxes_1_y
N_boxes_2 = Parameters.no_boxes_2_x * Parameters.no_boxes_2_y
# subImg examples to compile fft kernel
subImg1 = np.zeros((N_boxes_1, Parameters.box_size_1_y,
Parameters.box_size_1_x), dtype=np.complex64)
subImg2 = np.zeros((N_boxes_2, Parameters.box_size_2_y,
Parameters.box_size_2_x), dtype=np.complex64)
# Array of parameters data
peak_ratio = int(Parameters.peak_ratio * 1000) # trick to use as int
data1 = np.array((width, height,
Parameters.box_size_1_x, Parameters.box_size_1_y,
Parameters.no_boxes_1_x, Parameters.no_boxes_1_y,
Parameters.window_1_x, Parameters.window_1_y,
peak_ratio, Parameters.gaussian_size)).astype(np.int32)
data2 = np.array((width, height,
Parameters.box_size_2_x, Parameters.box_size_2_y,
Parameters.no_boxes_2_x, Parameters.no_boxes_2_y,
Parameters.window_2_x, Parameters.window_2_y,
peak_ratio, Parameters.gaussian_size)).astype(np.int32)
# Send mesh to gpu (only done once)
GPU.box_origin_x_1 = thr.to_device(grid.box_origin_x_1)
GPU.box_origin_y_1 = thr.to_device(grid.box_origin_y_1)
GPU.box_origin_x_2 = thr.to_device(grid.box_origin_x_2)
GPU.box_origin_y_2 = thr.to_device(grid.box_origin_y_2)
GPU.x1 = thr.to_device(grid.x_1)
GPU.y1 = thr.to_device(grid.y_1)
GPU.x2 = thr.to_device(grid.x_2)
GPU.y2 = thr.to_device(grid.y_2)
if Parameters.mask:
GPU.mask_1 = thr.to_device(grid.mask_1)
GPU.mask_2 = thr.to_device(grid.mask_2)
if Parameters.direct_calc:
GPU.box_origin_x_d = thr.to_device(grid.box_origin_x_d)
GPU.box_origin_y_d = thr.to_device(grid.box_origin_y_d)
box_size_x_d = round(Parameters.window_1_x / 2) * 2 + 1
box_size_y_d = round(Parameters.window_1_y / 2) * 2 + 1
directCorr = np.zeros((N_boxes_1, box_size_y_d,
box_size_x_d), dtype=np.float32)
GPU.size_direct = np.prod(directCorr.shape)
GPU.directCorr = thr.to_device(directCorr)
# Send PIV parameters to gpu (only done once)
GPU.data1 = thr.to_device(data1)
GPU.data2 = thr.to_device(data2)
GPU.median_limit = thr.to_device(np.float32(Parameters.median_limit))
# Temp manager to reduce use of memory
temp_manager = ZeroOffsetManager(
thr, pack_on_alloc=True, pack_on_free=False)
# Initialize all GPU variables for first iteration (only done once)
GPU.subImg1_1 = temp_manager.array([N_boxes_1,
Parameters.box_size_1_y, Parameters.box_size_1_x],
np.complex64)
GPU.subImg1_2 = temp_manager.array([N_boxes_1,
Parameters.box_size_1_y, Parameters.box_size_1_x],
np.complex64, dependencies=[GPU.subImg1_1])
GPU.subMean1_1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2])
GPU.subMean1_2 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1])
GPU.u1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2])
GPU.v1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1])
GPU.u1_f = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1])
GPU.v1_f = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1, GPU.u1_f])
GPU.du_dx_1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1, GPU.u1_f,
GPU.v1_f])
GPU.du_dy_1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1, GPU.u1_f,
GPU.v1_f, GPU.du_dx_1])
GPU.dv_dx_1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1, GPU.u1_f,
GPU.v1_f, GPU.du_dx_1, GPU.du_dy_1])
GPU.dv_dy_1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1, GPU.u1_f,
GPU.v1_f, GPU.du_dx_1, GPU.du_dy_1, GPU.dv_dx_1])
GPU.temp_dx_1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1, GPU.u1_f,
GPU.v1_f, GPU.du_dx_1, GPU.du_dy_1, GPU.dv_dx_1,
GPU.dv_dy_1])
GPU.temp_dy_1 = temp_manager.array(
[Parameters.no_boxes_1_y, Parameters.no_boxes_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1, GPU.u1_f,
GPU.v1_f, GPU.du_dx_1, GPU.du_dy_1, GPU.dv_dx_1,
GPU.dv_dy_1, GPU.temp_dx_1])
GPU.u_index_1 = temp_manager.array([N_boxes_1,
Parameters.box_size_1_y, Parameters.box_size_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1,
GPU.du_dx_1, GPU.du_dy_1, GPU.dv_dx_1,
GPU.dv_dy_1, GPU.temp_dx_1, GPU.temp_dy_1])
GPU.v_index_1 = temp_manager.array([N_boxes_1,
Parameters.box_size_1_y, Parameters.box_size_1_x],
np.float32, dependencies=[GPU.subImg1_1, GPU.subImg1_2,
GPU.subMean1_1, GPU.subMean1_2, GPU.u1, GPU.v1,
GPU.du_dx_1, GPU.du_dy_1, GPU.dv_dx_1,
GPU.dv_dy_1, GPU.temp_dx_1, GPU.temp_dy_1,
GPU.u_index_1])
# Create temp images for gaussian filter if needed
if (Parameters.gaussian_size):
kernel = DPIV.gaussian_kernel(Parameters.gaussian_size)
GPU.kernel = thr.to_device(np.float32(kernel))
# Images filtered to be used only on first iteration
GPU.img1_g = temp_manager.array([height, width],
np.float32, dependencies=[GPU.subImg1_1,
GPU.subImg1_2, GPU.u_index_1, GPU.v_index_1])
GPU.img2_g = temp_manager.array([height, width],
np.float32, dependencies=[GPU.subImg1_1,
GPU.subImg1_2, GPU.u_index_1, GPU.v_index_1,
GPU.img1_g])
# Initialize all GPU variables for second iteration (only done once)
GPU.subImg2_1 = temp_manager.array([N_boxes_2,
Parameters.box_size_2_y, Parameters.box_size_2_x],
np.complex64)
GPU.subImg2_2 = temp_manager.array([N_boxes_2,
Parameters.box_size_2_y, Parameters.box_size_2_x],
np.complex64, dependencies=[GPU.subImg2_1])
GPU.subMean2_1 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2])
GPU.subMean2_2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1])
GPU.u2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2])
GPU.v2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2])
GPU.u2_f = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2, GPU.u1])
GPU.v2_f = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2, GPU.u2_f, GPU.v1])
GPU.du_dx_2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2])
GPU.du_dy_2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2, GPU.du_dx_2])
GPU.dv_dx_2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2,
GPU.du_dx_2, GPU.du_dy_2])
GPU.dv_dy_2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2,
GPU.du_dx_2, GPU.du_dy_2, GPU.dv_dx_2])
GPU.temp_dx_2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2,
GPU.du_dx_2, GPU.du_dy_2, GPU.dv_dx_2,
GPU.dv_dy_2])
GPU.temp_dy_2 = temp_manager.array(
[Parameters.no_boxes_2_y, Parameters.no_boxes_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2,
GPU.du_dx_2, GPU.du_dy_2, GPU.dv_dx_2,
GPU.dv_dy_2, GPU.temp_dx_2])
GPU.u_index_2 = temp_manager.array([N_boxes_2,
Parameters.box_size_2_y, Parameters.box_size_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2,
GPU.du_dx_2, GPU.du_dy_2, GPU.dv_dx_2,
GPU.dv_dy_2, GPU.temp_dx_2, GPU.temp_dy_2])
GPU.v_index_2 = temp_manager.array([N_boxes_2,
Parameters.box_size_2_y, Parameters.box_size_2_x],
np.float32, dependencies=[GPU.subImg2_1, GPU.subImg2_2,
GPU.subMean2_1, GPU.subMean2_2, GPU.u2, GPU.v2,
GPU.du_dx_2, GPU.du_dy_2, GPU.dv_dx_2,
GPU.dv_dy_2, GPU.temp_dx_2, GPU.temp_dy_2,
GPU.u_index_2])
# Load mask if any
if Parameters.mask:
GPU.mask = thr.to_device(Parameters.Data.mask)
else:
temp = np.zeros(2)
GPU.mask = thr.empty_like(temp)
# Pre-allocate image buffers (reused every frame to avoid per-frame allocation)
GPU.img1 = thr.array((height, width), np.float32)
GPU.img2 = thr.array((height, width), np.float32)
# Compile SubMean here (not in compile_Kernels) because local memory size
# must be a compile-time constant and depends on box geometry.
max_wg = thr.device_params.max_work_group_size
submean_path = importlib_resources.files("dpivsoft") / "GPU_Kernels/SubMean.cl"
def _next_pow2(n):
"""Return the largest power of two <= n."""
p = 1
while p * 2 <= n:
p *= 2
return p
# rusticl/radeonsi reports a per-kernel CL_KERNEL_WORK_GROUP_SIZE
# (register/local-memory bound) below the device-wide
# CL_DEVICE_MAX_WORK_GROUP_SIZE. Honouring only the device limit makes
# clEnqueueNDRangeKernel raise INVALID_WORK_GROUP_SIZE on those drivers,
# so the local size must be clamped to this per-kernel value. reikna
# already queries it and exposes it as kernel.max_work_group_size (see
# reikna/cluda/ocl.py), so no private-attribute access is needed.
def _compile_submean(box_size_x, box_size_y):
"""Compile the SubMean reduction kernel for a given window size,
clamping the local size to the kernel's real CL_KERNEL_WORK_GROUP_SIZE
(recompiling if the driver's per-kernel limit is lower than requested).
Returns (kernel, local_size)."""
box_size = box_size_x * box_size_y
# box_size is a power of two, so a power-of-two local size always
# divides it evenly (required by the reduction) and is a valid
# workgroup size for the binary-tree reduction.
cap = max_wg
while True:
local_size = _next_pow2(min(box_size, cap))
items_per_thread = box_size // local_size
with open(submean_path, "r") as f:
prog = thr.compile(f.read(), render_kwds=dict(
box_size=box_size,
local_size=local_size,
items_per_thread=items_per_thread))
kernel = prog.SubMean
limit = kernel.max_work_group_size
if local_size <= limit:
return kernel, local_size
cap = limit # recompile clamped to the kernel's real limit
GPU.SubMean1, GPU.SubMean1_local = _compile_submean(
Parameters.box_size_1_x, Parameters.box_size_1_y)
GPU.SubMean2, GPU.SubMean2_local = _compile_submean(
Parameters.box_size_2_x, Parameters.box_size_2_y)
# Compile find_peak here for the same reason as SubMean: local memory
# arrays must have compile-time sizes (Mako template).
find_peak_path = importlib_resources.files("dpivsoft") / "GPU_Kernels/find_peak.cl"
def _compile_find_peak(box_size_x, box_size_y, window_x, window_y):
"""Compile the find_peak kernel for a given window/search size,
clamping the local size to the kernel's real CL_KERNEL_WORK_GROUP_SIZE
(same rusticl per-kernel limit as SubMean). Returns (kernel, local_size)."""
ds_x = max(2, (box_size_x - window_x) // 2)
ds_y = max(2, (box_size_y - window_y) // 2)
window_cols = box_size_x - 2 * ds_x
window_rows = box_size_y - 2 * ds_y
window_pixels = window_cols * window_rows
cap = max_wg
while True:
local_size = _next_pow2(min(window_pixels, cap))
items_per_thread = (window_pixels + local_size - 1) // local_size
with open(find_peak_path, "r") as f:
prog = thr.compile(f.read(), render_kwds=dict(
box_size_x=box_size_x, box_size_y=box_size_y,
box_pixels=box_size_x * box_size_y,
ds_x=ds_x, ds_y=ds_y,
window_cols=window_cols, window_rows=window_rows,
window_pixels=window_pixels,
local_size=local_size,
items_per_thread=items_per_thread,
lm=4))
kernel = prog.find_peak
limit = kernel.max_work_group_size
if local_size <= limit:
return kernel, local_size
cap = limit # recompile clamped to the kernel's real limit
GPU.FindPeak1, GPU.FindPeak1_local = _compile_find_peak(
Parameters.box_size_1_x, Parameters.box_size_1_y,
Parameters.window_1_x, Parameters.window_1_y)
GPU.FindPeak2, GPU.FindPeak2_local = _compile_find_peak(
Parameters.box_size_2_x, Parameters.box_size_2_y,
Parameters.window_2_x, Parameters.window_2_y)
# Initialize GPU computations for the cross-correlation
GPU.axes = (1, 2)
GPU.fft = FFT(subImg1, axes=GPU.axes).compile(thr)
GPU.fftshift = FFTShift(subImg1, axes=GPU.axes).compile(thr)
GPU.fft2 = FFT(subImg2, axes=GPU.axes).compile(thr)
GPU.fftshift2 = FFTShift(subImg2, axes=GPU.axes).compile(thr)
return 0
[docs]
def processing(img1_name, img2_name, thr):
"""
Perform the parallelized two-pass PIV algorithm with window
deformation on the GPU (OpenCL).
Developed by Jorge Aguilar-Cabello
The image pair to process must already be in GPU memory (GPU.img1,
GPU.img2). While the GPU queue is running, the *next* image pair is
loaded from disk and uploaded, so file I/O overlaps with the
computation.
Parameters
----------
img1_name, img2_name : str
Paths to the image pair of the next iteration, loaded
asynchronously during runtime.
thr : reikna.cluda Thread
Device thread returned by select_Platform.
Notes
-----
Processing options are taken from the Parameters class in
"Classes.py"; they can be set manually or loaded from a file with
the classmethod "readParameters".
Results are left in GPU memory, as attributes of the GPU class in
"Classes.py" (retrieve them with .get()):
GPU.x1, GPU.y1, GPU.u1, GPU.v1 : GPU 2d arrays
Grid and velocity field of the first pass.
GPU.x2, GPU.y2, GPU.u2_f, GPU.v2_f : GPU 2d arrays
Grid and median-filtered velocity field of the second pass
(the final result).
"""
if Parameters.weighting:
Parameters.weighting = False
print("=======================================================================================")
print("Warning: There is a bug in weighting function on the GPU. Parameter.weighting set to 0")
print("=======================================================================================")
N_boxes_1 = Parameters.no_boxes_1_x * Parameters.no_boxes_1_y
N_boxes_2 = Parameters.no_boxes_2_x * Parameters.no_boxes_2_y
N_pixels_1 = N_boxes_1 * Parameters.box_size_1_x * Parameters.box_size_1_y
N_pixels_2 = N_boxes_2 * Parameters.box_size_2_x * Parameters.box_size_2_y
# Zero the velocity-index buffers: the temp manager shares their memory
# with second-pass arrays, so they may hold stale data from a previous run
GPU.ini_index(GPU.u_index_1, GPU.v_index_1, local_size=None,
global_size=N_pixels_1)
# Mask images if required
if Parameters.mask:
GPU.masking(GPU.img1, GPU.img1, GPU.mask, local_size=None,
global_size=grid.pixels)
GPU.masking(GPU.img2, GPU.img2, GPU.mask, local_size=None,
global_size=grid.pixels)
# Apply Gaussian blur (first pass only) if required
if Parameters.gaussian_size:
GPU.gaussian_filter(GPU.img1_g, GPU.img1, GPU.kernel, GPU.data1,
local_size=None, global_size=grid.pixels)
GPU.gaussian_filter(GPU.img2_g, GPU.img2, GPU.kernel, GPU.data1,
local_size=None, global_size=grid.pixels)
thr.synchronize()
# Obtain SubImage
GPU.Slice(GPU.subImg1_1, GPU.img1_g, GPU.box_origin_x_1, GPU.box_origin_y_1,
GPU.data1, local_size=None, global_size=N_pixels_1)
GPU.Slice(GPU.subImg1_2, GPU.img2_g, GPU.box_origin_x_1, GPU.box_origin_y_1,
GPU.data1, local_size=None, global_size=N_pixels_1)
else:
# Obtain SubImage
GPU.Slice(GPU.subImg1_1, GPU.img1, GPU.box_origin_x_1, GPU.box_origin_y_1,
GPU.data1, local_size=None, global_size=N_pixels_1)
GPU.Slice(GPU.subImg1_2, GPU.img2, GPU.box_origin_x_1, GPU.box_origin_y_1,
GPU.data1, local_size=None, global_size=N_pixels_1)
for i in range(0, Parameters.no_iter_1):
# First iteration using direct cross correlation if needed
if not i and Parameters.direct_calc:
# Direct correlation
GPU.directCorrelation(GPU.img1, GPU.img2, GPU.directCorr,
GPU.box_origin_x_d, GPU.box_origin_y_d, GPU.data1,
local_size=None, global_size=int(GPU.size_direct))
# Find peak
GPU.find_peak_direct(GPU.v1, GPU.u1, GPU.directCorr, GPU.data1,
local_size=None, global_size=N_boxes_1)
if i or Parameters.direct_calc:
# Median Filter
GPU.Median_Filter(GPU.u1_f, GPU.v1_f, GPU.u1, GPU.v1,
GPU.median_limit, GPU.data1, local_size=None,
global_size=N_boxes_1)
# Velocity=0 inside mask to prevent bleeding from median filter
if Parameters.mask:
GPU.masking(GPU.u1_f, GPU.u1_f, GPU.mask_1,
local_size=None, global_size=N_boxes_1)
GPU.masking(GPU.v1_f, GPU.v1_f, GPU.mask_1,
local_size=None, global_size=N_boxes_1)
# Jacobian matrix
GPU.jacobian(GPU.temp_dx_1, GPU.temp_dy_1, GPU.u1_f, GPU.x1, GPU.y1,
GPU.data1, local_size=None, global_size=N_boxes_1)
GPU.box_blur(GPU.du_dx_1, GPU.temp_dx_1, GPU.data1,
local_size=None, global_size=N_boxes_1)
GPU.box_blur(GPU.du_dy_1, GPU.temp_dy_1, GPU.data1,
local_size=None, global_size=N_boxes_1)
GPU.jacobian(GPU.temp_dx_1, GPU.temp_dy_1, GPU.v1_f, GPU.x1, GPU.y1,
GPU.data1, local_size=None, global_size=N_boxes_1)
GPU.box_blur(GPU.dv_dx_1, GPU.temp_dx_1, GPU.data1,
local_size=None, global_size=N_boxes_1)
GPU.box_blur(GPU.dv_dy_1, GPU.temp_dy_1, GPU.data1,
local_size=None, global_size=N_boxes_1)
# Deformed image
GPU.deform_image(GPU.subImg1_1, GPU.subImg1_2, GPU.img1, GPU.img2,
GPU.box_origin_x_1, GPU.box_origin_y_1, GPU.u1_f,
GPU.v1_f, GPU.du_dx_1, GPU.du_dy_1, GPU.dv_dx_1,
GPU.dv_dy_1, GPU.u_index_1, GPU.v_index_1, GPU.data1,
local_size=None, global_size=N_pixels_1)
# Normalize
GPU.SubMean1(GPU.subMean1_1, GPU.subImg1_1,
local_size=GPU.SubMean1_local,
global_size=N_boxes_1 * GPU.SubMean1_local)
GPU.SubMean1(GPU.subMean1_2, GPU.subImg1_2,
local_size=GPU.SubMean1_local,
global_size=N_boxes_1 * GPU.SubMean1_local)
GPU.Normalize(GPU.subImg1_1, GPU.subMean1_1, GPU.data1,
local_size=None, global_size=N_pixels_1)
GPU.Normalize(GPU.subImg1_2, GPU.subMean1_2, GPU.data1,
local_size=None, global_size=N_pixels_1)
# Weighting if required
if Parameters.weighting:
GPU.Weighting(GPU.subImg1_1, GPU.data1, local_size=None,
global_size=N_pixels_1)
GPU.Weighting(GPU.subImg1_2, GPU.data1, local_size=None,
global_size=N_pixels_1)
# FFT2D
GPU.fft(GPU.subImg1_1, GPU.subImg1_1)
GPU.fft(GPU.subImg1_2, GPU.subImg1_2)
# Conjugate
GPU.conjugate(GPU.subImg1_1, local_size=None, global_size=N_pixels_1)
# Multiplication
GPU.multiply_them(GPU.subImg1_1, GPU.subImg1_1, GPU.subImg1_2,
local_size=None, global_size=N_pixels_1)
# Inverse transform
GPU.fft(GPU.subImg1_1, GPU.subImg1_1, inverse=True)
# FFTShift
GPU.fftshift(GPU.subImg1_1, GPU.subImg1_1)
# Find peak
GPU.FindPeak1(GPU.v1, GPU.u1, GPU.subImg1_1, GPU.u_index_1,
GPU.v_index_1, GPU.data1,
local_size=GPU.FindPeak1_local,
global_size=N_boxes_1 * GPU.FindPeak1_local)
# Interpolate velocity results from first mesh
GPU.interpolate(GPU.u2_f, GPU.u1, GPU.x2, GPU.y2, GPU.x1, GPU.y1,
GPU.data1, local_size=None, global_size=N_boxes_2)
GPU.interpolate(GPU.v2_f, GPU.v1, GPU.x2, GPU.y2, GPU.x1, GPU.y1,
GPU.data1, local_size=None, global_size=N_boxes_2)
for i in range(0, Parameters.no_iter_2):
# Jacobian matrix
GPU.jacobian(GPU.temp_dx_2, GPU.temp_dy_2, GPU.u2_f, GPU.x2, GPU.y2,
GPU.data2, local_size=None, global_size=N_boxes_2)
GPU.box_blur(GPU.du_dx_2, GPU.temp_dx_2, GPU.data2, local_size=None,
global_size=N_boxes_2)
GPU.box_blur(GPU.du_dy_2, GPU.temp_dy_2, GPU.data2, local_size=None,
global_size=N_boxes_2)
GPU.jacobian(GPU.temp_dx_2, GPU.temp_dy_2, GPU.v2_f, GPU.x2, GPU.y2,
GPU.data2, local_size=None, global_size=N_boxes_2)
GPU.box_blur(GPU.dv_dx_2, GPU.temp_dx_2, GPU.data2, local_size=None,
global_size=N_boxes_2)
GPU.box_blur(GPU.dv_dy_2, GPU.temp_dy_2, GPU.data2, local_size=None,
global_size=N_boxes_2)
# Deformed image
GPU.deform_image(GPU.subImg2_1, GPU.subImg2_2, GPU.img1, GPU.img2,
GPU.box_origin_x_2, GPU.box_origin_y_2, GPU.u2_f,
GPU.v2_f, GPU.du_dx_2, GPU.du_dy_2, GPU.dv_dx_2,
GPU.dv_dy_2, GPU.u_index_2, GPU.v_index_2, GPU.data2,
local_size=None, global_size=N_pixels_2)
# Normalize
GPU.SubMean2(GPU.subMean2_1, GPU.subImg2_1,
local_size=GPU.SubMean2_local,
global_size=N_boxes_2 * GPU.SubMean2_local)
GPU.SubMean2(GPU.subMean2_2, GPU.subImg2_2,
local_size=GPU.SubMean2_local,
global_size=N_boxes_2 * GPU.SubMean2_local)
GPU.Normalize(GPU.subImg2_1, GPU.subMean2_1, GPU.data2,
local_size=None, global_size=N_pixels_2)
GPU.Normalize(GPU.subImg2_2, GPU.subMean2_2, GPU.data2,
local_size=None, global_size=N_pixels_2)
# Weighting if required
if Parameters.weighting:
GPU.Weighting(GPU.subImg2_1, GPU.data2, local_size=None,
global_size=N_pixels_2)
GPU.Weighting(GPU.subImg2_2, GPU.data2, local_size=None,
global_size=N_pixels_2)
# FFT2D
GPU.fft2(GPU.subImg2_1, GPU.subImg2_1)
GPU.fft2(GPU.subImg2_2, GPU.subImg2_2)
# Conjugate
GPU.conjugate(GPU.subImg2_1, local_size=None, global_size=N_pixels_2)
# Multiplication
GPU.multiply_them(GPU.subImg2_1, GPU.subImg2_1, GPU.subImg2_2,
local_size=None, global_size=N_pixels_2)
# Inverse transform
GPU.fft2(GPU.subImg2_1, GPU.subImg2_1, inverse=True)
# FFTShift
GPU.fftshift2(GPU.subImg2_1, GPU.subImg2_1)
# Find peak
GPU.FindPeak2(GPU.v2, GPU.u2, GPU.subImg2_1, GPU.u_index_2,
GPU.v_index_2, GPU.data2,
local_size=GPU.FindPeak2_local,
global_size=N_boxes_2 * GPU.FindPeak2_local)
# Median Filter
GPU.Median_Filter(GPU.u2_f, GPU.v2_f, GPU.u2, GPU.v2,
GPU.median_limit, GPU.data2, local_size=None,
global_size=N_boxes_2)
if Parameters.mask:
# Check if inside mask to prevent bleeding from median filter
GPU.masking(GPU.u2_f, GPU.u2_f, GPU.mask_2,
local_size=None, global_size=N_boxes_2)
GPU.masking(GPU.v2_f, GPU.v2_f, GPU.mask_2,
local_size=None, global_size=N_boxes_2)
# Load the next iteration's images while the GPU queue is still running
Img1, Img2 = DPIV.load_images(img1_name, img2_name)
thr.synchronize()
# Send next iteration images to the GPU (overwrite pre-allocated buffers)
GPU.img1.set(Img1)
GPU.img2.set(Img2)
return 0