```{r init, , echo=FALSE}
#
#  Jaromir D.B. Nemec 19.6.2016
#
# library(knitr)
# knit2html("linear_clustering_3d.rmd")
#
time <- format(Sys.time(), "%a %b %d %X %Y")
#
library(lattice)
#
options(dec=".") ## decimal point
options(scipen=999) ## remove scientific format
options(digits=8)
```
# Linear Clustering

Started on `r Sys.time()`.
## The Task
Given data points in a three dimensional plane define a *k* lines so that a maximum number of points are covered with the lines.
A points is covered with a line if the distance of the points from the line is below a defined limit. 

We will use following parameters to define the line:

```{r conf, echo=FALSE}
line_width <- 0.6
min_line_length <- 140
max_line_length <- 175
min_line_points <- 18
```

Parameters:

**line width** - line_width/2 is the allowed distance of the point from the ideal line =   `r line_width`

**minimum line length**  - points with shorter distance are not considered to belong to the same line =  `r min_line_length`

**maximum line length**  - points with longer distance are not considered to belong to the same line =  `r max_line_length`

**minimum points on a line** - lines with less points are ignored = `r min_line_points`

### Sample Data

The sample data is stored in the csv file [linear_clustering_3d.csv](linear_clustering_3d.csv) with the following layout:

X,Y,Z   coordinates of the point
 
```{r csv, echo=FALSE}
dfpts <- read.csv(file="pts.csv")
summary(dfpts)
str(dfpts)
```

To speed up the processing we will use only 20% sample of the data.

```{r csvsmpl, echo=FALSE}
###### test with smaller number of lines
###### dfpts <- head(dfpts,600)

## sample
df <- dfpts[sample(nrow(dfpts), round(nrow(dfpts)*.20)), ]
summary(df)
str(df)
```

## Graph Data

The complete set of the data points is shown below in 3D.

```{r showtable2, echo=FALSE}
# 3D Scatterplot
library(scatterplot3d)
xy <-  scatterplot3d(dfpts$X,dfpts$Y,dfpts$Z, main="Generated Points")
``` 


## Approach
I developed the following algorithm using a dedicated distance measure.

Following steps are performed:

1) Define a **distance metric** returning: 

   *zero* - if the points do not belong to a line
   
   *Euclidian  distance of the points* - if the points constitute a line according to the defined parameters, i.e.
    
    - their distance is higher or equal than the *min_line_length* and
    
    - their distance is lower or equal than the *max_line_length* and
    
    - the line consists or at least *min_line_points* points with a distance lower that *line_width/2* from the line                        
  

2) Calculate distance matrix using this distance measure (use sample of the data for large data sets; adjust the line parameters accordingly)

3) Find the points A and B with maximum distance - break to step 5) if the distance is zero

   Note that if the distance is higher than zero the points A and B are building a line based on our definition 

4) Get all points belonging to the line AB and remove them from the distance matrix. Repeat the step 3) to find another line

5) Check the coverage of the point with the selected lines, if substantial number of points remains uncovered, repeat the whole algorithm with adjusted line parameters. 

6) In case that data sample was used - reassign all points to the lines and recalculate the boundary points.

## Define Distance Measure

Step 1 - Here we define the distance measure

```{r dist1, echo=TRUE}
## distance of point p from line (a,b) in 3D
## see
## http://stackoverflow.com/questions/35194048/using-r-how-to-calculate-the-distance-from-one-point-to-a-line


dist3d <- function(p,a,b) {
  ##print(paste ("in dist3d p ",  (p), "a ",  (a), "b ",  (b)))
  v1 <- a - b
  v2 <- p - a      
  v3 <- cross3d_prod(v1,v2)
  area <- sqrt(sum(v3*v3))/2
  d <- 2*area/sqrt(sum(v1*v1))
}

cross3d_prod <- function(v1,v2){
  v3 <- vector()
  v3[1] <- v1[2]*v2[3]-v1[3]*v2[2]
  v3[2] <- v1[3]*v2[1]-v1[1]*v2[3]
  v3[3] <- v1[1]*v2[2]-v1[2]*v2[1]
  return(v3)
}

euc.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2))
 

## Distance of two points a,b in data.frame df using line parameters
## returns 0 if the distance of a and b is below the min_line_length
## returns 0 if the number of points with the distance from the line a,b less than line_width/2 is less than min_line_points
## otherwise returns Euclidian distance of the points a and b


dist_line <- function(a, b, pts = dfm, l_w = line_width, m_l_l = min_line_length, m_l_p = min_line_points ) {
  ##print(paste ("in function a ", str(a), "b ", str(b)))
  ed <- euc.dist(a,b)
  if (ed < m_l_l) return(0)
  
  ## count points
  dd <- apply(pts,1,function(x) dist3d(x, a, b))
  ptcnt <- sum(dd < l_w / 2)
   
  if (ptcnt <  m_l_p)  
    {dist <- 0}
  else
    {dist <- ed}  
  return(dist)
}

## Line Points Count 
## used to report the number of points associated with a line 
line_pts_cnt <- function(a, b, pts = dfm, l_w = line_width, m_l_l = min_line_length, m_l_p = min_line_points ) {
  ##print(paste ("in function a ", str(a), "b ", str(b)))
  ed <- euc.dist(a,b)
  if (ed < m_l_l) return(0)
  
  ## count points
  dd <- apply(pts,1,function(x) dist3d(x, a, b))
  ptcnt <- sum(dd < l_w / 2)
   
  if (ptcnt <  m_l_p)  
    {dist <- 0}
  else
    {dist <- ed}  
  return(ptcnt)
}
```

## Calculate Distance Matrix

The following **r** code implements the described algorithm - steps 2 - 4 (Step 5 is not covered).
```{r clust1, echo=TRUE}
## convert the data frame to matrix and remove the line_id
dfm <- as.matrix(df[,c("X","Y","Z")])


## Step 2 - Calculate distance matrix using the distance measure
require(proxy)
d <- dist(dfm , dist_line)

dm <- as.matrix(d) 

## reset upper triangle
dm[upper.tri(dm)] <- 0

lines <- list()
i <- 1

repeat {
 ## Step 3 - Find the points A and B with maximum distance

 mx <- which(dm == max(dm), arr.ind = TRUE)

 if (dm[mx][1] == 0) {break}

 a <- as.numeric(df[mx[1,1],])
 b <- as.numeric(df[mx[1,2],])
 ## 
 
 print ( paste ("found line A=(",paste(round(a,2),collapse=','), ") B=(", paste(round(b,2),collapse=','),") Point Count= ",
   line_pts_cnt(a,b),
   sep = ''))
 lines[[i]]   <- list(A=a,B=b)
 i <- i + 1
 ### get all points on the line ab

 pts_dist <- apply(df,1,function(x) dist3d(x, a, b))
 
 pts_on_line <- df[pts_dist < line_width / 2,]

 ## eliminate from the distance matrix all columns of points on the line
     
 dm[,rownames(pts_on_line)] <- 0
} 

```

## Show Result

The assigned lines are shown on the diagram below.

```{r showtable3, echo=FALSE}
# found lines
##lines

# 3D Scatterplot
library(scatterplot3d)

xy <-  scatterplot3d(dfpts$X,dfpts$Y,dfpts$Z, main="Matched Lines")

for (i in (1 : length(lines)))  {
 ## add lines
 p2 <- xy$xyz.convert(lines[[i]]$A[1],lines[[i]]$A[2],lines[[i]]$A[3])
 p3 <- xy$xyz.convert(lines[[i]]$B[1],lines[[i]]$B[2],lines[[i]]$B[3])
 segments(p2$x,p2$y,p3$x,p3$y,lwd=2,col=2)
 } 
```

Note that the lines are matched based on sample data, so they do not cover teh full length.
An additional step considering all points that recalculate the boundary points of the line segemnts should be performed.
Step 6 is not shown.

Processing finished on `r Sys.time()` using knitr source file [linear_clustering_3d.rmd](linear_clustering_3d.rmd).




