T02: Ajuste de curvas paramétricas

[AntonioG]

Code: Select all

x<-runif(20,3,7)
y<-3*x-2+rnorm(20)
plot(x,y,type="p",pch=19)
xy<-cbind(x,y)

Code: Select all

# Minimos cuadrados
M<-lm(y~x)
summary(M)
abline(M,lwd=2,col="red")

# A "mano"

Code: Select all

f<-function(B,xy){
   B0<-B[1]
   B1<-B[2]
   x<-xy[,1]
   y<-xy[,2]
   e<-y-(B0+B1*x)
   r<-sum(e^2)
   return(r)
   }

res<-optim(par=c(0,1),fn=f,xy=xy)
str(res)
res$par
abline(res$par,lwd=2,col="red")

# De aqui en adelante es solamente para visualizar la funcion

Code: Select all

fxy <- function (B0, B1, xy){
   B<-cbind(B0,B1)
   r<-apply(B,1,f,xy=xy)
   return(r)
   }

B0i <- seq(-10, 10, length= 200)
B1i <- seq(-10, 10, length=200)
fi <- outer(B0i, B1i, fxy,xy=xy)

persp(B0i, B1i, fi, theta = 45, phi = 0, expand = 0.5, col = "lightblue")
contour(B0i, B1i, fi,nlevels = 100)

minij<-which(fi==min(fi),arr.ind = TRUE)
B0i[minij[1]]
B1i[minij[2]]

[AntonioG]

Code: Select all

library(jpeg)

plotIm<-function(x){
   x<-(x-min(x))/(max(x)-min(x))
   nr<-dim(x)[1]
   nc<-dim(x)[2]
   par(mai=c(0,0,0,0))
   plot(c(1,nc),c(1,nr), type='n', axes = FALSE, xlab = "", ylab = "")
   rasterImage(x, 0, 0, nc, nr,interpolate=FALSE)
   }

ruta<-file.choose()
img <- readJPEG(ruta)
R <- img[,,1]
plotIm(R)
[code]

[code]
nr<-dim(R)[1]
nc<-dim(R)[2]

dr<-diff(R)
dc<-t(diff(t(R)))
par(mfcol=c(1,2))
plotIm(dr)
plotIm(dc)

dr2<-diff(R, differences = 2)
dc2<-t(diff(t(R), differences = 2))
par(mfcol=c(1,2))
plotIm(dr2)
plotIm(dc2)

str(dr2)
str(dc2)
Drc<-dr2[,1:(nc-2)]^2+dc2[1:(nr-2),]^2
plotIm(Drc)
[code]

[code]
det<-Drc>0.015
plotIm(det)

xy<-which(det==1,arr.ind = TRUE)

Code: Select all

f<-function(par,xy){
   x0<-par[1]
   y0<-par[2]
   r<-par[3]
   x<-xy[,1]
   y<-xy[,2]
   e<-((x-x0)^2+(y-y0)^2)-r^2
   res<-sum(e^2)
   return(res)
   }

res<-optim(par=c(nc/2,nr/2,1),fn=f,xy=xy)
str(res)
res$par

Code: Select all

plotcircle<-function(x0,y0,r){
   theta=seq(0, 2*pi, length=200)
   x=r*cos(theta)+x0
   y=r*sin(theta)+y0
   lines(x, y, lwd=4, col="red")
   }
plotIm(R)
plotcircle(res$par[2],nr-res$par[1],res$par[3])