# Syntethic image Generation
import numpy as np
import cv2
import random
import os
import sys
# import vtk
import shapely.geometry as shapely
from scipy import interpolate
import dpivsoft.meshTools as mt
# Image Parameteres import. Needs to be changed on the class
from dpivsoft.Classes import Synt_Img
# Generate images from a given analytical flow field
[docs]
def Analytic_Syntetic(dirSave, Name, saveData=True, randTransform=False):
"""
Generates a pair of imnages where the trazers particles moves
accordly to an analytical flow velocity field. The parameters of
the generated images are defined on the class Synt_Img (see class
Synt_Img for more info).
"""
# Generate trazers
trazers_shine, trazers_D = Gen_trazers()
# Velocity profile function
x, y, u, v = Velocity_Profile()
if randTransform:
x, y, u, v = random_transformation(x, y, u, v)
img1, img2 = Img_Generation(u, v, x, y, trazers_shine, trazers_D)
# Save Images
SaveImg(img1, img2, dirSave, Name)
# Save Data
if saveData:
np.savez(dirSave + '/Py_profile_' + Name, x=x, y=y, u=u, v=v)
return u, v
[docs]
def Custom_Syntetic(xx, yy, uu, vv, scale, dirSave, Name, dt, limits=None, randTransform=False):
"""
Generates a pair of syntetic images where the trazers particles move
accordly a custom velocity field loaded from a numpy array.
The numpy array must be saved in colums as: | x | y | u | v |
"""
if limits is None:
x_min = np.min(xx)
x_max = np.max(xx)
y_min = np.min(yy)
y_max = np.max(yy)
else:
# Limits of the mesh
x_min = limits[0]
x_max = limits[1]
y_min = limits[2]
y_max = limits[3]
# Generate trazers
trazers_shine, trazers_D = Gen_trazers()
# Image pixels mesh in problem units to interpolate
xv = np.linspace(x_min, x_max, Synt_Img.width) / scale
yv = np.linspace(y_min, y_max, Synt_Img.height) / scale
xv, yv = np.meshgrid(xv, yv)
# Interpolate velocity into each pixel
u = interpolate.griddata((xx / scale, yy / scale), uu * dt / scale,
(xv, yv), method='linear')
v = interpolate.griddata((xx / scale, yy / scale), vv * dt / scale,
(xv, yv), method='linear')
# Get x,y mesh in pixel for image generation
xv, yv = np.meshgrid(np.arange(0, Synt_Img.width),
np.arange(0, Synt_Img.height))
if randTransform:
xv, yv, u, v = random_transformation(xv, yv, u, v)
# Generate Image
img1, img2 = Img_Generation(u, v, xv, yv, trazers_shine, trazers_D)
# Save Images
SaveImg(img1, img2, dirSave, Name)
return u, v
# Generation of the selected analytical velocity field to create the
# PIV image pairs
[docs]
def Velocity_Profile():
"""
Build the analytical velocity field used to advect the tracers.
The profile is selected by ``Synt_Img.vel_profile`` and evaluated on
the pixel grid: 'Constant', 'Couette', 'Poiseuille', 'Vortex'
(Lamb-Oseen-like) or 'Frequency' (sinusoidal shear). Exits if the
profile name is unknown.
Returns
-------
xv, yv : 2d np.ndarray
Pixel-coordinate meshgrid.
u, v : 2d np.ndarray
Velocity components on that grid.
"""
N_pixel = Synt_Img.width * Synt_Img.height
vel = Synt_Img.vel
# Mesh generation
x = np.linspace(0, Synt_Img.width - 1, Synt_Img.width)
y = np.linspace(0, Synt_Img.height - 1, Synt_Img.height)
xv, yv = np.meshgrid(x, y)
x = xv.reshape((N_pixel), order='C') # X axis write on 1-D mode
y = yv.reshape((N_pixel), order='C') # Y axis write on 1-D mode
if Synt_Img.vel_profile == 'Constant':
v = np.zeros([Synt_Img.height, Synt_Img.width])
u = vel * np.ones([Synt_Img.height, Synt_Img.width])
elif Synt_Img.vel_profile == 'Couette':
v = np.zeros([Synt_Img.height, Synt_Img.width])
u = vel * yv / np.max(np.max(yv))
elif Synt_Img.vel_profile == 'Poiseuille':
h = np.max(np.max(yv))
v = np.zeros([Synt_Img.height, Synt_Img.width])
u = vel * yv * (h - yv) / h / h
elif Synt_Img.vel_profile == 'Vortex':
m_xv = np.mean(np.mean(xv))
m_yv = np.mean(np.mean(yv))
xv = xv - m_xv
yv = yv - m_yv
R0 = 200
r = np.sqrt(xv**2 + yv**2)
u = -vel * yv / R0 / (1 + (r / R0)**2)
v = vel * xv / R0 / (1 + (r / R0)**2)
xv = xv + m_xv
yv = yv + m_yv
elif Synt_Img.vel_profile == 'Frequency':
u = vel * np.sin(2 * np.pi * yv / Synt_Img.amp)
v = np.zeros([Synt_Img.height, Synt_Img.width])
else:
sys.exit("Velocity profile not found")
return xv, yv, u, v
# Generation of the trazers
[docs]
def Gen_trazers():
"""
Randomly sample the brightness and diameter of the tracer particles.
Both are drawn from triangular distributions centred on the class
means (``Shine_m``, ``D_m``) with half-widths ``d_Shine`` / ``d_D``,
one value per particle.
Returns
-------
trazers_shine : 1d np.ndarray
Peak brightness of each particle.
trazers_D : 1d np.ndarray
Diameter of each particle (pixels).
"""
N_trazers = round(Synt_Img.width * Synt_Img.height * Synt_Img.trazers_density)
# particles light
trazers_shine = np.random.triangular(Synt_Img.Shine_m - Synt_Img.d_Shine,
Synt_Img.Shine_m, Synt_Img.Shine_m + Synt_Img.d_Shine, N_trazers)
# particles dimaeter
trazers_D = np.random.triangular(Synt_Img.D_m - Synt_Img.d_D, Synt_Img.D_m,
Synt_Img.D_m + Synt_Img.d_D, N_trazers)
return trazers_shine, trazers_D
# Create the image pair with particles from the velocity field in 1-D
[docs]
def Img_Generation(u, v, xv, yv, trazers_shine, trazers_D):
"""
Render the synthetic image pair from a velocity field and tracers.
Particles are seeded at random subpixel positions for the first
image, then advected to their second-image positions by integrating
the velocity field with RK4 over ``Synt_Img.n_steps`` substeps.
Particles are drawn as Gaussian spots; those leaving the frame
re-enter periodically from the opposite side. Triangular-distributed
background noise is added and intensities are clipped to [0, 255].
Parameters
----------
u, v : 2d np.ndarray
Per-pixel displacement field (pixels between the two exposures).
xv, yv : 2d np.ndarray
Pixel-coordinate meshgrid matching ``u``, ``v``.
trazers_shine, trazers_D : 1d np.ndarray
Per-particle brightness and diameter (from Gen_trazers).
Returns
-------
img1, img2 : 2d np.ndarray
The synthetic image pair.
"""
# Variables
width = Synt_Img.width
height = Synt_Img.height
trazers_density = Synt_Img.trazers_density
# randomly distributed trazers inside an image
trazer = np.array(random.sample(range(0, width * height),
round(width * height * trazers_density)))
# Add subpixel offsets so img1 and img2 have the same statistical
# distribution of Gaussian peak positions (avoids systematic brightness
# differences caused by img2 float positions vs img1 integer positions)
x_1 = (trazer) % width + np.random.uniform(0, 1, len(trazer))
y_1 = (trazer) // width + np.random.uniform(0, 1, len(trazer))
interp_u = interpolate.RegularGridInterpolator(
(yv[:, 0], xv[0, :]), u,
bounds_error=False, fill_value=None)
interp_v = interpolate.RegularGridInterpolator(
(yv[:, 0], xv[0, :]), v,
bounds_error=False, fill_value=None)
def velocity_at(x, y):
u = interp_u((y, x))
v = interp_v((y, x))
return u, v
def rungeKuta4(x, y, dt):
u1, v1 = velocity_at(x, y)
u2, v2 = velocity_at(x + 0.5 * dt * u1, y + 0.5 * dt * v1)
u3, v3 = velocity_at(x + 0.5 * dt * u2, y + 0.5 * dt * v2)
u4, v4 = velocity_at(x + dt * u3, y + dt * v3)
x_new = x + dt * (u1 + 2 * u2 + 2 * u3 + u4) / 6
y_new = y + dt * (v1 + 2 * v2 + 2 * v3 + v4) / 6
return x_new, y_new
dt = 1.0 / Synt_Img.n_steps
x_2 = x_1.astype(np.float64)
y_2 = y_1.astype(np.float64)
for _ in range(Synt_Img.n_steps):
x_2, y_2 = rungeKuta4(x_2, y_2, dt)
# Image inicializate
img1 = np.zeros([height, width])
img2 = np.zeros([height, width])
# Particles added to images
for i in range(0, len(trazers_D)):
x1_p = x_1[i]
y1_p = y_1[i]
xp = int(x1_p)
yp = int(y1_p)
td = int(trazers_D[i])
# Gaussian distribution bright for each particle on the first image
y0, y1, x0, x1 = max(0, yp-td), min(height, yp+td), max(0, xp-td), min(width, xp+td)
img1[y0:y1, x0:x1] = (img1[y0:y1, x0:x1] +
trazers_shine[i] * np.exp(-((xv[y0:y1, x0:x1]
- x1_p)**2 + (yv[y0:y1, x0:x1] - y1_p)**2)
* 4 / (trazers_D[i]**2)))
# Interpolated position of particle for second image
x2_p = x_2[i]
y2_p = y_2[i]
# Rounded center of particle on second image
xp = x_2[i].astype(int)
yp = y_2[i].astype(int)
# Detect if a particle left the image after displacement to
# enter throught the other side
if xp < width and yp < height and yp > 0 and xp > 0:
# Gaussian distribution bright for each particle on the
# second image
y0, y1, x0, x1 = max(0, yp-td), min(height, yp+td), max(0, xp-td), min(width, xp+td)
img2[y0:y1, x0:x1] = (img2[y0:y1, x0:x1] +
trazers_shine[i] * np.exp(-((xv[y0:y1, x0:x1]
- x2_p)**2 + (yv[y0:y1, x0:x1] - y2_p)**2)
* 4 / trazers_D[i]**2))
else:
# Add randomness
if x2_p >= width:
xp = xp - width
x2_p = x2_p - width
if y2_p >= height:
yp = yp - height
y2_p = y2_p - height
if x2_p < 0:
xp = xp + width
x2_p = x2_p + width
if y2_p < 0:
yp = yp + height
y2_p = y2_p + height
# Gaussian distribution bright for each particle on the
# second image
y0, y1, x0, x1 = max(0, yp-td), min(height, yp+td), max(0, xp-td), min(width, xp+td)
img2[y0:y1, x0:x1] = (img2[y0:y1, x0:x1] +
trazers_shine[i] * np.exp(-((xv[y0:y1, x0:x1]
- x2_p)**2 + (yv[y0:y1, x0:x1] - y2_p)**2)
* 4 / (trazers_D[i]**2)))
# Final image adding random noise distribution
img_noise = np.random.triangular(Synt_Img.noise_m - Synt_Img.d_noise,
Synt_Img.noise_m, Synt_Img.noise_m + Synt_Img.d_noise,
(Synt_Img.height, Synt_Img.width))
img1 = img1 + img_noise
img_noise = np.random.triangular(Synt_Img.noise_m - Synt_Img.d_noise,
Synt_Img.noise_m, Synt_Img.noise_m + Synt_Img.d_noise,
(Synt_Img.height, Synt_Img.width))
img2 = img2 + img_noise
img1[img1 > 255] = 255
img1[img1 < 0] = 0
img2[img2 > 255] = 255
img2[img2 < 0] = 0
return img1, img2
# This funciton transform velocity field used to generate the testing images
# to avarage value inside the deformation windows introduced from the PIV
[docs]
def Pix2PIV(Xv, Yv, Uv, Vv, no_boxes_x, no_boxes_y, box_size_1_x,
box_size_1_y, box_size_2_x, box_size_2_y):
"""
Downsample a per-pixel velocity field onto the coarse PIV grid.
Averages the pixel-level field over each interrogation window,
reproducing the window-averaging inherent to classical PIV so a
ground-truth field can be compared with a PIV result (or used as the
coarse prior in the Stage 2 PIV-conditioned DL model). The grid is
laid out with the same margin/spacing as ``grid.generate_mesh``.
Parameters
----------
Xv, Yv : 2d np.ndarray
Pixel-coordinate meshgrid.
Uv, Vv : 2d np.ndarray
Per-pixel velocity components.
no_boxes_x, no_boxes_y : int
Number of interrogation windows in x and y.
box_size_1_x, box_size_1_y, box_size_2_x, box_size_2_y : int
First- and second-pass window sizes (first pass sets the margin,
second pass sets the averaging window).
Returns
-------
Xi, Yi : 1d np.ndarray
PIV grid coordinates.
Ui, Vi : 2d np.ndarray
Window-averaged velocity components on the PIV grid.
"""
Ui = np.zeros([no_boxes_y, no_boxes_x])
Vi = np.zeros([no_boxes_y, no_boxes_x])
Xi = np.zeros([no_boxes_y, no_boxes_x])
Yi = np.zeros([no_boxes_y, no_boxes_x])
box_origin_x = np.zeros([no_boxes_y, no_boxes_x])
box_origin_y = np.zeros([no_boxes_y, no_boxes_x])
Height, Width = Xv.shape
x_margin = 3 / 2 * np.amax([box_size_1_x, box_size_2_x])
y_margin = 3 / 2 * np.amax([box_size_1_y, box_size_2_y])
Xi = 2 + np.round(np.arange(0, no_boxes_x) *
(Width - x_margin - 6) / (no_boxes_x - 1) + x_margin / 2)
Yi = 2 + np.round(np.arange(0, no_boxes_y) *
(Height - y_margin - 6) / (no_boxes_y - 1) + y_margin / 2)
box_origin_x = Xi - box_size_2_x / 2
box_origin_y = Yi - box_size_2_y / 2
for i in range(0, no_boxes_x):
for j in range(0, no_boxes_y):
# Define sub images to work
Sub_u = (Uv[int(box_origin_y[j]):int(box_origin_y[j]) + box_size_2_y,
int(box_origin_x[i]):int(box_origin_x[i]) + box_size_2_x])
Sub_v = (Vv[int(box_origin_y[j]):int(box_origin_y[j]) + box_size_2_y,
int(box_origin_x[i]):int(box_origin_x[i]) + box_size_2_x])
Ui[j, i] = np.mean(Sub_u)
Vi[j, i] = np.mean(Sub_v)
return Xi, Yi, Ui, Vi
[docs]
def SaveImg(Img1, Img2, dirSave, Name):
"""
Save the syntetic Images
"""
cv2.imwrite(dirSave + '/' + Name + '1' + Synt_Img.ext, Img1)
cv2.imwrite(dirSave + '/' + Name + '2' + Synt_Img.ext, Img2)