The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
SOCR ≫ | TCIU Website ≫ | TCIU GitHub ≫ |
The Laplace transform allows us to examine the relations between the space-time and space-kime representations of longitudinal data. The Fourier transformation is a linear operator that maps complex-valued functions of real variables (e.g., space, time domains) to complex valued functions of other real variables (e.g., frequency domain). The Laplace transform is similar. However, it sends complex-valued functions of positive real variables (e.g., time) to complex-valued functions defined on complex variables (e.g., kime).
In this vignette, we introduce the use of inv_kimesurface_transform() and kimesurface_transform() functions in our package. We demonstrate how to do Laplace transform and kimesurface transform with our functions.
require(TCIU)
require(doParallel)
require(cubature)
require(oro.nifti)
require(magrittr)
require(plotly)
require(ggplot2)
Here we first apply the discrete Laplace Transform (LT) function in our package on the sine function with the domain of \([0 : 2\pi]\) to see whether it has the same 2D surface function as the analytic form of LT of sine. The analytic form is \(1/(1+z^2)\).
# For this part of code, we comment it out and import the output plot already generated before to reduce time
# But this part of code can be run successfully. If you are interested, you can try it on your computer!
# discrete Laplace Transform of sine
= 2
range_limit = seq(from = 0, to = range_limit, length.out = 50)[2:50]
x2 # drop the first row to avoid real part value of 0
= seq(from = 0, to = range_limit, length.out = 50)[2:50]
y2 # drop the first column to avoid imaginary part value of 0
# Recompute the LT(sin) discretized on lower-res grid
= array(dim=c(length(x2), length(y2)))# x2 %o% y2
z2_grid
<- function(t) { sin(t) }
f_sin
# kime surface transform
# use parallel computing to speed up code
= 2 # please choose the ncors according to the number of cores your PC has
ncors # it is better that you increase the number of cores used for parallel computing if your computer allows
<- makeCluster(ncors)
cl registerDoParallel(cl)
= list()
F for (i in 1:length(x2) ){
=
F[[i]] foreach(j = 1:length(y2),
.export='cubintegrate',
.packages='cubature') %dopar% {
::LT(FUNCT=f_sin, complex(real=x2[i], imaginary = y2[j]))
TCIU
}
}
stopCluster(cl)
= lapply(F, unlist)
F_vec = unlist(do.call(rbind, F_vec))
z2_grid
# explicit form of Laplace Transform of sine
= function(p) { 1/(p^2 + 1) } # Exact Laplace transform of sin(x), continuous function
laplace_sine
= expand.grid(X=x2, Y=y2)
XY = mapply(complex, real=XY$X,imaginary=XY$Y)
complex_xy =laplace_sine(complex_xy)
sine_z dim(sine_z) = c(length(x2), length(y2))# dim(sine_z) # [1] 49 49
# make the two plots in the same plot to compare
= plot_ly(hoverinfo="none", showscale = FALSE)%>%
lt_ilt_plotly add_trace(z=Re(sine_z)-1, type="surface", surfacecolor=Im(sine_z)) %>%
add_trace(z = Re(z2_grid), type="surface", opacity=0.7, surfacecolor=Im(z2_grid) )%>%
layout(title =
"Laplace Transform, LT(sin()), Height=Re(LT(sin())), Color=Re(LT(sin())) \n Contrast Exact (Continuous) vs.
Approximate (Discrete) Laplace Transform", showlegend = FALSE)
lt_ilt_plotly
1]] sample_save[[
From the plot, we can easily tell that our discrete LT function generate the same 2D surface as the analytic LT function does.
Here we first apply the discrete Laplace Transform (LT) function in our package on the sine function with the domain of \([0, 2\pi]\), and then use the inverse Laplace Transform (ILT) function in our package to prove the Laplace Transformation has been converted back to sine.
# For this part of code, we comment it out and import the output plot already generated before to reduce time
# But this part of code can be run successfully. If you are interested, you can try it on your computer!
# discrete Laplace Transform of sine
= function(t) { sin(t) }
f_sin = function(z) TCIU::LT(f_sin, z)
lt_sine
# inverse Laplace Transform on the lt_sine
# using parallel computing to speed up code
<- seq(0, pi*2, length.out = 20)
tvalsn <- makeCluster(ncors)
cl registerDoParallel(cl)
<- foreach(t=1:length(tvalsn),
sinvalsn .export='cubintegrate',
.packages='cubature') %dopar% {
::ILT(FUNCT=lt_sine, t=tvalsn[t])
TCIU
}stopCluster(cl)
= unlist(sinvalsn)
sinvalsn
# make the plot of the result from ILT
# to see whether it still looks like sine
<- as.data.frame(cbind(Re=Re(sinvalsn),Im=Im(sinvalsn),
sinvalsn_df2 Sin=sin(tvalsn), time_points=tvalsn))
= ggplot(sinvalsn_df2, aes(x=time_points))+
lt_ilt_sine geom_line(aes(y=Re, color="Real"), linetype=1, lwd=2) +
geom_line(aes(y = Sin, color="Sin"), linetype=2, lwd=1) +
scale_color_manual(name="Index",
values = c("Real"="steelblue", "Sin"="darkred"))+
labs(title = "Original fMRI Time-series f(t)=sin(t) and \n Reconstructed f'(t)=ILT(F)=ILT(discrete LT(f))",
subtitle = bquote("F" ~ "=" ~ "discrete LT(sine)")) +
xlab("Time") + ylab("fMRI Image Intensities (f and f')") +
theme_grey(base_size = 16) +
theme(legend.title = element_text(size=14, color = "black", face="bold"),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank(),
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
lt_ilt_sine
2]] sample_save[[
In this part, we still apply the Kimesurface transform on our sine function, and then check whether the inverse Kimesurface transform can convert it back to our original form of sine. Kimesurface transform mainly does a Laplace Transform and then applies a height and surface transformation to have a better view of the plot.
Below, we apply the Kimesurface transform on the sine function, and make the plot of the 2D function using the magnitude of the complex value as y.
= seq(0, 2, length.out=50)[2:50]; y = seq(0, 2, length.out=50)[2:50];
x # do Kimesurface transform on sine function
= kimesurface_transform(FUNCT = function(t) { sin(t) },
z2_grid real_x = x, img_y = y)
# make the plot after Kimesurface transformation
<- atan2(Im(z2_grid), Re(z2_grid))
surf_color = cbind(seq(0, 1, by=1/(length(x) - 1)), rainbow(length(x)))
colorscale <- (sqrt( Re(z2_grid)^2+ Im(z2_grid)^2))
magnitude
<- plot_ly(hoverinfo="none", showscale = FALSE) %>%
p add_trace(z = magnitude,
surfacecolor=surf_color, colorscale=colorscale, #Phase-based color
type = 'surface', opacity=1, visible=T) %>%
layout(title = "fMRI Kime-Surface, F=LT(fMRI) \n Height=Mag(F), Color=Phase(F)", showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0) ) ) # 1:1:1 aspect ratio
p
3]] sample_save[[
Also, you can try to apply the Kimesurface transformation on our fMRI data for a brain. Each voxel of the brain contains a set of time series data that can be analyzed using our transform. However, due to the number of cores limited for parallel computing in the vignette, we did not run this chuck. You can try it if you are interested.
# load the fMRI data
<- "http://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
fMRIURL <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
fMRIFile download.file(fMRIURL, dest=fMRIFile, quiet=TRUE)
<- readNIfTI(fMRIFile, reorient=FALSE)
fMRIVolume # dimensions: 64 x 64 x 21 x 180 ; 4mm x 4mm x 6mm x 3 sec
# extract a set of time series data from a voxel of fMRI data
<- fMRIVolume[20, 20, 11, ]
xA_fMRI_1D_x20_y20_z11
# get the smooth function of this time series data
# Instead of using the extremely noisy fMRI data and avoiding integration problems,
# smooth "f" and use the **smooth version, f**
<- seq(0+0.001, 2*pi, length.out = 180)
time_points <- smooth.spline(ceiling((180*time_points)/(2*pi)),
f df = 10)$y
xA_fMRI_1D_x20_y20_z11,
# Define the f(t)=smooth(fMRI)(t) signal as a function of real time 0<t<=2*pi
<- function(t) {
fmri_funct if (t < 0+0.001 || t > 2*pi) { return ( 0 ) } else {
return ( f[ceiling((180*t)/(2*pi))] )
}
}
# do the Kimesurface transform
= 8
ncors # please choose the ncors according to the number of cores your PC has
= seq(0, 2, length.out=50)[2:50]; y = seq(0, 2, length.out=50)[2:50];
x # do Kimesurface transform on sine function
= kimesurface_transform(FUNCT = fmri_funct,
z2_grid_fmri glb_para="f",
real_x = x, img_y = y,
parallel_computing = TRUE,
ncor=ncors)
# make the plot of function after Kimesurface transformation
<- atan2(Im(z2_grid_fmri), Re(z2_grid_fmri))
surf_color = cbind(seq(0, 1, by=1/(length(x) - 1)), rainbow(length(x)))
colorscale <- (sqrt( Re(z2_grid_fmri)^2+ Im(z2_grid_fmri)^2))
magnitude
<- plot_ly(hoverinfo="none", showscale = FALSE) %>%
p_fmri add_trace(z = magnitude,
surfacecolor=surf_color, colorscale=colorscale, # Phase-based color
type = 'surface', opacity=1, visible=T) %>%
layout(title = "fMRI Kime-Surface, F=LT(fMRI) \n Height=Mag(F), Color=Phase(F)", showlegend = FALSE,
scene = list(aspectmode = "manual", aspectratio = list(x=1, y=1, z=1.0) ) ) # 1:1:1 aspect ratio
p_fmri
4]] sample_save[[
After the Kimesurface transformation, we apply the inverse Kimesurface transformation, and see that we can get a curve that captures the most trend of the sine function.
<- seq(0+0.001, 2*pi, length.out = 180)
time_points = inv_kimesurface_transform(time_points, z2_grid)
inv_data = inv_kimesurface_transform(time_points, z2_grid,num_length = 23,
inv_data m=1, msg=TRUE)
<- as.data.frame(cbind(Re=scale(Re(inv_data$Smooth_Reconstruction)),
time_Intensities_ILT_df2 Im=scale(Re(inv_data$Raw_Reconstruction)),
fMRI=scale(Re(sin(time_points))),
time_points=time_points))
colnames(time_Intensities_ILT_df2) = c("Smooth Reconstruction",
"Raw Reconstruction",
"Original Sin", "time_points")
= reshape2::melt(time_Intensities_ILT_df2, id.var = "time_points")
df <-ggplot(df, aes(x = time_points, y = value, colour = variable)) +
ppppgeom_line(linetype=1, lwd=3) +
ylab("Function Intensity") + xlab("Time") +
theme(legend.position="top")+
labs(title= bquote("Comparison between" ~ "f(t)=Smooth(Sin)(t)" ~ "and Smooth(ILT(LT(Sin)))(t); Range [" ~ 0 ~":"~ 2*pi~"]"))
5]] sample_save[[
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.