########################### ## FRACTAL INTERPOLATION ## ########################### #Points, X- and Y-coordinates. X=c(0,20,60,95,110) Y=c(0,60,30,10,20) # Coefficients for affine transformation. A=c() B=rep(0,(length(Y)-1)) C=c() # D is for vertical scaling, free parameter, 0 < D < 1. D=c(0.3,0.2,0.4,0.1) # Check the number of D's: n-1. E=c() F=c() N=length(X) for(n in 1:(length(X)-1)) { A[n]=(X[n+1]-X[n])/(X[N]-X[1]) E[n]=((X[N]*X[n])-(X[1]*X[n+1]))/(X[N]-X[1]) C[n]=(Y[n+1]-Y[n])/(X[N]-X[1])-(D[n]*((Y[N]-Y[1])/(X[N]-X[1]))) F[n]=((X[N]*Y[n]-X[1]*Y[n+1])/(X[N]-X[1]))-(D[n]*((X[N]*Y[1]-X[1]*Y[N])/(X[N]-X[1]))) } BMat=cbind(A,B,C,D,E,F) #-----------------------------------# x=0;y=0;x1=c();y1=c() # This loop is for finding the limits of the picture: for(i in 1:3000) { x2=x i=sample(1:length(BMat[,1]),1) x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5] y=BMat[i,3]*x+BMat[i,4]*y+BMat[i,6] x1=c(x1,x);y1=c(y1,y) } X11() plot(0,0,xlim=c(min(x1),max(x1)),ylim=c(min(y1),max(y1)+10),col="white",xlab="x", ylab="f(x)") points(X,Y,pch=20,type="b") for(j in 1:40000) { x3=x i=sample(1:length(BMat[,1]),1) x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5] y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6] points(x,y,pch=".",cex=1) } title(main="Fractal Interpolation")