distance<-function(x,y,m,q) { abs((y-m*x-q)/sqrt(m^2+1)) } 
aligned <- function(x,y,e)
{
  rl <- lm(y~x);
  m <- as.numeric(rl[[1]][2]);
  q <- as.numeric(rl[[1]][1]);
  d <- distance(x,y,m,q)<e;
  tmp <- TRUE;
  for(i in d) tmp <- tmp & i;
  tmp
}

db<-read.csv(file="C:\\coordinates.csv",sep=";",header=T) 
x<-db$x[1:n_plots]
y<-db$y[1:n_plots]

n_plots <- 40 # Limits the points to the first 40
e<-50         # Tolerance of 50 meters

plot(x,y,col="light gray")
a<-list()
v<-vector(mode="list",length=n_plots)

#---------------------------------------------------------------------
# Records alignments with 2 points (ALL!)
#---------------------------------------------------------------------

for(i in 1:n_plots) for(j in 1:n_plots) if(j>i) a<-c(a,list(p=c(i,j)))
v[[2]]<-a
j<-2

#---------------------------------------------------------------------
# Records alignments with j points
#---------------------------------------------------------------------

while(length(v[[j]])>0)
{
	a<-list()
	for(c in v[[j]])
		if(c[j]<n_plots)
			for(i in (c[j]+1):n_plots)
				if(aligned(c(x[c],x[i]),c(y[c],y[i]),e))
					a<-c(a,list(p=c(c,i)))
	v[[(j+1)]]<-a
	j<-j+1
}

#---------------------------------------------------------------------
# Prints the result
#---------------------------------------------------------------------

best <- (j-1)
freq <- c(rep(0,times=best+1))
for(i in 2:best) freq[i] <- length(v[[i]])
for(i in v[[best]])
{
	points(x[i],y[i])
	lr<-lm(y[i]~x[i])
	abline(lr)
}
freq