To reproduce the error
Module 1
"""
module ImageEnhancement
This module contains various compensation methods that can be applied to the
creation of PAT images.
"""
module ImageEnhancement
"""
find_relative_c0_map(args)
This function will take a speed of sound map, which is generally formed from
a UST system, and find the relative speed of sound between a transducer and
a given grid point. It does so by fitting a line y=m*x+b between the
transducer and the grid point of interest. It then uses these fitted values
and integrates along the path length to find the relative/mean speed of sound.
Input:
Output:
"""
function find_relative_c0_map(c0_map::Array{Float64,2},
grid_index::Array{CartesianIndex{2},1},
ImagingGrid)
elements = grid_index
relative_c0_map = zeros(ImagingGrid.Nx, ImagingGrid.Ny)
for element in elements
x1 = element[1]
y1 = element[2]
for x2 in 1:ImagingGrid.Nx
for y2 in 1:ImagingGrid.Ny
if x2 == x1 && y2 == y1
#=
If the currently selected pixel is at the same location
as the receiver then skip it.
=#
continue
else
point1 = CartesianIndex(x1,y1)
point2 = CartesianIndex(x2,y2)
mean_val = c0_mean_path(c0_map, point1, point2)
relative_c0_map[x2,y2] = mean_val
end
end
end
end
return relative_c0_map
end
"""
c0_mean_path(args)
This fucntion integrates along the path length between the receiver and the
pixel of interest.
Inputs:
c0_map:
point1:
point2:
Outputs:
mean_val:
Type: Float64
"""
function c0_mean_path(c0_map, point1::CartesianIndex{2},
point2::CartesianIndex{2})
include("PATRecon.jl")
x1 = point1[1]; y1 = point1[2]
x2 = point2[1]; y2 = point2[2]
m = PATRecon.Utilities.find_slope(x1,y1,x2,y2)
b = PATRecon.Utilities.find_y_inter(x1, y1, m)
#=
Find the number of points that we have to step to find the mean path.
The addition of +1 is because we are looking for the number of points
contained in this distance.
=#
num_x = x1-x2
num_y = y1-y2
num_points = â(num_x^2 + num_y^2)+1
mean_value = 0
#=
Start selection tree to fit the line between the points (x1,y1)
and (x2, y2).
=#
if m == Inf || m == -Inf
#=
When the slope tends towards +â
Which occurs when x1 == x2 and y2-y1 > 0
or
When the slope tends towards -â
Which occurs when x1 == x2 and y2-y1 < 0
=#
for step in range(y1,y2,length=num_points)
mean_val += c0_map[x2, step]
end
elseif x2-x1 < 0 || x2-x1 > 0
#=
x2-x1 < 0
tells us that the pixel is to the left of the receiver.
or
x2-x1 > 0 tells us that the pixel is to the right of the receiver.
=#
for step in range(x1,x2,length=num_points)
y_cal = PATRecon.Utilities.find_y(step,m,b)
#=
Round the calculated y value down unless it is zero in which
case set the value to one, since indexing starts at one.
=#
y_cal = Int64(floor(y_cal))
if y_cal == 0
y_cal = 1
end
mean_val += co_map[step, y_cal]
end
else
print("ERROR: Unknown condiation meet when fitting a line to points
(x1:$x1, y1:$y1) and (x2:$x2, y2:$y2)")
end
# Divide by how many steps we took to get the average value.
mean_val /= num_points
return mean_val::Float64
end
end # module CompMethods
Module 2
"""
module Utilities
This module holds all utilies that can be shared through the code base.
"""
module Utilities
"""
cart2grid(args)
This function converts a set of cartesian coordinates to a matrix of grid
points.
Inputs:
positions:
Type: Array{Float64,2}
The array of position for each element in 'mm' where...
position[element#, posX==1 or posY==2]
i.e. position[1,1] will give you the position for the first
element in the x direction.
dx:
The spacing in mm/grid point for the x-direction.
dy:
The spacing in mm/grid point for the y-direction.
Nx:
The number of total grid points in the x-direction.
Ny:
The number of total grid points in the y-direction.
Outputs:
grid_data:
Type: Array{Int64,2}
A matrix that contains the transducer positions, where the location
of the transducer is its number i.e. transducer one has the value 1
within the matrix grid_data.
grid_index:
Type: Array{CartesianIndex{2},1}
An array that has the indexes of the transducers, where transducer
one is index one.
example:
grid_index[1] = CartesianIndex(0, 0)
This is the position of the first transducer.
"""
function cart2grid(positions::Array{Float64,2}, ImagingGrid)
# Get positional data.
x_data = positions[:,1] # in [mm]
y_data = positions[:,2] # in [mm]
# Convert to grid points to the nearest integer.
x_data = Int64.(round.(x_data./ImagingGrid.dx)) # in [grid points]
y_data = Int64.(round.(y_data./ImagingGrid.dy)) # in [grid points]
# Shift the grid point indexes so they aline with the imaging grid.
x_data .+= floor(ImagingGrid.Nx/2) + 1
y_data .+= floor(ImagingGrid.Ny/2) + 1
# Verfiy that all points lie within the grid.
if maximum(x_data) > ImagingGrid.Nx || maximum(y_data) > ImagingGrid.Ny ||
minimum(x_data) < 1 || minimum(y_data) < 1
print("ERROR: The transducer positions must lie within the
imaging grid.")
end
# Assign transducer numbers.
transducer_nums = 1:size(x_data,1)
# Allocate matrix for transducer element position.
grid_data = zeros(Int64,ImagingGrid.Nx,ImagingGrid.Ny)
# Allocate array for transducer element indexes.
grid_index = zeros(CartesianIndex{2},length(transducer_nums))
# Assign transducer number to grid.
for idx in range(1,length=size(x_data,1))
# Mark transducer position in grid.
grid_data[x_data[idx],y_data[idx]] = transducer_nums[idx]
# Save position of transducer.
grid_index[idx] = findall(grid_data.==idx)[1]
end
# Rotate the grid to the left 90° so that outputted matrixes upper
# right side is +x and +y while the lower left is -x and -y.
# Just like a normal cartesian coordinate system.
grid_data = rotl90(grid_data)
return grid_data::Array{Int64,2}, grid_index::Array{CartesianIndex{2},1}
end
"""
find_slope(x1,y1,x2,y2)
Finds the slope for two sets of points (x1,y1) and (x2,y2).
Inputs:
x1,y1:
The coordinates of the first grid point in the matrix. Generally
this is where the transducer is located. This should be given as
the indexes of that matrix.
x2,y2:
The coordinates of the second grid point in the matrix. Generally
this is the pixel whose relative speed of sound we would like to
know.
Outputs:
m:
The calculated slope.
"""
function find_slope(x1,y1,x2,y2)
return m = (y2-y1)/(x2-x1)
end
"""
find_y_inter(y,x,m)
This function will find the y intercept for a pair of points and the
associated slope that was calculated with "find_slope(...)".
Inputs:
x:
The x value from a pair of points. Given (x1,y1) this would be x1.
y:
The y value from a pair of points. Given (x1,y1) this would be y1.
m:
The slope that was calculated from "find_slope(...)" with a set of
points (x1,y1) and (x2,y2).
Outputs:
b:
The y-intercept.
"""
function find_y_inter(x,y,m)
return b = y-m*x
end
"""
find_y(m,x,b)
This function calculates the linear equation y=m*x+b. This function can take
an array of inputs for x in order to avoid looping through x-values.
Inputs:
m:
The slope that was calculated from "find_slope(...)" with a set of
points (x1,y1) and (x2,y2).
x:
The input x that we would like to find the associated y for.
b:
The y-intercept that was found from "find_y_inter(...)".
Outputs:
y:
The given value y for the input x according to the equation
y=m*x+b.
"""
function find_y(m,x,b)
return y = m.*x.+b
end
"""
A structure that holds information regarding the discrete imaging grid.
Values:
dx:
Type: Int64
Unit: mm/grid point
The spacing in mm/grid point for the x-direction.
dy:
Type: Int64
Unit: mm/grid point
The spacing in mm/grid point for the y-direction.
Nx:
Type: Int64
Unit: grid point(s)
The number of total grid points in the x-direction.
Ny:
Type: Int64
Unit: grid point(s)
The number of total grid points in the y-direction.
"""
struct ImagingGrid
dx::Int64
dy::Int64
Nx::Int64
Ny::Int64
end
end # module Utilities
Module 3
"""
module PATRecon
This module imports all of the sub-modules for this library and acts like
the parent module.
"""
module PATRecon
include("Utilities.jl")
include("ImageEnhancement.jl")
using Reexport
@reexport using .Utilities, .ImageEnhancement
end # module PATRecon
Test Script
include("PATRecon.jl")
c0_map = Array{Float64,2}([1400 1400 1400;
1400 1400 1400;
1500 1500 1500])
Nx = size(c0_map,1)
Ny = size(c0_map,2)
dx = dy = 1
mygrid = PATRecon.Utilities.ImagingGrid(dx,dy,Nx,Ny)
grid_index = [CartesianIndex(1,1)]
relative_c0_map = PATRecon.ImageEnhancement.find_relative_c0_map(c0_map, grid_index, mygrid)
print(relative_c0_map)