Skip to content
Permalink
master
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
156 lines (129 sloc) 3.49 KB
#
# above-below-light.R, 30 Oct 30
# Data from:
# Example
#
# Example from:
# Evidence-based Software Engineering: based on the publicly available data
# Derek M. Jones
#
# TAG example_optical-illusion
source("ESEUR_config.r")
library("plyr")
par(mar=MAR_default-c(1.0, 2.0, 0.7, 0.7))
# # step size of 0.005 and 60 shades of grey produce acceptable size of pdf
# To reduce the size of the image use (does not seem to impact the effect):
# step size of 0.010 and 70 shades of grey produce acceptable size of pdf
shade_col=sequential_hcl(70, c=0, l=c(30, 90), power=2.2)
#shade_col=rev(shade_col)
# Got the algorithm below from rosetta code.
#
# The following ideas produces images that did not look realistic.
# Calculate distance from a point y_offset from the
# circle, i.e., a light source.
# Simulate a boundary using a logistic equation, numbers
# fitted so boundary occurs roughly in the middle of the circle.
ld_boundary=function(x_vals, y_vals)
{
dist=sqrt(1.001-x_vals^2-y_vals^2)
#dist=sqrt(x_vals^2+(y_offset+y_vals)^2)
col=shade_col[22+30*(dist-(x_vals+y_vals)/2.0)]
if (length(y_vals) != length(col)) print(x_vals[1])
return(col)
}
# for (X in c(seq(-1, 1, 0.2), 0.387, 0.445))
# {
# y_vals=seq(-sqrt(1-X^2), sqrt(1-X^2), by=0.001)
# x_vals=rep(X, length(y_vals))
# dist=sqrt(1.0001-x_vals^2-y_vals^2)
# c=22+30*(dist-(x_vals+y_vals)/2.0)
# print(c(min(c), max(c), which(c < 1)))
# }
# plot(c)
#ld=ld_boundary(x_vals, y_vals)
mk_pt=function(X)
{
y_vals=seq(-sqrt(1-X^2), sqrt(1-X^2), by=0.010)
x_vals=rep(X, length(y_vals))
return(data.frame(x=x_vals,
y=y_vals,
col=ld_boundary(x_vals, y_vals)))
}
plot_light_above=function(x_offset, y_offset=3.2)
{
points(x_offset+shade_pts$x, y_offset-shade_pts$y, col=shade_pts$col, pch=".")
}
plot_light_below=function(x_offset, y_offset=2.2)
{
# Reverse x-points so light source comes from the same direction
# as plot_light_above
points(x_offset+rev(shade_pts$x), y_offset+shade_pts$y, col=shade_pts$col, pch=".")
}
x_vals=seq(-1, 1, by=0.010)
shade_pts=adply(x_vals, 1, mk_pt)
shade_pts$col=as.character(shade_pts$col)
plot(0, type="n", axes=FALSE,
xlim=c(0, 9.2), ylim=c(0, 9.2),
xlab="", ylab="")
polygon(c(0, 20, 20, 0, 0),
c(0, 0, 20, 20, 0), col=shade_col[25], border=NA)
plot_light_above(1.7, 7.0)
plot_light_above(4.7, 7.0)
plot_light_above(7.7, 7.0)
# plot_light_above(10.7, 7.0)
plot_light_below(1.7)
plot_light_below(4.7)
plot_light_below(7.7)
# plot_light_below(10.7)
#
# Light shining on a sphere.
# Algorithm below from rosetta code.
#
# dot=function(x, y, z)
# {
# d=light[1]*x + light[2]*y + light[3]*z
# return(ifelse(d < 0, -d, 0))
# }
#
#
# draw_line=function(line)
# {
# x_pos=line
# y=x
# y_2=y*y
#
# z=sqrt(r_2 - line*line - y_2)
#
# # Normalise the vector
# len = sqrt(line*line+y_2+z*z)
# line=line/len
# y=y/len
# z=z/len
# b = dot(line, y, z)^k + ambient
# intensity = round((1 - b)*length(shades))
# intensity[intensity < 0]=NA
# intensity[intensity == 0]=1 # black
# intensity[intensity > length(shades)]=length(shades) # white
# # print(length(shades[intensity]))
# points(rep(x_pos, length(x)), x, col=shades[intensity], pch=".");
# }
#
#
# shades=sequential_hcl(50, c=0, l=c(30, 90), power=1.2)
# radius=400
# r_2=radius^2
# ambient=0.00
# k=2
# x=floor(-radius):ceiling(radius)
#
# light=c(00, -120, -20)
# light=c(00, 120, -20)
# light=light/sqrt(sum(light*light)) # normalise
#
# plot(0, type="n", axes=FALSE,
# xlim=c(-radius, radius), ylim=c(-radius, radius),
# xlab="", ylab="")
#
# im=a_ply(as.array(x), 1, draw_line)
#
#