| SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
Spacekime Analytics (Time Complexity and Inferential Uncertainty)
Image Processing
SOCR Team
01/17/2022
# Get this figure: fig <- get_figure("MattSundquist", 4064)
# Get this figure's data: data <- get_figure("MattSundquist", 4064)$data
# Add data to this figure: p <- add_trace(p, x=c(4, 5), y=c(4, 5), kwargs=list(filename="Klein bottle", fileopt="extend"))
# Get y data of first trace: y1 <- get_figure("MattSundquist", 4064)$data[[1]]$y
library(plotly)1 1D Waves and 2D Waves
Waves of constant wavelength propagate maintaining its shape but shifting its location like in the image below. However, this oversimplified representation does not work well for computing, as adding wave amplitudes does not work as expected. It simply shifts the wave amplitude up and doesn’t move the wave in the direction of motion. Similar problem arises with wave multiplication by a constant. Wave multiplication by a number increases the height of the peaks and valleys without moving them in space.
# install.packages("deSolve")
times <- seq(from = 0, to = 40, by = 0.01)
y1 <- sin(times/2)
y2 <- sin(times/2+1)
y3 <- sin(times/2+2)
plot(times, y1, type="l", lwd=2, col="lightblue")
lines(times, 0*times, type="l", lwd=2, col="black")
lines(times, y2, type="l", lwd=2, col="darkblue")
lines(times, y3, type="l", lwd=2, col="blue")
These problems are solved by an alternative wave representation where the wave has a complex amplitude, and complex multiplication by a number works well.
The wave phase represents the direction of the wave amplitude in complex plane. Differences of phases yield interference phenomena. The interference of two waves of amplitude \(i\) will produce an increased amplitude of \(2i=i + i\). When the second wave is of amplitude \(-i\), then their interfere destroys the aggregate wave, i.e., produces a zero amplitude wave \(i - i = 0\).
Let’s start with a wave of constant wavelength that represents a particle of definite momentum. For simplicity, we’ll consider just a one space dimension representing the horizontal axis, i.e., the spatial extension of the wave. The remaining two dimensions will represent the phase of the wave complex amplitude. The vertical axis represents the imaginary \(+i\) and \(-i\) directions, and depth represents the real \(+1\) and \(-1\) directions. In this representation, the wave magnitude is constant, however, the 2D phase determines the uniform distance advancement in space.
As the wave propagates, the amplitude cycles through all directions in the 2D complex plane perpendicular to the 1D space dimension (x-axis). The result is a wave function represented by a helical (“corkscrew”) shape. In the previous example, the illusion that the wave amplitude periodically vanishes is just an artifact of an incomplete representation of the wave.
The Complex wave representation clarifies this illusion arising by slicing through the 3D “corkscrew” wave shape and showing just the 1D values in the space of complex amplitudes.
This complex wave representation also addresses the above multiplication problem. Multiplying the 3D wave curve by \(i\) will relocate the amplitudes yielding spatial propagation of the wave, which explains the basic rule of time evolution of quantum theory. The effect of multiplying the wave amplitude by \(i\) at each point is effectively a rotation of the amplitude by one quarter turn around the space of values. Thus, the combination of these point-wise effects results in the wave advancing forward in space (propagation).
library(plotly)
t <- seq(-6*pi, 6*pi, length.out=400)
r <- 3
x <- r * cos(t)
y <- r * sin(t)
z <- 2*t
c <- t%%(2*pi)
data1 <- data.frame(x, y, z)
p <- plot_ly(data1, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', showlegend = F,
line = list(width = 12, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')))) %>%
# trace the z-axis
add_trace(data1, x = 0, y =0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 10, color = ~c, colorscale = list(c(0,'red'), c(1,'blue'))),
name="Z", hoverinfo="none") %>%
# add projection on plane x=-3
add_trace(data1, x = -3, y = ~y, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add projection on plane y=-3
add_trace(data1, x = ~x, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add a z-line at x=-3,y=0
add_trace(data1, x = -3, y = 0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a z-line at x=0,y=-3
add_trace(data1, x = 0, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a x-line at y=0, z=-40
add_trace(data1, x = ~x, y = 0, z = -40, type="scatter3d", mode="lines",
line = list(width = 15, dash="solid", color = "gray"),
name="Z", hoverinfo="none")
# add a few annotations/arrows
pRecall that for a single particle, the Hamiltonian operator corresponding to the total energy of the system is defined by:
\[\hat {H}={\underbrace{\frac{\hat {p}^{2}}{2m}}_{\text{particle kinetic energy}}}+ {\underbrace{\frac {1}{2} k \hat {x}^2}_{\text{particle potential energy}}}= {\frac {{\hat {p}}^{2}}{2m}}+{\frac {1}{2}}m\omega ^{2}{\hat {x}}^{2} ,\]
where \(m\) is the particle mass, \(k\) is the force constant, \(w =\sqrt {\frac {k}{m}}\) is the angular frequency of the oscillator, \(\hat {x}\) is the position operator (corresponding to the spatial location \(x\)), and \(\hat {p}\) is the momentum operator (corresponding to the momentum), which is expressed by \(\hat {p}= -i\hbar \frac{\partial}{\partial x}\). Then, the time-independent Schr?dinger equation can be expressed as:
\[\hat {H}\left|\psi \right\rangle =E\left|\psi \right\rangle ,\]
where the eigenvalue \(E\) denotes a real number specifying a time-independent energy level and the wavefunction \(\psi\) represents the solution corresponding to that level’s energy eigenstate.
Another way of visualizing the propagation of a wavefunction over a 1D space in terms of its complex amplitude (\(\psi =Ae^{i(kx-wt)}\)), real and imaginary parts, where \(t\) is time, \(w = kc\), wave’s angular frequency which equals \(w=2\pi/T\), where \(T\) is the period of the wave, the position \(x\), and momentum \(p\) are conjugate variables. The wavefunction magnitude represents the probability of finding the particle at a given point \(x\in R^1\) is spread out like a waveform. There is no unique definite position of the particle. As the amplitude oscillates (increases and decreases), the result is a wave with an alternating amplitude.
\(k\), the wave’s radian spatial frequency (angular wave number), measures how rapidly the disturbance changes over a given distance at a particular point in time, \(k=2\pi/\lambda\), where \(\lambda\) is the wavelength of the wave and measured in radians per unit distance. Note that \((x,y,z)\), are not part of the equation because the wave’s magnitude and phase are the same at every point.
library(plotly)
t <- seq(-6*pi, 6*pi, length.out=400) # repeat the periodic wave 6 times
r <- 3 # Wave Ammplitude (radius)
x <- r * cos(t) # X coordinate
y <- r * sin(t) # Y coordinate
z <- 2*t # stretched (*2) time-axis (z)
c <- t%%(2*pi) # wave color
data1 <- data.frame(x, y, z)
scene = list(camera = list(eye = list(x = -1.25, y = 1.25, z = 1.25)),
aspectmode = "manual", aspectratio = list(x=1, y=1, z=5))
# define wave opacity as a parabola
const <- 1/(z[400]^2)
wave_opacity <- 0.3 -(0.7*const)*(z-z[1])*(z-z[400]); plot(wave_opacity, type = "l")
# bgcolor = toRGB("white", alpha = 0),
#cols = sample(grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)], 400)
#cols_func <- colorRampPalette(c("blue", "red")); cols <- cols_func(400)
#plot_ly(data1, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', showlegend = F,
# line = list(width = 40, color = toRGB(cols, alpha = wave_opacity)))
p <- plot_ly(data1, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', showlegend = F,
line = list(width = 40, color = ~c, colorscale = list(c(0,'red'), c(1,'blue')),
alpha = wave_opacity)) %>%
# trace the z-axis
add_trace(data1, x = 0, y =0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 10, color = ~c, colorscale = list(c(0,'red'), c(1,'blue'))),
name="Z", hoverinfo="none") %>%
# add projection on plane x=-3
add_trace(data1, x = -3, y = ~y, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add projection on plane y=-3
add_trace(data1, x = ~x, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 2, dash="solid", color = ~c,
colorscale = list(c(0,'gray'), c(1,'black'))),
name="Z", hoverinfo="none") %>%
# add a z-line at x=-3,y=0
add_trace(data1, x = -3, y = 0, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a z-line at x=0,y=-3
add_trace(data1, x = 0, y = -3, z = ~z, type="scatter3d", mode="lines",
line = list(width = 5, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
# add a x-line at y=0, z=-40
add_trace(data1, x = ~x, y = 0, z = -40, type="scatter3d", mode="lines",
line = list(width = 15, dash="solid", color = "gray"),
name="Z", hoverinfo="none") %>%
layout(scene = scene)
# add a few annotations/arrows
p- Comment: The probability density of eigenfunctions of time-dependent Schrodinger equation are time independent whereas more general wavefunctions representing combinations of eigenfunctions are typically time-dependent.
This argument is based on the definition of Schrodinger equation eigenfunctions. All eigenfunctions (\(\psi\)) of the Hamiltionian operator (\(\hat{H}\)) satisfy the equation:
\[\hat{H} \psi(x) = E\times \psi(x) ,\]
where eigenvalue \(E\) is the energy of the quantum state associated with the eigenfunction (wavefunction) \(\psi(x)\), which depends only on the spatial location (\(x\)) and not on time \(t\). Whereas the total eigenfunction is a function of both, space and time:
\[\Psi(x, t) = \psi(x) e^{-iEt}. \]
Then, the probability density is: \[P(x, t) =||\Psi||^2=\langle \Psi | \Psi \rangle =\Psi(x, t) \Psi^* (x, t)=\psi(x) e^{-iEt} \psi(x) e^{iEt}=||\psi(x)||^2.\]
Thus, the probability density is a function of a single eigenfunction and only depends on the spatial position.
However, combinations of eigenfunctions are typically time-dependent. A linear combination of a pair of eigenfunctions: \[\phi(x)= a_1\psi_1(x) +a_2\psi_2(x)\] has the following probability density representing the likelihood that the particle, given a specific location and time:
\[P(x, t) = \langle \phi | \phi \rangle = \left( a_1\psi_1(x)e^{-iE_1 t} + a_2\psi_2(x)e^{-iE_2 t} \right) \left( a_1\psi_1(x)e^{-iE_1 t} + a_2\psi_2(x)e^{-iE_2 t} \right)^*=\]
\[ \left( a_1\psi_1(x)e^{-iE_1 t} + a_2\psi_2(x)e^{-iE_2 t} \right) \left( a_1\psi_1(x)e^{iE_1 t} + a_2\psi_2(x)e^{iE_2 t} \right)=\] \[ a_1^2\psi_1^2(x) + a_2^2\psi_2^2(x) + a_1a_2\psi_1(x)\psi_2(x) \left( e^{-it(E_1 - E_2)} + e^{-it(E_2 - E_1)} \right)=\] \[ a_1^2\psi_1^2(x) + a_2^2\psi_2^2(x) + 2a_1a_2\psi_1(x)\psi_2(x) \cos(t(E_1 - E_2)).\]
Thus, unless \(E_1 = E_2\), the the probability density of the linear mixture (\(P(x, t)\)) would be time-dependent.
2 Fourier Transform of images
Let’s demonstrate the use of R’s fft() to calculate the 2D FT of images. Information on data acquisition frequency and block length (in sec) cannot be included into the fft() call.
R generates a single 2D vector of the same dimensions as the data containing a list of complex numbers. The function Mod(fft()) is used to extract the magnitudes of the Fourier coefficients, which are computed by: \[magnitude = \sqrt{(real * real + imaginary * imaginary)},\] where the real and imaginary components are extracted by functions Re(fft()) and Im(fft()), respectively.
The method fft() generates only meaningful frequency up to half the sampling frequency. The FT returns values of the discrete Fourier transform for both positive and negative frequencies. Although, as sampling a signal in discrete time intervals causes aliasing problems, R yields all frequencies up to the sampling frequency. For instance, sampling a 50 Hz sine wave and 950 Hz sine wave with 1000 Hz will generate identical results, as the FT cannot distinguish between the two frequencies. Hence, the sampling frequency must always be at least twice as high as the expected signal frequency. For each actual frequency in the signal, the FT will give 2 peaks (one at the “actual” frequency and one at sampling frequency minus “actual” frequency). This will make the second half of the magnitude vector a mirror image of the first half.
As long as the sampling frequency is at least twice as high as the expected signal frequency, all meaningful information will be contained in the the first half of the magnitude vector. However, a peak in the low frequency range might be present when high “noise” frequency is present in the signal (or image). At this point, the vector of extracted magnitudes is only indexed by the frequencies but has no associated frequencies. To calculate the corresponding frequencies, the FT simply takes (or generates) the index vector (1, 2, 3, …, length(magnitude vector)) and divides it by the length of the data block (in sec).
In 1D, the phases would represent a vector of the same length as the magnitude vector with the phases (0 to \(2\pi\) or \(-\pi\) to \(+\pi\)) of each frequency. Phase shifts are translations in space (e.g., x-axis) for a given wave component that are measured in angles (radians). For instance, shifting a wave \(f(x)=0.5sin(3wt)+0.25sin(10wt)\) by \(\frac{\pi}{2}\) would produce the following Fourier series:
\[f(t)=0.5sin\left (3wt+\frac{\pi}{2}\right )+0.25sin\left (10wt+\frac{\pi}{2}\right ).\]
In 2D, The Fourier transform (FT/IFT) for images is defined by: \[\hat{f}(u,v)=F(u,v)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{f(x,y)e^{-i2\pi(ux+vy)}dxdy}},\] \[f(x,y)=\hat{\hat{f}}(x,y)=\hat{F}(x,y)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{F(u,v)e^{i2\pi(ux+vy)}dudv}},\] where \(u\) and \(v\) are the spatial frequencies, \(F(u,v)=F_R(u,v)+iF_I(u,v)\) is a complex number for each pair of arguments, \[|F(u,v)|=\sqrt{F_R^2(u,v)+F_I^2(u,v)}\] is the magnitude spectrum, and \[arctan\left (\frac{F_I(u,v)}{F_R(u,v)}\right )\] is the phase angle spectrum.
The complex exponential \[e^{-i2\pi(ux+vy)}=cos(2\pi(ux+vy)) +i\ sin(2\pi(ux+vy))\] represents the real and imaginary (complex) sinusoidal terms in the 2D plane. The extrema of its real part (\(cos(2\pi(ux+vy))\)) occur at \(2\pi(ux+vy)=n\pi\). Using vector notation, \[2\pi(ux+vy)=2\pi \langle U, X \rangle =n\pi,\] where the extrema points \(U=(u,v)^T\) and \(X=(x,y)^T\) represent sets of equally spaced parallel lines with normal \(U\) and wavelength \(\frac{1}{\sqrt{u^2+v^2}}\).
Let’s define the index shifting paradigm associated with the discrete FT, which is simply used for convenience and better visualizaiton. It has no other relevance to the actual calculation of the FT and its inverse, IFT.
When applying the forward or reverse generalized discrete FT it is possible to shift the transform sampling in time and frequency domain by some real offset values, \(a,b\). Symbolically,
\[\hat{f}(k) = \sum_{n=0}^{N-1} f(n) e^{-\frac{2 \pi i}{N} (k+b) (n+a)} \quad \quad k = 0, \dots, N-1.\]
Note: Remember that in R, the array indices start with 1, not 0, as in some other languages.
The function fftshift() is useful for visualizing the Fourier transform with the zero-frequency component in the middle of the spectrum. Its inverse counterpart, ifftshift(), is needed to rearrange again the indices appropriately after the IFT is employed, so that the image is correctly reconstructed in spacetime. The FT only computes half of the frequency spectrum corresponding to the non-negative (positive and zero if the length(f) is odd) frequencies in order to save computation time. To preserve the dimensions of the output \(\hat{f}=FT(f)\), the second half of the frequency spectrum (the complex conjugate of the first half) is just added at the end of this vector. In a 1D setting, the result of fft() is:
\(0\ 1\ 2\ 3\ ...\ (freq\ bins > 0)\ ... {\frac{Fs}{2}}\) and \(-{\frac{Fs}{2}}\ ... \ (freq\ bins < 0)\ ...\ -3\ -2\ -1.\)
where \(F_s\) is the frequency sample. The fftshift method sets the zero-frequency component in the center of the array, i.e., it just shifts (offsets) the second part with the negative frequency bins to the beginning and the first part to the end of the resulting FT vector, or matrix. Thus, the shifted discrete FT can be nicely plotted in the center covering the frequency spectrum from \(-{\frac{Fs}{2}}\) on the left to \({\frac{Fs}{2}}\) on the right. This is not necessary, but is used for better visualization aesthetics. To synthesize back the correct image, after using fftshift on the FT signal, we always have to undo that re-indexing by using ifftshift() on the inverse-FT.
# FFT SHIFT
fftshift <- function(img_ff, dim = -1) {
rows <- dim(img_ff)[1]
cols <- dim(img_ff)[2]
# planes <- dim(img_ff)[3]
swap_up_down <- function(img_ff) {
rows_half <- ceiling(rows/2)
return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
}
swap_left_right <- function(img_ff) {
cols_half <- ceiling(cols/2)
return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
}
#swap_side2side <- function(img_ff) {
# planes_half <- ceiling(planes/2)
# return(cbind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], img_ff[1:rows, 1:cols, 1:planes_half]))
#}
if (dim == -1) {
img_ff <- swap_up_down(img_ff)
return(swap_left_right(img_ff))
}
else if (dim == 1) {
return(swap_up_down(img_ff))
}
else if (dim == 2) {
return(swap_left_right(img_ff))
}
else if (dim == 3) {
# Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
# install.packages("abind")
library(abind)
planes <- dim(img_ff)[3]
rows_half <- ceiling(rows/2)
cols_half <- ceiling(cols/2)
planes_half <- ceiling(planes/2)
img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)],
img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)],
img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)],
img_ff[1:rows, 1:cols, 1:planes_half], along=3)
return(img_ff)
}
else {
stop("Invalid dimension parameter")
}
}
ifftshift <- function(img_ff, dim = -1) {
rows <- dim(img_ff)[1]
cols <- dim(img_ff)[2]
swap_up_down <- function(img_ff) {
rows_half <- floor(rows/2)
return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
}
swap_left_right <- function(img_ff) {
cols_half <- floor(cols/2)
return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
}
if (dim == -1) {
img_ff <- swap_left_right(img_ff)
return(swap_up_down(img_ff))
}
else if (dim == 1) {
return(swap_up_down(img_ff))
}
else if (dim == 2) {
return(swap_left_right(img_ff))
}
else if (dim == 3) {
# Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
# install.packages("abind")
library(abind)
planes <- dim(img_ff)[3]
rows_half <- floor(rows/2)
cols_half <- floor(cols/2)
planes_half <- floor(planes/2)
img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)],
img_ff[1:rows, 1:cols, 1:planes_half], along=3)
img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)],
img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)],
img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
return(img_ff)
}
else {
stop("Invalid dimension parameter")
}
}2.1 FT and IFT of a pair of synthetic 2D images (a square and a disk)
# library(RCurl)
library(imager)
# install EBImage
# source("https://bioconductor.org/biocLite.R")
# biocLite("EBImage")
library(EBImage)
# Define two synthetic images
square_arr <- matrix(nrow=256, ncol=256)
circle_arr <- matrix(nrow=256, ncol=256)
for (i in 1:256) {
for (j in 1:256) {
if ( abs(i-128) < 30 && abs(j-128) < 30)
square_arr[i,j]=1 # sqrt((i-128)^2+(j-128)^2)/30
else square_arr[i,j]=0
if ( sqrt((i-128)^2 + (j-128)^2)<30)
circle_arr[i,j]=1 # 1-sqrt((i-128)^2+(j-128)^2)/30
else circle_arr[i,j]=0
}
}
# image(square_arr); image(circle_arr)
EBImage::display(square_arr, method = "raster"); display(circle_arr, method = "raster")
#download.file("http://wiki.socr.umich.edu/images/e/ea/BrainCortex2.png",
# paste(getwd(),"./image.png", sep="/"), mode = 'wb')
#im <- load.image(paste(getwd(),"./image.png", sep="/"))
# plot(im)
# Grayscaled
# img_gray <- im[ , , 1]
# display(img_gray, title='Image')
# FFT
# img_ff <- fft(img_gray) #fftw2d
ft_square <- fft(square_arr) #fftw2d # Display Re(FT): display(fftshift(ft_square))
ft_circle <- fft(circle_arr) # display(fftshift(ft_circle))
# Magnitude and Phase
mag_ft_square <- sqrt(Re(ft_square)^2+Im(ft_square)^2)
mag_ft_circle <- sqrt(Re(ft_circle)^2+Im(ft_circle)^2)
# Phase <- atan(Im(img_ff)/Re(img_ff))
phase_ft_square <- atan2(Im(ft_square), Re(ft_square))
phase_ft_circle <- atan2(Im(ft_circle), Re(ft_circle))
# FFT SHIFT
shift_ft_square <- fftshift(mag_ft_square)
shift_ft_circle <- fftshift(mag_ft_circle)
# Display FT
EBImage::display(log(shift_ft_square),title="FT Magnitude")
EBImage::display(log(shift_ft_circle),title="FT Phase")
# Magnitude and Phase
mag_shift_ft_square <- sqrt(Re(shift_ft_square)^2+Im(shift_ft_square)^2)
# phase <- atan(Im(img_ff)/Re(img_ff))
phase_ft_square <- atan2(Im(ft_square), Re(ft_square)); EBImage::display(fftshift(phase_ft_square))
phase_ft_circle <- atan2(Im(ft_circle), Re(ft_circle)); EBImage::display(fftshift(phase_ft_circle))
# ImplicitlyInvert the FT (IFT)
fftinv <- function( x ) { fft( x, inverse=TRUE ) / length( x ) }
EBImage::display(Re(fftinv(fft(square_arr))),title="(IFT o FT) Magnitude")
EBImage::display(Re(fftinv(fft(circle_arr))),title="(IFT o FT) Magnitude")
# FT of Square # Display the FT with a shift or not: display(fftshift(Re(X1))
X1 = fft(square_arr); display(fftshift(Re(X1)), method = "raster")
X1_mag <- sqrt(Re(X1)^2+Im(X1)^2); display(fftshift(X1_mag), method = "raster") # magnitude only
X1_phase <- atan2(Im(X1), Re(X1)); display(fftshift(X1_phase), method = "raster") # phase only
## Implicit Automated IFT
hat_X1 = Re(fft(X1, inverse = T)/length(square_arr)); display(hat_X1, method = "raster")
## Manually invert the FT (IFT) using the magnitudes and phases
Real1 = X1_mag * cos(X1_phase)
Imaginary1 = X1_mag * sin(X1_phase)
man_hat_X1 = Re(fft(Real1+1i*Imaginary1, inverse = T)/length(X1)); EBImage::display(man_hat_X1, method = "raster")
# FT of Circle # No shift applied here (perhaps should be consistent or just show the difference?)
X2 = fft(circle_arr); display(Re(X2), method = "raster")
X2_mag <- sqrt(Re(X2)^2+Im(X2)^2); display(X2_mag, method = "raster") # magnitude only
X2_phase <- atan2(Im(X2), Re(X2)); display(X2_phase, method = "raster") # phase only
## Implicit Automated IFT
# hat_X2 = Re(fft(X2, inverse = T)/length(circle_arr)); display(hat_X2, method = "raster")
## Manually invert the FT (IFT) using the magnitudes and phases
Real2 = X2_mag * cos(X2_phase)
Imaginary2 = X2_mag * sin(X2_phase)
man_hat_X2 = Re(fft(Real2+1i*Imaginary2, inverse = T)/length(X1)); EBImage::display(man_hat_X2, method = "raster")
# IFT Square-Magnitude and Circle-Phase
Real = X1_mag * cos(X2_phase)
Imaginary = X1_mag * sin(X2_phase)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
# IFT Circle-Magnitude and Square-Phase
Real = X2_mag * cos(X1_phase)
Imaginary = X2_mag * sin(X1_phase)
ift_X1phase_X2mag = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1phase_X2mag, method = "raster")
# IFT Circle-Magnitude and Nil-Phase
Real = X2_mag * cos(0)
Imaginary = X2_mag * sin(0)
ift_NilPhase_X2mag = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(X2)))
display(ift_NilPhase_X2mag, method = "raster")
# IFT Square-Magnitude and Nil-Phase
Real = X1_mag * cos(0)
Imaginary = X1_mag * sin(0)
ift_NilPhase_X1mag = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(X1)))
display(ift_NilPhase_X1mag, method = "raster")
# IFT Square-Magnitude and Phase = Avg(Square, Circle)
Real = X1_mag * cos((X2_phase+X1_phase)/2)
Imaginary = X1_mag * sin((X2_phase+X1_phase)/2)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
# IFT Magnitude == Avg(Square, Circle) and Phase = Avg(Square, Circle)
Real = ((X1_mag +X2_mag)/2)* cos((X2_phase+X1_phase)/2)
Imaginary = ((X1_mag +X2_mag)/2) * sin((X2_phase+X1_phase)/2)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
library(spatstat)
# IFT Square-Magnitude and Phase = Blur(Square, sigma=0.5)
Real = X1_mag * cos(as.matrix(blur(as.im(X1_phase), sigma=0.5)))
Imaginary = X1_mag * sin(as.matrix(blur(as.im(X1_phase), sigma=0.5)))
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
# IFT Magnitude= Blur(Square, sigma=0.5) and Phase = Square
Real = as.matrix(blur(as.im(X1_mag), sigma=3)) * cos(X1_phase)
Imaginary = as.matrix(blur(as.im(X1_mag), sigma=3)) * sin(X1_phase)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
### Note on phase transformations ##################
# Fourier surrogate realizations are used in time-series analysis and other applications to synthetically
# generate simulated data from the approximate type/distribution as some observed data
# these realizations may be constructed by manipulating the phases, e.g., by adding random noise to phases,
# deriving new highly correlated phases, or via other phase-transformations.
# Below we examine these phase-manipulation effects on the surrogate realizations
# using the 2D square and circle images; See Ref=https://doi.org/10.1016/S0167-2789(00)00043-9
# IFT Magnitude= Square and Phase = Square + IID noise
set.seed(1234)
IID_noise <- matrix(runif(prod(dim(X1_phase)), min = -1, max = 1), nrow = dim(X1_phase)[1]); dim(IID_noise) # 256 256## [1] 256 256# IID_noise <- matrix(rnorm(prod(dim(X1_phase)), mean=0, sd=1), nrow = dim(X1_phase)[1]); dim(IID_noise) # 256 256
# plot(density(IID_noise), xlim=c(-4,4), col="blue", lwd=2); lines(density(X1_phase), col="red", lwd=2)
Real = X1_mag * cos(X1_phase + IID_noise)
Imaginary = X1_mag * sin(X1_phase + IID_noise)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
# image(ift_X1mag_X2phase)
# IFT Magnitude= Circle and Phase = Square + IID noise (Uniform(-2,2))
set.seed(1234)
IID_noise <- matrix(runif(prod(dim(X1_phase)), min = -2, max = 2), nrow = dim(X1_phase)[1]); dim(IID_noise) # 256 256## [1] 256 256plot(density(IID_noise), xlim=c(-4,4), col="blue", lwd=2); lines(density(X1_phase), col="red", lwd=2)
Real = X2_mag * cos(X1_phase + IID_noise)
Imaginary = X2_mag * sin(X1_phase + IID_noise)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
windows(width=10, height=8)
plot(density(X1_phase + IID_noise), xlim=c(-5,5), col="blue", lwd=2,
main="Phase Distributions: Raw Square, Square+IID U(-2,2)",
xlab = "Phase", ylab = "Density")
lines(density(X1_phase), col="red", lwd=2)
legend("center", legend=c("Square + U(-2,2)", "(Raw) Square Phases"),
col=c("blue", "red"), lty=1, lwd=3, cex=1.0, y.intersp=1.0,
x.intersp=1.0, title = "Phases", bty = "n")
# randomly sample 10 indices to pairwise plot
ind1 <- sample(dim(X1_phase)[1], 10); ind2 <- sample(dim(X1_phase)[2], 10)
plot((X1_phase)[ind1,ind2], (X1_phase + IID_noise)[ind1,ind2],
main=sprintf("Square Image: Raw vs. IID Noise-corrupted Phases: Corr=%s",
round(cor(as.vector(X1_phase), as.vector(X1_phase + IID_noise)), digits=3)),
xlab = "Raw", ylab = "Phase + U(-2,2)")
plot(X1_mag, X2_mag, xlim=c(0,20), ylim=c(0, 20),
main="Fourier Magnitudes", xlab = "Square", ylab = "Circle")
# Square vs. Circle Phases
plot((X1_phase), (X2_phase),
main=sprintf("Square vs. Circle Phases: Corr=%s",
round(cor(as.vector(X1_phase), as.vector(X2_phase)), digits=3)),
xlab = "Square", ylab = "Circle")
# Scramble/Swap/Permute the Phases, IFT Magnitude= Square and Phase = permuted(Square)
set.seed(123)
X1_phase_swapped <-X1_phase[sample(nrow(X1_phase)), sample(ncol(X1_phase))] # Shuffle Phases row-wise & col-wise
# X1_phase_swapped <- matrix(sample(X1_phase), nrow = dim(X1_phase)[1], ncol = dim(X1_phase)[2])
# image(X1_phase); image(X1_phase_swapped) # compare the phase images (raw vs. purmuted phases)
Real = X1_mag * cos(X1_phase_swapped)
Imaginary = X1_mag * sin(X1_phase_swapped)
ift_X1mag_X1phase_purmuted = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X1phase_purmuted, method = "raster")
2.2 Compare the effects of the phases on reconstructing 2D images of the Earth and Saturn
# library(RCurl)
library(imager)
# install EBImage
# source("https://bioconductor.org/biocLite.R")
# biocLite("EBImage")
library(EBImage)
# Load Eath and Saturn 2D images
earth = readImage("E:/Ivo.dir/Eclipse_Projects/HTML5_WebSite/TCIU/Chapter6/Earth_Globe_2D_img.png")
saturn = readImage("E:/Ivo.dir/Eclipse_Projects/HTML5_WebSite/TCIU/Chapter6/Saturn_2D_img.png")
dim(earth); dim(saturn)## [1] 1142 1189 4## [1] 1143 1148 4EBImage::display(earth, method="raster")
EBImage::display(saturn, method="raster")
# text(x = 20, y = 20, label = "Earth", adj = c(0,1), col = "orange", cex = 2)
# resize to 256^2
# library(imager)
earth_256 = EBImage::resize(earth, w=256, h=256)
EBImage::display(earth_256, method="raster"); dim(earth_256)
## [1] 256 256 4saturn_256 = EBImage::resize(saturn, w=256, h=256)
EBImage::display(saturn_256[ , , 4], method="raster"); dim(saturn_256)
## [1] 256 256 4earth_arr <- imageData(earth_256[ , , 4]); dim(earth_arr); str(earth_arr)## [1] 256 256## num [1:256, 1:256] 0 0 0 0 0 0 0 0 0 0 ...EBImage::display(Image(earth_arr), method="raster")
saturn_arr <- imageData(saturn_256[ , , 4]); dim(saturn_arr); str(saturn_arr)## [1] 256 256## num [1:256, 1:256] 0.0283 0.0283 0.0283 0.0283 0.0283 ...EBImage::display(Image(saturn_arr), method="raster")
# FFT
# img_ff <- fft(img_gray) #fftw2d
ft_earth <- fft(earth_arr) #fftw2d # Display Re(FT): display(fftshift(ft_earth))
ft_saturn <- fft(saturn_arr) # display(fftshift(ft_saturn))
# Magnitude and Phase
mag_ft_earth <- sqrt(Re(ft_earth)^2+Im(ft_earth)^2)
mag_ft_saturn <- sqrt(Re(ft_saturn)^2+Im(ft_saturn)^2)
# Phase <- atan(Im(img_ff)/Re(img_ff))
phase_ft_earth <- atan2(Im(ft_earth), Re(ft_earth))
phase_ft_saturn <- atan2(Im(ft_saturn), Re(ft_saturn))
# FFT SHIFT
shift_ft_earth <- fftshift(mag_ft_earth)
shift_ft_saturn <- fftshift(mag_ft_saturn)
# Display FT
EBImage::display(log(shift_ft_earth),title="FT Magnitude (Earth)")
EBImage::display(log(shift_ft_saturn),title="FT Magnitude (Saturn)")
# Magnitude and Phase
mag_shift_ft_earth <- sqrt(Re(shift_ft_earth)^2+Im(shift_ft_earth)^2)
# phase <- atan(Im(img_ff)/Re(img_ff))
phase_ft_earth <- atan2(Im(ft_earth), Re(ft_earth)); EBImage::display(fftshift(phase_ft_earth))
phase_ft_saturn <- atan2(Im(ft_saturn), Re(ft_saturn)); EBImage::display(fftshift(phase_ft_saturn))
# ImplicitlyInvert the FT (IFT)
fftinv <- function( x ) { fft( x, inverse=TRUE ) / length( x ) }
display(Re(fftinv(fft(earth_arr))),title="(IFT o FT) Magnitude (Earth)")
display(Re(fftinv(fft(saturn_arr))),title="(IFT o FT) Magnitude (Saturn)")
# FT of Earth # Display the FT with a shift or not: display(fftshift(Re(X1))
X1 = fft(earth_arr); display(fftshift(Re(X1)), method = "raster")
X1_mag <- sqrt(Re(X1)^2+Im(X1)^2); display(fftshift(X1_mag), method = "raster") # magnitude only
X1_phase <- atan2(Im(X1), Re(X1)); display(fftshift(X1_phase), method = "raster") # phase only
## Implicit Automated IFT
hat_X1 = Re(fft(X1, inverse = T)/length(earth_arr)); display(hat_X1, method = "raster")
## Manually invert the FT (IFT) using the magnitudes and phases
Real1 = X1_mag * cos(X1_phase)
Imaginary1 = X1_mag * sin(X1_phase)
man_hat_X1 = Re(fft(Real1+1i*Imaginary1, inverse = T)/length(X1)); display(man_hat_X1, method = "raster")
# FT of Saturn # No shift applied here (perhaps should be consistent or just show the difference?)
X2 = fft(saturn_arr); display(Re(X2), method = "raster")
X2_mag <- sqrt(Re(X2)^2+Im(X2)^2); display(X2_mag, method = "raster") # magnitude only
X2_phase <- atan2(Im(X2), Re(X2)); display(X2_phase, method = "raster") # phase only
## Implicit Automated IFT
# hat_X2 = Re(fft(X2, inverse = T)/length(circle_arr)); display(hat_X2, method = "raster")
## Manually invert the FT (IFT) using the magnitudes and phases
Real2 = X2_mag * cos(X2_phase)
Imaginary2 = X2_mag * sin(X2_phase)
man_hat_X2 = Re(fft(Real2+1i*Imaginary2, inverse = T)/length(X1)); EBImage::display(man_hat_X2, method = "raster")
# IFT Earth-Magnitude and Saturn-Phase
Real = X1_mag * cos(X2_phase)
Imaginary = X1_mag * sin(X2_phase)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase, method = "raster")
# IFT Saturn-Magnitude and Earth-Phase
Real = X2_mag * cos(X1_phase)
Imaginary = X2_mag * sin(X1_phase)
ift_X1phase_X2mag = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
EBImage::display(ift_X1phase_X2mag, method = "raster")
# IFT Earth-Magnitude and Nil-Phase
Real = X1_mag * cos(0)
Imaginary = X1_mag * sin(0)
ift_NilPhase_X1mag = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(X1)))
EBImage::display(ift_NilPhase_X1mag, method = "raster")
# IFT Saturn-Magnitude and Nil-Phase
Real = X2_mag * cos(0)
Imaginary = X2_mag * sin(0)
ift_NilPhase_X2mag = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(X2)))
EBImage::display(ift_NilPhase_X2mag, method = "raster")
# IFT Earth-Magnitude and Phase = Avg(Earth, Saturn)
Real = X1_mag * cos((X2_phase+X1_phase)/2)
Imaginary = X1_mag * sin((X2_phase+X1_phase)/2)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
EBImage::display(ift_X1mag_X2phase, method = "raster")
# IFT Magnitude == Avg(Earth, Saturn) and Phase = Avg(Earth, Saturn)
Real = ((X1_mag +X2_mag)/2)* cos((X2_phase+X1_phase)/2)
Imaginary = ((X1_mag +X2_mag)/2) * sin((X2_phase+X1_phase)/2)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
EBImage::display(ift_X1mag_X2phase, method = "raster")
library(spatstat)
# IFT Earth-Magnitude and Phase = Blur(Earth, sigma=0.5)
Real = X1_mag * cos(as.matrix(blur(as.im(X1_phase), sigma=0.5)))
Imaginary = X1_mag * sin(as.matrix(blur(as.im(X1_phase), sigma=0.5)))
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
EBImage::display(ift_X1mag_X2phase, method = "raster")
# IFT Magnitude= Blur(Earth, sigma=0.5) and Phase = Earth
Real = as.matrix(blur(as.im(X1_mag), sigma=3)) * cos(X1_phase)
Imaginary = as.matrix(blur(as.im(X1_mag), sigma=3)) * sin(X1_phase)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
EBImage::display(ift_X1mag_X2phase, method = "raster")
### Note on phase transformations ##################
# Fourier surrogate realizations are used in time-series analysis and other applications to synthetically
# generate simulated data from the approximate type/distribution as some observed data
# these realizations may be constructed by manipulating the phases, e.g., by adding random noise to phases,
# deriving new highly correlated phases, or via other phase-transformations.
# Below we examine these phase-manipulation effects on the surrogate realizations
# using the 2D square and circle images; See Ref=https://doi.org/10.1016/S0167-2789(00)00043-9
# IFT Magnitude= Square and Phase = Square + IID noise
set.seed(1234)
IID_noise <- matrix(runif(prod(dim(X1_phase)), min = -1, max = 1), nrow = dim(X1_phase)[1]); dim(IID_noise) # 256 256## [1] 256 256# IID_noise <- matrix(rnorm(prod(dim(X1_phase)), mean=0, sd=1), nrow = dim(X1_phase)[1]); dim(IID_noise) # 256 256
# plot(density(IID_noise), xlim=c(-4,4), col="blue", lwd=2); lines(density(X1_phase), col="red", lwd=2)
Real = X1_mag * cos(X1_phase + IID_noise)
Imaginary = X1_mag * sin(X1_phase + IID_noise)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
EBImage::display(ift_X1mag_X2phase, method = "raster")
# image(ift_X1mag_X2phase)
# IFT Magnitude= Circle and Phase = Square + IID noise (Uniform(-2,2))
set.seed(1234)
IID_noise <- matrix(runif(prod(dim(X1_phase)), min = -2, max = 2), nrow = dim(X1_phase)[1]); dim(IID_noise) # 256 256## [1] 256 256plot(density(IID_noise), xlim=c(-4,4), col="blue", lwd=2); lines(density(X1_phase), col="red", lwd=2)
Real = X2_mag * cos(X1_phase + IID_noise)
Imaginary = X2_mag * sin(X1_phase + IID_noise)
ift_X1mag_X2phase = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
EBImage::display(ift_X1mag_X2phase, method = "raster")
windows(width=10, height=8)
plot(density(X1_phase + IID_noise), xlim=c(-5,5), col="blue", lwd=2,
main="Phase Distributions: Raw Square, Square+IID U(-2,2)",
xlab = "Phase", ylab = "Density")
lines(density(X1_phase), col="red", lwd=2)
legend("center", legend=c("Square + U(-2,2)", "(Raw) Square Phases"),
col=c("blue", "red"), lty=1, lwd=3, cex=1.0, y.intersp=1.0,
x.intersp=1.0, title = "Phases", bty = "n")
3 The effects of Kime-Magnitudes and Kime-Phases
Jointly, the amplitude spectrum (magnitudes) and the phase spectrum (phases) uniquely describe the spacetime representation of a signal. However, the importance of each of these two spectra is not equivalent. In general, the effect of the phase spectrum is more important compared to the corresponding effects of the amplitude spectrum. In other words, the magnitudes are less susceptible to noise or the accuracy of their estimations. The effects of magnitude perturbations are less critical relative to proportional changes in the phase spectrum. For instance, particularly in terms of spacetime locations where the signal is zero, the signal can be reconstructed (by the IFT) relatively accurately using incorrect magnitudes solely by using the correct phases REF. For a real valued signal \(f\), suppose the amplitude of its Fourier transform, \(FT(f)=\hat{f}\), is \(A(\omega) > 0, \forall \omega\), then: \[f(x)=IFT(\hat{f})=Re\left (\frac{1}{2\pi}\int_{R} \underbrace{A(\omega)e^{i\phi(\omega)}}_{\hat{f}(\omega)}\ e^{i\omega x}d\omega \right)= Re\left (\frac{1}{2\pi}\int_{R}A(\omega)e^{i(\phi(\omega)+\omega x)}d\omega\right) = \frac{1}{2\pi}\int_{R} {A(\omega) \cos(\phi(\omega)+\omega x)}d\omega.\]
Thus, the zeros of \(f(x)\) occur for \(\omega x+ \phi(\omega)=\pm k\frac{\pi}{2}\), \(k= 1,2,3,.\).
A solely amplitude driven reconstruction \(\left ( f_A(x)=IFT(\hat{f})=\frac{1}{2\pi}\int_{R}\underbrace{A(\omega)}_{no\ phase}\ e^{i\omega x}d\omega \right)\) would yield worse results than a solely-phase based reconstruction \(\left ( f_{\phi}(x)=IFT(\hat{f})=\frac{1}{2\pi} \int_{R}\underbrace{e^{i\phi(\omega)}}_{no\ amplitude}\ e^{i\omega x}d\omega\right )\). The latter would have a different total energy from the original signal, however, it would include some signal recognizable features as the zeroth-level curves of the original \(f\) and the phase-only reconstruction \(f_{\phi}\) signals will be preserved. This suggests that the Fourier phase of a signal is more informative than the Fourier amplitude, i.e., the magnitudes are robust to errors or perturbations.
In X-ray crystallography, crystal structures are bombarded by particles/waves, which are diffracted by the crystal to yield the observed diffraction spots or patterns. Each diffraction spot corresponds to a point in the reciprocal lattice and represents a particle wave with some specific amplitude and a relative phase. Probabilistically, as the particles (e.g., gamma-rays or photons) are reflected from the crystal, their scatter directions are proportional to the square of the wave amplitude, i.e., the square of the wave Fourier magnitude. X-rays capture these amplitudes as counts of particle directions, but miss all information about the relative phases of different diffraction patterns.
Spacekime analytics are analogous to X-ray crystallography, DNA helix modeling, and other applications, where only the Fourier magnitudes (time), i.e., power spectrum, is only observed, but not the phases (kime-directions), which need to be estimated to correctly reconstruct the intrinsic 3D object structure REF, in our case, the correct spacekime analytical inference. Clearly, signal reconstruction based solely on either the amplitudes or the phases is an ill-posed problem, i.e., there will be many alternative solutions. In practice, such signal or inference reconstructions are always application-specific, rely on some a priori knowledge on the process (or objective function), or depend an information-theoretic criteria to derive conditional solutions. Frequently, such solutions are obtained via least squares, maximum entropy criteria, maximum a posteriori distributions, Bayesian estimations, or simply by approximating the unknown amplitudes or phases using prior observations, similar processes, or theoretical models.
4 Solving the Missing Kime-Phase Problem
There are many alternative solutions to the problem of estimating the unobserved kime-phases. All solutions depend on the quality of the data (e.g., noise), the signal energy (e.g., strength of association between covariates and outcomes), and the general experimental design. There can be rather large errors in the phase reconstructions, which will in turn effect the final spacekime analytic results. Most phase-problem solutions are based on the idea that having some prior knowledge about the characteristics of the experimental design (case-study phenomenon) and the desired inference (spacekime analytics). For instance, if we artificially load the energy of the case-study (e.g., by lowering the noise, increasing the SNR, or increasing the strength of the relation between explanatory and outcome variables), the phases computed from the this stronger-signal dataset will be more accurate representations than the original phase estimates. Examples of phase-problem solutions include energy modification and fitting and refinement methods.
4.1 Energy Modification Strategies
In general, energy modification techniques rely on prior knowledge, testable hypotheses, or intuition to modify the dataset by strengthening the expected relation we are trying to uncover using spacekime analytics.
4.1.1 Kime-phase noise distribution flattening
In many practical applications, part of the dataset (including both cases and features) include valuable information, whereas the rest of the data may include irrelevant, noisy, or disruptive information.
Clearly, we can’t explicitly untangle these two components, however, we do expect that the irrelevant data portion would yield uninformative/unimportant kime-phases, which may be used to estimate the kime-phase noise-level and noise-distribution. Intuitively, if we modify the dataset to flatten the irrelevant kime-phases, the estimates of the corresponding true-signal kime-phases may be more accurate or more representative. We can think of this process as using kime-phase information from some known strong features to improve the kime-phase information of other particular features. Kime-phase noise distribution flattening requires that the kime-phases be good enough to detect the boundaries between the strong-features and the rest.
4.1.2 Multi-sample Kime-Phase Averaging
It’s natural to assume that multiple instances of the same process would yield similar analytics and inference results. For large datasets, we can use ensemble methods (e.g., SuperLearner, and CBDA) to iteratively generate independent samples, which would be expected to lead to analogous kime-phase estimated and analytical results. This, we expect that when salient features are extracted by spacekime analytics based on independent samples, their kime-phase estimates should be highly associated (e.g., correlated), albeit perhaps not identical. However, weak features would exhibit exactly the opposite effect - their kime-phases may be highly variable (noisy). By averaging the kime-phases, noisy-areas in the dataset may cancel out, whereas, patches of strong-signal may preserve the kime-phase details, which would lead to increased kime forecasting accuracy and reproducibility of the kime analytics.
4.1.3 Histogram equalization
As common experimental designs and similar datasets exhibit analogous characteristics, the corresponding spacekime analytics are also expected to be synergistic. Spacekime inference that does not yield results in some controlled or expected range, may be indicative of incorrect kime-phase estimation. We can use histogram equalization methods to improve the kime-phase estimates. This may be accomplished by altering the distribution of kime-phases to either match the phase distribution of other similar experimental designs or generate more expected spacekime analytical results.
4.1.4 Fitting and refinement
Related to energy modificaiton strategies, the fitting and refinement technique capitalizes on the fact that strong energy datasets tend to have a smaller set of salient features. So, if we construct case-studies with some strong features, the corresponding kime-phases will be more accurate, and the resulting inference/analytics will be more powerful and highly reproducible. Various classification, regression, supervised and unsupervised methods, and other model-based techniques allow us to both fit a model (estimate coefficients and structure) as well as apply the model for outcome predictions and forecasting. Such models permit control over the characteristics of individual features and multi-variate inter-relations, which can be can be exploited to gather valuable kime-phase information. Starting with a reasonable guess (kime-phase prior), the fitting and refinement technique can be applied iteratively to (1) reconstructing the data into spacetime using the kime-phase estimates, (2) fit or estimate the spacekime analytical model, (3) compare the analytical results and inference to expected outcomes, and (4) refine the kime-phase estimator aiming to gain better outcomes (#3). Indeed, other energy modificaiton strategies (e.g., averaging or flattening) can be applied before a new iteration to build a new model is initiated (#1 and #2).
# Take 2: IFT Magnitude= Square and Phase = Square + IID noise (N(0,3/2))
set.seed(1234)
IID_noise <- matrix(rnorm(prod(dim(X1_phase)), mean=0, sd=1.5), nrow=dim(X1_phase)[1]); dim(IID_noise) # 256 256## [1] 256 256plot(density(IID_noise), xlim=c(-8,8), col="blue", lwd=2); lines(density(X1_phase), col="red", lwd=2)
Real = X1_mag * cos(X1_phase + IID_noise)
Imaginary = X1_mag * sin(X1_phase + IID_noise)
ift_X1mag_X1phase_Noise = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X1phase_Noise, method = "raster")
windows(width=10, height=8)
mixed_density <- density(X1_phase + IID_noise)
mixed_density_mod2pi <- density((X1_phase + IID_noise)%%(2*pi) -pi)
plot(mixed_density, xlim=c(-8,8), ylim=c(0,0.17), col="blue", lwd=2, # should density be modulo %%(2*pi)?
main="Phase Distributions: Raw Square, Square+IID N(m=0,s=3/2), Mixed Phases mod 2*Pi",
xlab = "Phase", ylab = "Density")
lines(density(X1_phase), col="red", lwd=4)
lines(mixed_density_mod2pi, col="green", lwd=2)
text(x=3.2, y=-0.005, expression(pi)); text(x=-3.2, y=-0.005, expression(-pi))
legend("center", legend=c("(Raw) Square Phases", "Square + N(0,1.5)", expression(paste("Square + N(0,1.5) mod 2*", pi))),
col=c("red","blue", "green"), lty=1, lwd=c(4,2,2), cex=1.0, y.intersp=1.0,
x.intersp=1.0, title = "Phases", bty = "n")
# randomly sample 20 indices to pairwise plot
ind1 <- sample(dim(X1_phase)[1], 20); ind2 <- sample(dim(X1_phase)[2], 20)
corr <- cor.test((X1_phase)[ind1,ind2], (X1_phase + IID_noise)[ind1,ind2], method = "pearson", conf.level = 0.99)
plot((X1_phase)[ind1,ind2], (X1_phase + IID_noise)[ind1,ind2],
main=sprintf("Square Image: Raw vs. IID N(m=0,s=3/2) Noise-corrupted Phases: Corr=%s CI=(%s,%s)",
round(cor(as.vector(X1_phase[ind1,ind2]), as.vector(
(X1_phase + IID_noise)[ind1,ind2])), digits=3),
# round(cor(as.vector(X1_phase), as.vector(X1_phase + IID_noise)), digits=3),
round(corr$conf.int[1], digits=2),
round(corr$conf.int[2], digits=2)),
xlab = "Raw", ylab = "Phase + N(0,1.5)")
abline(lm(as.vector((X1_phase + IID_noise)[ind1,ind2]) ~ as.vector((X1_phase)[ind1,ind2])), col="red", lwd=2)
# Take 2: Same level of noise on Amplitudes:
# IFT Magnitude= Square + 50*IID noise (N(0,3/2)) and Phase = Square
set.seed(1234)
IID_noise <- matrix(rnorm(prod(dim(X1_mag)), mean=0, sd=3/2), nrow=dim(X1_mag)[1]); dim(IID_noise) # 256 256## [1] 256 256plot(density(X1_mag+50*IID_noise), xlim=c(0,40), ylim=c(0,1.5), col="blue", lwd=2)
lines(density(X1_mag), col="red", lwd=2)
# randomly sample 20 indices to pairwise plot
set.seed(1234); ind1 <- sample(dim(X1_phase)[1], 20); ind2 <- sample(dim(X1_phase)[2], 20)
corr <- cor.test((X1_mag+50*IID_noise)[ind1,ind2], (X1_mag)[ind1,ind2], method = "pearson", conf.level = 0.99)
plot(X1_mag[ind1,ind2], abs(X1_mag+50*IID_noise)[ind1,ind2], xlim=c(0,10), ylim=c(0,200),
main=sprintf("Square Image: Mag + 50*IID N(m=0,s=3/2) Noise-corrupted vs. Raw Amplitudes: Corr=%s CI=(%s,%s)",
round(cor(as.vector(X1_mag[ind1,ind2]), as.vector(
(X1_mag + 50*IID_noise)[ind1,ind2])), digits=3),
round(corr$conf.int[1], digits=2),
round(corr$conf.int[2], digits=2)),
xlab = "Raw", ylab = "Amplitude + 50*N(0,1.5)")
abline(lm(as.vector(abs(X1_mag+50*IID_noise)[ind1,ind2]) ~ as.vector((X1_mag)[ind1,ind2])), col="red", lwd=2)
Real = (X1_mag+ 50*IID_noise) * cos(X1_phase)
Imaginary = (X1_mag+ 50*IID_noise) * sin(X1_phase)
ift_X1mag_X1phase_Noise = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X1phase_Noise, method = "raster")
# Magnitude distributions
mixed_density <- density(abs(X1_mag + 50*IID_noise))
windows(width=6, height=5)
plot(mixed_density, xlim=c(0,200), col="blue", lwd=2,
main="Magnitude Distributions: Raw Square, Square+50*N(m=0,s=3/2)",
xlab = "Magnitude", ylab = "Density")
lines(density(X1_mag), col="red", lwd=4)
legend("topright", legend=c("(Raw) Square Magnitudes", "Square + 50*N(0,1.5)"),
col=c("red","blue"), lty=1, lwd=c(4,2), cex=1.0, y.intersp=1.0,
x.intersp=1.0, title = "Magnitudes", bty = "n")
#### Take 3: Linear Transform of the Phases: SquarePhase ~ CirclePhase
# Take 2: IFT Magnitude= Square and Phase = LM(CirclePhase)
lm_Squ_Cir <- lm(as.vector(X1_phase) ~ as.vector(X2_phase))
plot(as.vector(X2_phase), as.vector(X1_phase), col="blue", lwd=2, xlab = "Circle", ylab = "Square",
main=sprintf("Linear Phase Transformation (SquarePhase ~ CirclePhase), Corr(Cir, Squ)=%s",
round(cor(as.vector(X1_phase) , as.vector(X2_phase)), digits=3)))
abline(lm_Squ_Cir, col="red", lwd=2)
Real = X1_mag * cos(lm_Squ_Cir$coefficients[1] + lm_Squ_Cir$coefficients[2]*X2_phase)
Imaginary = X1_mag * sin(lm_Squ_Cir$coefficients[1] + lm_Squ_Cir$coefficients[2]*X2_phase)
ift_X1mag_X2phase_LM = Re(fft(Real+1i*Imaginary, inverse = T)/length(X1))
display(ift_X1mag_X2phase_LM, method = "raster")
5 2D FT and Image Synthesis
Let’s try to perform the same image synthesis to reconstruct the Cyrillic and English alphabet images.
library(EBImage)
library(imager)
img_CyrAlpha <- load.image("E:/Ivo.dir/Research/UMichigan/Publications_Books/2019/DataScience_Book_Value_Uncertainty_Kime_2019/other/BulgAlpha.jpg")
# plot(img_CyrAlpha, axes = F)
# Grayscaled # img_gray <- im[ , , 1, 1]
img_CyrAlpha <- matrix(img_CyrAlpha, nrow = dim(img_CyrAlpha)[1], ncol = dim(img_CyrAlpha)[2])
# EBImage::display(img_CyrAlpha, title='Image', method = "raster")
img_EngAlpha <- load.image("E:/Ivo.dir/Research/UMichigan/Publications_Books/2019/DataScience_Book_Value_Uncertainty_Kime_2019/other/EngAlpha.jpg")
img_EngAlpha <- matrix(img_EngAlpha, nrow = dim(img_EngAlpha)[1], ncol = dim(img_EngAlpha)[2])
# add an extra zero column at the end to make the 2D Alphabet images homologous: 411 353
img_EngAlpha <- cbind(img_EngAlpha, rep(0, dim(img_EngAlpha)[1]))
dim(img_CyrAlpha); dim(img_EngAlpha)## [1] 411 353## [1] 411 353# FFT
ft_img_CyrAlpha <- fft(img_CyrAlpha) # fftw2d # Display Re(FT): display(fftshift(ft_img_CyrAlpha))
ft_img_EngAlpha <- fft(img_EngAlpha) # display(fftshift(ft_img_EngAlpha))
# Magnitude and Phase
mag_ft_img_CyrAlpha <- sqrt(Re(ft_img_CyrAlpha)^2+Im(ft_img_CyrAlpha)^2)
mag_ft_img_EngAlpha <- sqrt(Re(ft_img_EngAlpha)^2+Im(ft_img_EngAlpha)^2)
# Phase <- atan(Im(img_ff)/Re(img_ff))
phase_ft_img_CyrAlpha <- atan2(Im(ft_img_CyrAlpha), Re(ft_img_CyrAlpha))
phase_ft_img_EngAlpha <- atan2(Im(ft_img_EngAlpha), Re(ft_img_EngAlpha))
# FFT SHIFT
shift_ft_img_CyrAlpha <- fftshift(mag_ft_img_CyrAlpha)
shift_ft_img_EngAlpha <- fftshift(mag_ft_img_EngAlpha)
# Display FT
EBImage::display(log(shift_ft_img_CyrAlpha), title="FT Magnitude (Cyrillic Alphabet)")
EBImage::display(log(shift_ft_img_EngAlpha), title="FT Phase (English Alphabet)")
# ImplicitlyInvert the FT (IFT)
fftinv <- function( x ) { fft( x, inverse=TRUE ) / length( x ) }
EBImage::display(Re(fftinv(ft_img_CyrAlpha)),title="(IFT o FT) Magnitude (Cyrillic Alphabet)")
EBImage::display(Re(fftinv(ft_img_EngAlpha)),title="(IFT o FT) Magnitude (English Alphabet)")
############## FT of img_CyrAlpha
#X1 = fft(img_CyrAlpha); display(fftshift(Re(X1)), method = "raster")
#X1_mag <- sqrt(Re(X1)^2+Im(X1)^2); display(fftshift(X1_mag), method = "raster") # magnitude only
#X1_phase <- atan2(Im(X1), Re(X1)); display(fftshift(X1_phase), method = "raster") # phase only
##### Implicit Automated IFT
#hat_X1 = Re(fft(X1, inverse = T)/length(square_arr)); display(hat_X1, method = "raster")
###### Manually invert the FT (IFT) using the magnitudes and phases
#Real1 = X1_mag * cos(X1_phase)
#Imaginary1 = X1_mag * sin(X1_phase)
#man_hat_X1 = Re(fft(Real1+1i*Imaginary1, inverse = T)/length(X1)); display(man_hat_X1, method = "raster")
############### FT of img_EngAlpha
#X2 = fft(circle_arr); display(Re(X2), method = "raster")
#X2_mag <- sqrt(Re(X2)^2+Im(X2)^2); display(X2_mag, method = "raster") # magnitude only
#X2_phase <- atan2(Im(X2), Re(X2)); display(X2_phase, method = "raster") # phase only
###### Implicit Automated IFT
# hat_X2 = Re(fft(X2, inverse = T)/length(circle_arr)); display(hat_X2, method = "raster")
###### Manually invert the FT (IFT) using the magnitudes and phases
#Real2 = X2_mag * cos(X2_phase)
#Imaginary2 = X2_mag * sin(X2_phase)
#man_hat_X2 = Re(fft(Real2+1i*Imaginary2, inverse = T)/length(X1)); display(man_hat_X2, method = "raster")
# IFT Magnitude=mag_ft_img_CyrAlpha AND Phase=phase_ft_img_EngAlpha
Real = mag_ft_img_CyrAlpha * cos(phase_ft_img_EngAlpha)
Imaginary = mag_ft_img_CyrAlpha * sin(phase_ft_img_EngAlpha)
ift_MagCyr_PhaseEng = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_img_CyrAlpha))
EBImage::display(ift_MagCyr_PhaseEng, method = "raster", title="(IFT o FT) Magnitude=Cyr | Phase=Eng")
# IFT Magnitude=mag_ft_img_EngAlpha and Phase=phase_ft_img_CyrAlpha
Real = mag_ft_img_EngAlpha * cos(phase_ft_img_CyrAlpha)
Imaginary = mag_ft_img_EngAlpha * sin(phase_ft_img_CyrAlpha)
ift_MagEng_PhaseCyr = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_img_CyrAlpha))
EBImage::display(ift_MagEng_PhaseCyr, method = "raster", title="(IFT o FT) Magnitude=Eng | Phase=Cyr")
# IFTMagnitude=mag_ft_img_CyrAlpha and Phase=Nil
Real = mag_ft_img_CyrAlpha * cos(0)
Imaginary = mag_ft_img_CyrAlpha * sin(0)
ift_MagCyr_PhaseNil = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_img_CyrAlpha)))
EBImage::display(ift_MagCyr_PhaseNil, method = "raster")
# IFT Magnitude=mag_ft_img_EngAlpha and Phase=Nil
Real = mag_ft_img_EngAlpha * cos(0)
Imaginary = mag_ft_img_EngAlpha * sin(0)
ift_MagEng_PhaseNil = Re(ifftshift(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_img_CyrAlpha)))
EBImage::display(ift_MagEng_PhaseNil, method = "raster")